临界梯度模型的优化及集成模拟中高能量粒子模块的搭建*

邹云鹏 陈锡熊 陈伟

1) (核工业西南物理研究院,成都 610041)

2) (通用原子能,加州 CA 92186-5608)

3) (中国科学技术大学核科学与技术学院,合肥 230026)

基于临界梯度模型,使用TGLFEP 和EPtran 两个程序可以模拟分析阿尔芬本征模引起的高能量粒子输运问题.本文在原有模型的基础上,加入了两点改进使模拟结果更接近实验.其一,考虑阈值剖面的演化过程.判断阈值的物理量由密度梯度(dn/dr)改为归一化的密度梯度((dn/dr)/(n/a)),并且使用TGLFEP 模拟证明阈值与高能量粒子密度成反比例关系,也就是说,当密度降低时,阈值会增大.第二,考虑有限轨道宽度效应.使用OBRIT 程序计算高能量粒子的损失锥,并输入到EPtran 程序中,从而增加了一种高能量粒子的损失通道.利用DIII-D#142111 和#153071 进行实验验证,结果表示改进后的模型更接近实验.最后,利用神经网络代替TGLFEP 计算临界梯度,并实现EPtran 的并行计算以缩短运行时间.以此建立一个EP 模块并加入到OMFIT 集成模拟中,模拟结果表示当阿尔芬本征模驱动高能量粒子输运,会导致芯部的压强和电流降低,从而提升当地的安全因子,改变平衡位形.

1 引言

中国聚变工程实验堆(CFETR)[1,2]的主要目的之一是实现自持燃烧,因此高能量粒子(EP)分布的预测成为了一个重要课题.在过去的研究中,通常使用经典慢化分布(classical slowing down)[3,4]模型,但是DIII-D 的实验表明在不稳定的阿尔芬本征模(AE)存在时,诊断得到的EP 剖面与经典慢化分布预测的剖面相差甚远[5,6].类似的现象在其他实验装置上[7-9]也同样被观测到.在理论上,EP 在实空间的梯度会驱动AE,而不稳定的AE会反过来增强EP 的径向输运,从而导致EP 密度剖面变平.当考虑多个不稳定的AE 时,在实空间会产生岛状结构的重叠,这会导致粒子的随机扩散[10],当然这种现象需要AE 的振幅足够大时才会发生,因此本文重点讨论的是单模影响EP 分布的情况.此外,AE 还会改变EP 的运动轨迹从而增加轨道损失,当AE 驱动EP 向外输运时会将粒子推进损失锥中,从而增大EP 的损失份额.因此,在实际托卡马克中,EP 分布和AE 稳定性构成了一个复杂的负反馈系统.

为了解释经典慢化模型与实验结果的区别,并使模拟结果更接近实验,国际上研发(或改进)了,例如ORBIT[11],MEGA[12]和kick model[13]等模型.然而自洽的数值模拟需要大量的计算资源,因此本文选择一个简化的模型来预测EP 输运和再分布,即临界梯度模型.该模型的主要思路为AE的不稳定性主要由EP 压强梯度驱动,那么EP 的压强梯度应该存在一个阈值使AE 达到临界稳定(这个阈值的存在已经在理论[14]和实验中被证实[15]),当EP 压强梯度超过这个阈值时,AE 的振幅快速上升,EP 输运系数快速增大(stiff transport),直到EP 压强梯度降回到阈值.基于这个理论,本文使用TGLFEP 和EPtran 的组合来计算EP 的输运和再分布,其中TGLFEP 用来计算临界梯度剖面,而EPtran 用来计算EP 的再分布.并在此基础上,增加两个改进点: 第一,优化临界阈值的计算方法.EP 梯度的阈值不再是一个固定值,而是与EP 密度成反比例关系,当EP 密度降低时,阈值会随之增大;第二,考虑有限轨道宽度效应(FOW)造成的EP 损失.旧模型仅考虑EP 在最外层闭合磁面(LCFS)处,由于扩散而造成的损失,本文使用ORBIT 程序计算EP 的运动轨迹并确定其相空间的损失锥.

采用两次安全因子剖面不同DIII-D 放电实验来验证改进后的临界梯度模型,其中#142111 为反磁剪切[16],#153071 为单调递增[17].EP 的初始分布使用TRANSP[18]无输运的模型计算得到.使用临界梯度模型模拟环向阿尔芬本征模(toroidal Alfvén eigenmode,TAE)引起的EP 输运,并将输运后的损失份额和沉积剖面与实验进行对比.最后,利用神经网络模型替代TGLFEP 计算临界梯度,并配合EPtran 程序的并行化,缩短计算时间.基于OMFIT 集成模拟框架,搭建一个高能量粒子模块(module of EP,MOE),从而实现集成模拟中高能量粒子压强自洽的迭代.

2 临界梯度模型

2.1 临界梯度计算

TGLF 程序基于回旋朗道流体模型,主要用来计算湍流物理,但是程序本身采用局域的流管(flux-tube)模型,因此同样可以用于AE 的计算.由于TGLF 使用准线性模型,因此其计算速度要远大于回旋动理学程序.基于DIII-D 的H 模放电实验,TGLF 程序已经用GYRO 程序线性和非线性模拟校准.TGLF 程序主要求解6 个矩方程,即密度、平行速度、平行压强、总压强、平行能量流和总能量流.本征方程采用系列厄密函数表示.过去研究结果表示,TGLF 的模拟结果对基函数中的高斯宽度θw的设置非常敏感[19,20],通常选择对应增长率最大的θw.因为程序自动寻找该值的方法仅对湍流适用,并不适用于AE 的计算,所以改为扫描一定区间内的增长率并选择对应最大值的θw.以#142111 为例,扫描高斯宽度0.4—2.0 区间.另外还有一个原则是模的频率和增长率在一定区间内有良好的鲁棒性.图1 展示了3 个环向模数n=3的模的增长率和频率随高斯宽度的变化.因为橙色和紫色曲线均在一定区间内表现出良好的-鲁棒性,所以被认为是TAE.而蓝色曲线的波动过大,且不满足输出边界条件(outgoing boundary condition)[21],因此被认为是一种数值噪声.TGLF EP 是TGLF 的并行的封装程序,可以自动求解多个磁面的高斯宽度并计算AE 的频率和增长率.

图1 (a)增长率和(b)频率与高斯宽度的关系,其中橙色和紫色的曲线表示TAE,蓝色线被认为是一种数值噪声,作为参考,黑色实线为MEGA 的模拟结果Fig.1.(a) Growth rate and (b) frequency as a function of Gaussian width.The orange and purple curves represent acceptable TAE calculated using TGLFEP.Because the blue curve has a strong perturbation,the mode is treated as polluted by numerical noise.For comparison,the MEGA results are depicted by the black lines.

2.2 高能量粒子输运

在临界梯度剖面确定下来之后,就可以使用EPtran 程序来计算EP 的再分布.EP 在(r,E,λ,σ)空间的输运方程如下:

空间坐标r,E,λ,σ 分别为径向位置、能量、投掷角和运动方向(沿着磁场方向为正,逆着磁场方向为负).方程(1)等号右边包括慢化项、源项和螺旋角空间的散射项[22].其中fEP表示高能量粒子分布函数,νd表示螺旋角散射率,S0表示EP 源,E0为初始能量(α 粒子为3.5 MeV,中性束则为注入能量),τs为慢化时间,Ec为临界能量,TEP为等效麦克斯韦温度,a是小半径.V '=∂rV,V E==v///v,λ=μB/E,这里v表示粒子的速度,v//表示平行磁场方向的速度,μ为磁矩,B为磁感应强度.在外中平面,0 ≤ λ ＜ λTP表示通行粒子,λTP≤ λ ≤ 1 表示捕获粒子,其中λTP(r)=B(0)/B(π)[23].

其中,当x＜ 0 时,[x]＞0=0,为局域的临界梯度;nEP为高能量粒子密度;DITG/TEM为湍流输运系数;DAE并不是从物理推导的,通常将其设置得足够大以驱动EP 快速的输运.在过去的研究中,ITER 模拟设置为0.3[24],DIII-D 模拟设置为10[25].实际上,再分布剖面对DAE的设置并不敏感.因此根据已经发表的文献 [24,25],在计算DIII-D 的实验时,设置DAE=10,在计算CFETR 时,设置DAE=0.3.

图2 能量相关(橙)/不相关(蓝) Angioni 模型对比Fig.2.Comparison between energy dependent (orange) and energy independent (blue) Angioni diffusion coefficients.

根据AE 增长率表达式[27],方程(4)改写为

其中,LnEP为EP 密度的特征长度,将在3.1 节详细说明此改动.

利用方程(1)计算EP 的分布后,对其积分可以得到EP 的密度和等效麦克斯韦温度

EP 的初始分布采用经典慢化分布:

根据方程(6)可知,EP 密度剖面为nsd(r)=S0τsI2(vc/v0),温度为Tsd=(2I4/3I2)E0.其中积分函数被定义为.对于α 粒子,S0由背景DT 密度和相应的反应截面计算;对中性束注入来说,S0由NUBEAM 程序[28]计算.

3 临界梯度模型的优化

3.1 阈值剖面演化

在已经发表的研究中,Waltz等[23-25]讨论了多种临界稳定的判据,本文使用判据γAE+ITG/TEM=γITG/TEM来判断是否达到临界稳定(γAE+ITG/TEM指考虑背景等离子体时AE 的增长率,γITG/TEM指湍流的增长率).在旧模型中,EP密度的特征长度(LnEP)保持不变,不断降低EP 密度(密度梯度成比例降低),直到达到临界稳定的判据,此时的密度梯度则被视为临界梯度.在这个假设下,密度剖面的形状是不变的,显然这不符合实际情况.因此使用TGLFEP 计算一个二维矩阵(a/LnEP,nEP)来确定临界的归一化梯度与EP 密度的对应关系,利用这个矩阵,可以在EP 密度下降时,使a/LnEP上升.为了说明我们的改进,选择一个解析的AE 增长率的表达式[27]:

其中,ωi为增长率,k//为平行波矢,vA为阿尔芬速度,n为环向模数,βc和βα分别为背景粒子和高能量粒子的比压.G和H是与粒子速度相关的拟合函数,具体表达式见文献 [27].等式右边的第1项为背景离子和电子的朗道阻尼,当假设平衡不变时,其为常数.第2 项包括了EP 对AE 的阻尼和驱动,其中βα和δα(∝dpα/(pαdr),这里pα为高能量粒子压强)两个量会在EP 输运时改变.由于EP 的驱动中同时包含压强和压强梯度,因此应随着EP 剖面变化而变化.令方程(8)等于0,则可以得到如下形式:

其中k1包含了背景粒子的朗道阻尼,k2包含了EP 本身的阻尼和温度梯度驱动,在忽略温度变化时,这两个系数在每个磁面上都是固定的.到此为止,我们证明了与nEP成反比例关系.图3(a)中,红色虚线表示对应初始EP 密度的a/LnEP,黑色实线表示经典慢化分布的a/LnEP,黑色实线高于红色虚线的区域表示AE 是不稳定的.针对这个不稳定区域,使用TGLFEP 扫描a/LnEP与nEP,图3(b)展示了3 个不同位置,TGLFEP 的模拟结果也验证了这种反比例关系.当EP 输运时,芯部的EP 密度是降低的,那么阈值则会相应的上升.以ρ=0.4 为例(图3(c)),红色的叉表示初始时刻,黑色和红色箭头分别表示旧模型和改进后的模型下随EP 密度的变化轨迹,这可以明显看到两个模型的差别.

图3 (a)输运前的AE 不稳定区间;(b)在ρ=0.4,0.5,0.6 处,a/LnEP与nEP 的反比例关系;(c) ρ=0.4 处,旧模型和改进模型模拟的密度演化轨迹Fig.3.(a) Unstable AE region of n=3 without transport.Critical a/LnEP is depicted by dash red curve,and a/LnEP of classical slowing down distribution is depicted by solid black curve;(b) inverse proportional function between critical a/LnEP and nEP at ρ=0.4,0.5,0.6;(c) density evolution trajectories of previous CGM (black) and improved CGM (red) at ρ=0.4.

3.2 有限轨道宽度(FOW)效应

磁场的不均匀性导致粒子导心轨道的漂移被称为有限轨道宽度(FOW)效应,在等离子体边界附近,能量越高的粒子就越容易在FOW 效应的作用下从LCFS 损失掉,这种现象在实验中是普遍存在的[16].然而旧模型并没有考虑这种损失机制,因此只有很少的EP 会损失,而大量的粒子都沉积在等离子体边界附近.本文使用ORBIT 程序计算EP在相空间中的损失锥,并将其加入到输运模型中.

在ORBIT 计算中,初始EP 分布选择在(ρ,E,ξ)空间的均匀分布,且仅考虑外中平面的粒子.由于ORBIT 程序采用二维空间坐标,因此这里仅考虑外中平面.当然,由于粒子更容易在低场侧损失,因此这会导致过度的估计EP 的损失份额.此外,在AE 扰动下,一部分原本可以约束的粒子会转变为损失粒子,即AE 会增大EP 的损失锥.以#142111 为例,图4(a)和图4(b)分别展示了捕获粒子和通行粒子在无扰动(蓝色)和有AE 扰动(红色)下的运动轨迹.其中通行粒子在高场侧的轨迹被放大到图4(c)中.在每幅图中,黑色叉表示粒子的初始位置,红色叉表示粒子从LCFS 逃逸的位置,黑色箭头表示运动方向.

图4 (a)捕获粒子和(b)通行粒子在AE 扰动下的运动轨迹示例;(c)图(b)在高场侧的放大图Fig.4.Representative trajectories of (a) trapped and (b) counter passing particles by including AE perturbation;the high magnetic field side of (b) is enlarged in (c) to reveal trajectory variations.

图5(a)展示了无AE 扰动时不同能量的EP的损失边界,在该边界的右侧表示粒子的损失区域.由于初始粒子均选择在外中平面上,因此所有的损失粒子都具有负的螺旋角.损失锥的面积随粒子能量增大而增大.对于能量为注入能量80 keV的粒子,其损失锥的峰靠近磁轴.图5(b)展示了AE 扰动下损失锥的进一步增大.根据文献 [5,16],AE 的扰动幅度设为dB/B≈O(10-4).将得到的损失锥输入到EPtran 程序中,令损失锥内的EP 分布函数为零,从而等效的考虑EP 由于FOW 效应的损失.

图5 (a)无扰动时,不同能量EP 在(ρ,v///v)空间的损失锥;(b)有无AE 扰动下,80 keV 的EP 损失锥对比图Fig.5.(a) Loss boundary in (ρ,v///v) space with different energies;the loss region is on the right side of the curve.(b) loss cone without/with AE perturbation.The blue area represents the loss cone without AE,and the additional loss by including AE is highlighted by the red area.

4 临界梯度模型的实验验证

4.1 #142111 放电实验

实验中利用ECE 测量到525 ms 时刻频率为70—90 kHz 的低nTAE,这些模在空间上重叠,径向位置均在ρ ≈ 0.4.MEGA 程序的模拟表示,n=1—5 的TAE 的能量演化如图6(a)所示,n=3的TAE 增长率最大且最先达到饱和.该TAE 在线性阶段((a)中红色虚线对应的位置)的模结构如图6(b)所示.图7(a)展示了TGLFEP 程序计算的n=1—5 的模的增长率在径向上的分布情况,可以看到从ρ ≈ 0.1 到0.6 这个范围内,n=1—4 分别占主导,因此不同径向位置的临界梯度均考虑当地增长率最大的模(图7(b)).

图6 MEGA 模拟结果 (a) n=1—5 的TAEs 的能量演化过程;(b) n=3 的TAE 的模结构,vr,cos 为径向扰动速度的cos 分量,不同的极向谐波用不同的颜色区分Fig.6.(a)Evolutionof energy with n=1-5 TAEs by MEGA;(b) cosine part of radial velocity for the most unstable n=3 TAE.

图7 (a)n=1—5 的TAEs在不同磁面的增长率;(b)单n和多n计算的临界梯度剖面的对比Fig.7.(a) Growt hratein eachfluxsurfaceofn=1-5TAEsbyTGLFEP;(b)comparisonbetweencritical a/LnEP profiles with n=1-4 and n=3.

输运后的EP 径向密度剖面如图8(a)中的紫色曲线所示,为了可以更清晰比较,图8 中还加入了经典慢化分布,即输运前的剖面(黑色)、实验结果(红色)、旧临界梯度模型的模拟结果(蓝色)、改进后的临界梯度模型但不考虑AE 对损失锥影响的模拟结果(绿色)以及MEGA 的结果[29].从图8(a)可以看出,改进后的临界梯度模型的模拟结果比旧模型更接近实验结果.这里定义EP 的损失份额为损失的粒子数与经典慢化分布预测的约束的总粒子数之比,那么旧模型的损失份额为15%,而改进后模型的模拟结果超过40%,这与实验的约50%的损失份额非常接近.而且AE 区域梯度过小和EP 堆积在边界附近的问题也都有所改善.在临界梯度模型中,EP 的输运过程可以理解为: 不稳定的AE 不仅增强EP 的径向输运,而且会轻微增大损失锥,EP 向外输运进入损失锥中后,由于FOW 效应而直接从LCFS 逃逸出去.最后,由于考虑了负螺旋角的分布,模型可以计算EP 在螺旋角空间的分布(图8(b)),黑色曲线对应60 keV 的EP 的分布,红色曲线则为0—80 keV全部EP 的分布.

图8 (a) EP 密度剖面对比图,经典慢化分布(黑色)、实验(红色)、旧临界梯度模型(蓝色)、考虑(紫色)/不考虑(绿色)AE 扰动对损失锥影响的改进的临界梯度模型、MEGA(黄色);(b) EP 在螺旋角空间的分布Fig.8.(a) Density profile comparison: Black curve represents classical slowing down;the red curve is inferred from experiment data;blue curve represents original CGM without loss cone effect;purple/green curve is improved CGM with loss cone from AE perturbed/unperturbed orbits;yellow curve represents MEGA results.(b) EP redistribution in pitch angle space.

4.2 #153071 放电实验

在#153071 放电实验中,频率在100—200 kHz的低nTAE 在靠外侧的位置被激发,但是忽略等离子体旋转,模拟得到了一个在ρ=0.2—0.6 范围内的一个较宽的模.增长率最大的n=4 的TAE如图9 所示,图中实线和虚线分别为扰动速度的余弦和正弦分量.FOW 效应引起的损失锥如图10所示,由于本次实验AE 振幅较弱,且AE 的径向宽度与损失锥的重叠较小,因此在计算损失锥时没有考虑AE 扰动的影响.

图9 n=4 的TAE 的模结构,实线和虚线分别表示扰动速度的余弦和正弦分量Fig.9.Spatial profile of n=4 TAE,where cosine and sine part of radial velocity are depicted by solid and dash curve,respectively.

图10 不同能量的EP 在(ρ,v///v)空间的损失锥Fig.10.Loss boundary in (ρ,v///v) space with different energies.

将模拟得到的阈值剖面和损失锥代入到EPtran 程序中后,计算得到的EP 再分布剖面如图11 中的红线所示,同样,为了便于比较,图中还给出了经典慢化分布(黑色)、实验分布(绿色)及旧模型的模拟结果(蓝色).改进后的模拟结果呈类似高斯分布的形式,与实验结果非常吻合,模拟曲线几乎都在误差棒内.

图11 压强剖面对比图,经典慢化分布为黑色曲线,旧临界梯度模型为蓝色曲线,改进的临界梯度模型为红色曲线,实验结果为绿色三角并配有误差棒Fig.11.Pressure profile of classical slowing down (black),previous (blue) and improved (red) CGM.For comparison,experimental data is depicted by green triangles with error bar.

4.3 基于神经网络构建HL-3 的高能量粒子模拟模块

目前OMFIT 集成模拟中,EP 的计算主要是用ONETWO (α 粒子)和NUBEAM (NBI)程序实现.但这两个程序使用的是经典慢化模型,模拟结果与实际分布相差很大.因此本章重点讨论如何将临界梯度模型嵌入集成模拟中.根据第2 节的模型,可知利用临界梯度模型预测EP 分布分为两步,第1 步是利用TGLFEP 计算临界梯度,第2步是利用EPtran 程序计算EP 的再分布.目前EPtran 程序已经改写为并行程序,～5 min 就可以完成计算,因此,只有缩短第1 步的运行时间,才有可能将临界梯度模型加入到集成模拟的迭代中,这里就选用神经网络的方法来实现.

搭建神经网络第1 步需要合适的数据库,样本的输入太多会导致模型预测不准确,因此从物理出发,尽量减少输入变量.这里不考虑杂质的影响,因此有效电荷数Zeff=1,而EP 的密度是扫描量,因此离子的密度可以通过准中性条件来计算.为了方便,认为离子温度与电子温度相同.参考TGLFEP的输入后,共有18 个输入变量 (表1).模型的输入则为3.1 节中方程(9)中的两个系数k1和k2.根据输入输出,参照中国环流3 号(HL-3)的参数范围,利用EFIT 生成150 个不同的平衡,每个平衡随机选取8 个磁面的参数作为样本,因此总共生成1200 个样本.

表1 神经网络输入变量Table 1.Variables for NN input.

神经网络模型采用最简单的线性结构,主要有两个功能: 1)预测是否存在不稳定的AE;2)对于存在不稳定AE 的磁面,计算k1和k2.流程图如图12 所示.

图12 神经网络流程图,第1 步判断AE 是否被激发,第2步计算系数k1和k2Fig.12.Flow chart of the neutral network (NN).The NN estimates if the AE can be excited at first,and calculates the two coefficients k1 and k2 for AE unstable location.

将1200 个样本进行随机洗牌,80%用来训练模型,剩下的20%用来对模型进行检验.第一阶段是判断是否存在不稳定的AE,图13 表示在模型学习200 次之后,均方差大约收敛为0.1,准确率收敛到0.9 左右.继续训练会发生过拟合.对于存在AE 的样本,同样选择80%作为训练集,20%为验证集.图14 展示了预测值与目标值之间的差距((a)k1和(b)k2),圆点越接近对角线(虚线)表示预测值越接近目标值.从图14 可以看到,迭代次数越高,预测值就越接近,但总体变化不大.

图13 (a)训练集和验证集的均方差随学习次数的变化;(b)准确度随学习次数的变化Fig.13.(a) Loss (mean square error) and (b) accuracy for AE stability estimation versus training epoch.

图14 预测值与目标值的对比图 (a) k1;(b) k2Fig.14.Predicted (a) k1 and (b) k2 compared with the targets.

此外,还额外制作了3 个HL-3 的平衡来验证这个神经网络模型,3 个平衡具有不同的安全因子剖面,分别为单调递增、芯部弱反磁剪切、强反磁剪切,这3 个平衡的压强和安全因子剖面如图15所示.此外,还根据高能量粒子临界梯度与高能量粒子密度成反比这一物理,采用了一个自定义的损失函数作为对比.图16为k1和k2的预测值和目标值(TGLFEP 计算值)的径向分布,可以看到预测值与目标值相比误差不大,判断的AE 的范围基本相同.将二者代入EPtran 中计算的高能量粒子剖面如图17 所示,这说明训练的神经网络可以替代TGLFEP.

图15 用来验证神经网络的3 个平衡的(a)压强剖面和(b)安全因子剖面Fig.15.(a) Pressure and (b) safety factor profile of three additional equilibria for NN validation.

图16 HL-3 神经网络预测值(MSE 为蓝色,自定义损失函数为红色)与TGLFEP 计算值(黑色)对比图 (a),(b) Case 1;(c),(d) Case 2;(e),(f) Case 3Fig.16.Coefficients of k1 and k2 predicted by NN with loss function of MSE (blue) and custom loss function (red): (a),(b) Case 1;(c),(d) Case 2;(e),(f) Case 3.For comparison,TGLFEP results are depicted by black curve.

图17 具有(a)单调递增、(b)芯部弱反磁剪切和(c)强反磁剪切安全因子剖面的平衡位形下,高能量粒子剖面对比图.绿色曲线为经典慢化模型计算的初始EP 剖面,黑色、蓝色、红色曲线分别为根据TGLFEP、MSE 和自定义损失函数得到的临界梯度计算的EP 剖面Fig.17.EP profile comparison for the equilibrium with (a) monotonic,(b) weak and (c) strong shear q-profile.In each panel,green curve depicts initial EP profile with classical slowing down distribution and black curve depicts EP profile with the critical gradient calculated by TGLFEP.The blue and red curves depict EP profiles by NN with loss function of MSE and custom loss function,respectively.

最后采用神经网络模型代替TGLFEP 程序,与并行的EPtran 程序组合,构建MOE.为了提高运行速度,MOE 中使用解析公式[30]代替ORBIT程序计算EP 的损失锥.图18 为OMFIT 工作流,右侧的TGYRO,ONETWO 和EFIT 三个程序构成了目前常用的集成模拟(不考虑台基区).TGYRO计算各粒子的密度、温度剖面;ONETWO 计算电流和粒子源;EFIT 构建等离子体平衡.我们搭建的MOE 则读取TGYRO 计算的粒子剖面信息、EFIT 生成的平衡的磁场信息,以及ONETWO 计算的中性束粒子源(α 粒子源则根据DT 的剖面进行计算),计算得到EP 的压强和电流剖面传递给ONETWO,然后计算得到总压强和总电流(P'为压强的导数,即压强梯度,F是一个与极向电流相关的函数,F'为F的导数)再传递给EFIT 完成迭代.

图18 OMFIT 集成模拟流程图,蓝色框里是旧OMFIT 迭代流程,红色框中的是高能量粒子模块Fig.18.OMFIT workflow with MOE.The typical iteration is in the blue border,and MOE is in the red border.

利用带有MOE 的集成模拟计算的等离子体平衡如图19 所示,其中绿线仅考虑EP 输运对总压强的影响,红线考虑EP 输运对压强和电流的影响,蓝线则为不包含MOE 的集成模拟的结果.可以看到EP 输运导致芯部的总压强和总电流下降,从而使得芯部的安全因子增大,在一定程度上有利于反磁剪切位形的形成.

图19 利用带有MOE 的集成模拟计算的(a)磁面、(b)总压强、(c)安全因子剖面,其中绿线仅考虑了EP 输运对总压强的影响,红线考虑了EP 输运对压强和电流的影响,蓝线则为不包含MOE 的集成模拟的结果Fig.19.(a) Flux surface,(b) total pressure and (c) safety factor profile calculated by the OMFIT integrated simulation with MOE.Green curve only considers pressure modification,red curve considers both pressure and current modification,and blue curve is calculated without MOE.

5 结论

本文介绍了临界梯度模型以及使用TGLFEP和EPtran 的组合模拟AE 导致EP 再分布的方法.两点改进使模拟结果更接近实验结果: 第一,根据解析理论,将判断临界稳定的物理量从密度梯度改为归一化的密度梯度,并修改相应的扩散系数表达式.同时优化阈值扫描的数值方法,令临界的归一化密度梯度成为EP 密度的反比例函数,系数由TGLFEP 计算得到,从而使阈值可以随着EP的变化而变化.第二,使用ORBIT 计算EP 在相空间的损失锥,并输入到EPtran 程序中,从而增加一种由FOW 效应导致的EP 损失通道.改进后的临界梯度模型使用#142111 和#153071 两次实验进行验证,模拟得到的EP 剖面与MEGA 结果类似,与旧模型相比更接近实验结果.EP 的损失率也达到42%,这与实验观测的结果很接近.

在此基础上,还利用神经网络代替TGLFEP,实现了EPtran 程序的并行计算,加快了整个流程的计算速度,并为HL-3 的集成模拟搭建了一个高能量粒子模块(MOE).集成模拟的结果表示当AE驱动EP 输运,使得芯部压强和电流减小,从而影响整个平衡.