本发明属于磁约束受控核聚变领域，涉及磁约束受控核聚变领域托卡马克装置放电的数值模拟，特别涉及一种用于分析托卡马克中高能量粒子约束性能的测试粒子模拟方法。
背景技术：
随着人类社会的高速发展，能源问题成为全人类面临的重要难题。传统化石能源日益枯竭，在使用过程中会产生大量有毒有害物质；而风能、水能、太阳能等新型能源由于其对地理环境等因素的苛刻要求，并不能很好地替代化石能源。相比之下，受控核聚变能因其原料来源广泛、环境友好、固有安全性高等诸多优势被人们认为是人类理想的新能源。托卡马克装置被认为是最有希望解决人类能源问题的磁约束聚变装置，其内部由于发生氘氚聚变或因辅助加热(如中性束注入、离子回旋共振加热)会产生大量的高能量粒子，这些高能量粒子可以加热等离子体，促进聚变反应的持续进行。但是，如果这些高能量粒子没有被很好地约束，它们将会打到托卡马克器壁上，对装置造成严重的损害。在未来聚变反应堆中(如中国聚变工程实验堆cfetr)，高能量粒子将占有很大的份额，实验操作不当会造成巨大的经济损失，因此在设计阶段对托卡马克内部高能量粒子的约束性能进行数值模拟显得尤为重要。
对高能量粒子的数值研究有多种方法，而测试粒子方法因其模型相对简单、物理图像清晰，在高能量粒子的模拟研究中扮演重要的角色。本发明的粒子推动部分采用导心轨道方程与全轨道方程，通过相空间坐标的相互转化，实现两种轨道的混合计算。为更加准确完整地描述高能量粒子在托卡马克装置内部的运动特征，本模型同时包含芯部主等离子体区与边界等离子体区。本发明可以载入实验和解析模拟得到的多种扰动电磁场，并实现多场耦合；通过划分三角形网格使电磁场的空间分布更易被表达为任意函数，并实现由磁面坐标到几何坐标的转化。本发明可以结合实际实验装置，快速准确完整地描述高能量粒子的输运与损失，计算效率高，数值稳定性强，是一种精确高速的数值模拟方法。
技术实现要素：
为了实现上述目的，本发明提供一种能够结合实际实验装置快速分析高能量粒子在不同三维扰动场存在条件下的测试粒子模拟方法。
本发明的技术方案：
一种用于分析托卡马克中高能量粒子约束性能的测试粒子模拟方法，包括以下步骤：
首先，对平衡拟合程序(efit)生成的平衡文件gfile进行重构，匹配生成非结构化三角形网格，建立三维柱坐标系(r,ζ,z)，得到平衡磁场。其中，r、ζ和z分别是径向、环向和轴向坐标。然后，根据所考虑的高能量粒子来源，载入粒子初始相空间分布信息，并对粒子进行编号。载入初始磁场b(r,ζ,z,t0)与初始电场e(r,ζ,z,t0)信息，其中，t0为初始时刻。接着，根据粒子的全轨道方程与导心轨道方程，考虑高能量粒子与背景等离子体之间的库仑碰撞，计算每一时刻粒子信息并对其进行诊断；若粒子超出最外磁面，则判定该粒子为损失粒子，记录其编号、位置信息、速度信息、损失时刻等信息，该粒子直接进入后处理模块；若粒子没有成为损失粒子，则继续推动粒子进入下一时刻t δt。对于随时间演化的扰动电磁场，在t δt时刻重新加载扰动电磁场，在更新后的电磁场中，重复上述步骤直至时间t达到设定的最大计算时间t，输出约束粒子与损失粒子的全部信息。最后，经过后处理过程，输出粒子分布、粒子能量、壁面热负荷、电流等信息。
本发明的有益效果：本发明提供一种能够结合实际实验装置快速分析高能量粒子在不同三维扰动场存在条件下的测试粒子模拟方法。本发明所述的粒子推动部分采用导心轨道方程与全轨道方程，通过相空间坐标的相互转化，实现两种轨道的混合计算；为更加准确完整地描述高能量粒子在托卡马克装置内部的运动特征，本模型同时包含芯部主等离子体区与边界等离子体区。本发明可以载入实验和解析模拟得到的多种扰动电磁场，并实现多场耦合；通过划分三角形网格使电磁场的空间分布更易被表达为任意函数，并实现由磁面坐标到几何坐标的转化。本发明可以结合实际实验装置，快速、准确、完整地描述高能量粒子的输运与损失，能够准确、高效地分析高能量粒子的约束性能，计算效率高，数值稳定性强，是一种精确高速的数值模拟方法。
附图说明
图1为本发明所适用的托卡马克实验装置的三维示意图及高能量粒子在装置内部的运行轨迹图。
图2为本发明用于模拟计算采用的网格二维网格图。
图3为本发明模拟的初始粒子密度、能量、安全因子随归一化小半径的变化关系图。
图4(a)为本发明计算的阿尔法粒子第一轨损失率随时间变化关系图。
图4(b)为本发明计算的损失阿尔法粒子极向分布图。
图5为本发明诊断得到的粒子在相空间分布图像。
图6是本发明的方法流程图。
具体实施方式
下面结合附图和技术方案，进一步说明本发明的具体实施方式。
托卡马克装置是一种环形磁约束装置，其结构如图1所示，高能量粒子在装置内部高速运动。首先，根据真实磁场位形，建立三维柱坐标系(r,ζ,z)，构建计算网格，如图2所示，计算网格既包含芯部主等离子体区域，也包含等离子体与器壁之间的真空区。图3为初始粒子密度、能量、安全因子随归一化小半径的变化关系图，按照实际情况，载入初始粒子分布、初始平衡磁场与扰动电磁场。接下来，根据粒子导心轨道方程、全轨道方程与碰撞方程，计算每一时刻的粒子信息，统计损失粒子信息，如图4(a)、图4(b)所示。最后，完成计算过程，输出粒子全部信息。经后处理模块处理，可以得到粒子在μ-pφ相空间分布图，从而判断粒子属性，如图5所示。
具体的实施步骤如下：
步骤1：对平衡拟合程序(efit)生成的平衡文件gfile进行重构，匹配生成非结构化三角形网格，建立三维柱坐标系(r,ζ,z)，得到平衡磁场。其中，r、ζ和z分别是径向、环向和轴向坐标；
所述的非结构化三角形网格可根据托卡马克装置几何信息进行网格划分。区别于一般的网格划分方法，其在几何信息的基础上，还可以根据磁面的物理信息进行网格划分，具体表现为：网格三角形顶点分布在相互嵌套的平衡磁面上，可利用磁面弧长信息划分网格，保证三角形网格的对称性，保证在极向截面内网格有相同的拓扑结构。这种方法易于划分并行区域，易于根据研究物理问题的实际情况在局部对网格进行加密，提高计算精度。
步骤2：根据所考虑的高能量粒子来源，载入粒子初始相空间分布信息，并对粒子进行编号。相空间分布信息包括每一个粒子的位置信息(r,ζ,z)与速度信息(vr,vζ,vz)，其中，(vr,vζ,vz)为粒子在三维柱坐标系下的速度坐标值。
步骤3：载入初始电磁场信息：
b(r,ζ,z,t0)＝b0(r,ζ,z) b1(r,ζ,z,t0)(1)
e(r,ζ,z,t0)＝e1(r,ζ,z,t0)(2)
其中，b(r,ζ,z,t0)为初始磁场，e(r,ζ,z,t0)为初始电场，b0(r,ζ,z)为平衡磁场，b1(r,ζ,z,t0)为扰动磁场，e1(r,ζ,z,t0)为扰动电场，t0为初始时刻。扰动电磁场数据来自于解析模式或来自于实验数据重构，可以添加纵场波纹(tfripple)扰动、磁流体力学(mhd)扰动、阿尔芬扰动、多模协同扰动等多种扰动模式。
步骤4：采用三维柱坐标系下全轨道方程求解粒子运动轨迹，全轨道方程如下：
其中，(er,eζ,ez)与(br,bζ,bz)分别为电场和磁场在三维柱坐标系下的分量；ze为粒子电荷量；m为粒子质量；e为元电荷。
同时，采用漂移轨道求解粒子的导心运动轨迹，导心轨道方程如下：
其中，x为粒子导心坐标，v||为粒子速度沿磁场方向的分量，μ为粒子磁矩，
其中，b＝b/b，为沿磁场方向的单位矢量。
全轨道方程与导心轨道方程为一阶微分方程，为保证计算精度，采用四阶龙格-库塔法求解，解法如下：
对于一阶常微分方程：
第n 1步与第n步的y值之间的关系为
其中
已知t＝0时刻的y值，通过式(8)-(10)求出y随t的变化关系；
步骤5：高能量粒子与背景等离子体之间的库仑碰撞采用福克-普朗克方程来描述，并通过蒙特卡洛方法进行数值计算，方程如下：
其中，fα为粒子的分布函数，fi为拖拽参数，dik、dii为扩散参数，ξi为满足方差为1的高斯分布随机数，vsd为等效慢化碰撞频率，v||为等效散射频率，δtec为碰撞计算时间步长，n1、n2,3为高斯随机数，v为碰撞后速度，v0为碰撞前速度，
考虑聚变反应，单位时间单位体积的聚变率dr/dv表达式如下：
其中，si分别为第i种粒子的数量，sj为第j种粒子的数量，δij为克罗内克符号，<σvr>表达式如下：
其中，σ为聚变反应截面，m为粒子质量，vr为质心速度，∈为质心能量，kb为玻尔兹曼常数，t为热力学温度。
步骤6：诊断模块对当前每一个粒子的位置信息进行诊断，若粒子超出最外磁面，则判定该粒子为损失粒子，记录其编号、位置信息、速度信息、损失时刻等信息，同时该粒子直接进入步骤8的后处理模块；若粒子没有成为损失粒子，则进入步骤7。
步骤7：对于随时间演化的扰动电磁场，在t δt时刻重新加载扰动电磁场：
b(r,ζ,z,t δt)＝b0(r,ζ,z) b1(r,ζ,z,t δt)(18)
e(r,ζ,z,t δt)＝e1(r,ζ,z,t δt)(19)
在更新后的电磁场中，重复步骤4-步骤7直至时间t达到设定的最大计算时间t，即可输出约束粒子与损失粒子的全部信息。
步骤8：经过后处理过程，输出粒子分布、粒子能量、壁面热负荷、电流等信息。
1.一种用于分析托卡马克中高能量粒子约束性能的测试粒子模拟方法，其特征在于，所述的方法具体包括以下步骤：
步骤1：对平衡拟合程序efit生成的平衡文件gfile进行重构，匹配生成非结构化三角形网格，建立三维柱坐标系(r，ζ，z)，得到平衡磁场；其中，r、ζ和z分别是径向、环向和轴向坐标；
步骤2：根据所考虑的高能量粒子来源，载入粒子初始相空间分布信息，并对粒子进行编号；所述相空间分布信息包括每一个粒子的位置信息(r，ζ，z)与速度信息(vr，vζ，vz)，其中，(vr，vζ，vz)为粒子在三维柱坐标系下的速度坐标值；
步骤3：载入初始电磁场信息：
b(r，ζ，z，t0)＝b0(r，ζ，z) b1(r，ζ，z，t0)(1)
e(r，ζ，z，t0)＝e1(r，ζ，z，t0)(2)
其中，b(r，ζ，z，t0)为初始磁场，e(r，ζ，z，t0)为初始电场，b0(r，ζ，z)为平衡磁场，b1(r，ζ，z，t0)为扰动磁场，e1(r，ζ，z，t0)为扰动电场，t0为初始时刻；
步骤4：采用三维柱坐标系下全轨道方程求解粒子运动轨迹，全轨道方程如下：
其中，(er，eζ，ez)与(br，bζ，bz)分别为电场和磁场在三维柱坐标系下的分量；ze为粒子电荷量；m为粒子质量；e为元电荷；
同时，采用漂移轨道求解粒子的导心运动轨迹，导心轨道方程如下：
其中，x为粒子导心坐标，v||为粒子速度沿磁场方向的分量，μ为粒子磁矩，
其中，b＝b/b，为沿磁场方向的单位矢量；
步骤5：高能量粒子与背景等离子体之间的库仑碰撞采用福克-普朗克方程来描述，并通过蒙特卡洛方法进行数值计算，方程如下：
其中，fα为粒子的分布函数，fi为拖拽参数，dik、dii为扩散参数，ξi为满足方差为1的高斯分布随机数，νsd为等效慢化碰撞频率，v||为等效散射频率，δtec为碰撞计算时间步长，n1、n2，3为高斯随机数，v为碰撞后速度，v0为碰撞前速度，
考虑聚变反应，单位时间单位体积的聚变率dr/dv表达式如下：
其中，si分别为第i种粒子的数量，sj为第j种粒子的数量，δij为克罗内克符号，<σvr>表达式如下：
其中，σ为聚变反应截面，m为粒子质量，vr为质心速度，∈为质心能量，kb为玻尔兹曼常数，t为热力学温度；
步骤6：诊断模块对当前每一个粒子的位置信息进行诊断，若粒子超出最外磁面，则判定该粒子为损失粒子，记录其编号、位置信息、速度信息、损失时刻，同时该粒子直接进入步骤8的后处理模块；若粒子没有成为损失粒子，则进入步骤7；
步骤7：对于随时间演化的扰动电磁场，在t δt时刻重新加载扰动电磁场：
b(r，ζ，z，t δt)＝b0(r，ζ，z) b1(r，ζ，z，t δt)(18)
e(r，ζ，z，t δt)＝e1(r，ζ，z，t δt)(19)
在更新后的电磁场中，重复步骤4-步骤7直至时间t达到设定的最大计算时间t，输出约束粒子与损失粒子的全部信息；
步骤8：经过后处理过程，输出粒子分布、粒子能量、壁面热负荷、电流信息。
2.根据权利要求1所述的方法，其特征在于，步骤1中所述的非结构化三角形网格根据托卡马克装置几何信息进行网格划分，在几何信息的基础上，同时根据磁面的物理信息进行网格划分，具体为：网格三角形顶点分布在相互嵌套的平衡磁面上，利用磁面弧长信息划分网格，保证三角形网格的对称性，保证在极向截面内网格有相同的拓扑结构。
3.根据权利要求1或2所述的方法，其特征在于，步骤4中所述的全轨道方程与导心轨道方程为一阶微分方程，采用四阶龙格-库塔法求解，解法如下：
对于一阶常微分方程：
第n 1步与第n步的y值之间的关系为
其中
已知t＝0时刻的y值，通过式(8)-(10)求出y随t的变化关系。
技术总结