HL-2A 上NBI 加热H 模实验的集成模拟分析*

罗一鸣 王占辉† 陈佳乐 吴雪科 付彩龙 何小雪 刘亮杨曾辰 李永高 高金明 杜华荣 昆仑集成模拟设计组

1) (核工业西南物理研究院,聚变科学技术所,成都 610041)

2) (中国科学院等离子体物理研究所,合肥 230031)

3) (西南交通大学物理科学与技术学院,成都 611756)

托卡马克等离子体物理过程时空尺度跨度大,不同空间区域(如芯部、台基区、刮削层、靶板区)的主要物理过程不同,因此需要采用系统集成方法开展全域多时空尺度物理问题分析.为了更加深入地研究托卡马克等离子体放电实验的稳态运行及爬升期间的输运与约束过程,通常采用多种物理程序开展集成模拟研究,对放电实验结果进行集成模拟对照,相互验证并进一步开展物理分析.本文基于OMFIT 平台,结合HL-2A装置第37012 炮高比压放电实验结果完成了集成模拟验证与分析,验证了程序的可靠性与适用性.在该流程中,通过选取适当的模型,对实验参数进行了校核与补充,经演化后模拟结果与实验结果比较吻合.在此基础上,本文进一步采用TGLF 模型开展了芯部静电漂移波线性不稳定性分析,结果显示NBI 离轴加热导致H 模约束改善的原因是,该实验在NBI 功率沉积位置的ETG 不稳定性处于被抑制的状态,输运由ITG 不稳定性占据主导,同时输运水平降低至新经典水平.

1 引言

目前以托卡马克为代表的可控磁约束核聚变是核聚变研究的一种主要方式.托卡马克中等离子体具有复杂的物理现象,其时间尺度、空间尺度的跨度很大[1],且在不同区域的物理机制也不尽相同,因此无法使用单一理论模型或者物理程序进行准确且快速的模拟.对不同尺度的物理采用不同的模块分别计算,将平衡、输运、不稳定性、加热及电流驱动等程序通过适当的流程设计耦合起来,即可在尽可能快的时间内,得到尽可能精确的计算结果.

目前较为常见的集成模拟平台有OMFIT[2],CRONOS[3],METIS[4]和IMAS[5]等,OMFIT 作为目前国际上最为成熟、最为全面的集成模拟平台,集成了多个等离子物理计算程序,可以涵盖输运[6,7]、电流演化[8]、平衡[9]、不稳定性分析[10]等物理.目前在其他装置上,已采用OMFIT 开展了实验结果的集成模拟分析与验证[11-15].最早在DⅢ-D 的L 模放电[11]中进行了集成模拟分析与验证,而后在H 模放电实验中,实现了芯部与台基的自洽耦合集成模拟,对实验结果实现了很好的模拟验证[12].在EAST 上已经完成了对稳态长脉冲放电[13]、高比压放电[14]和高自举电流放电[15]实验的模拟验证.然而,在国内HL-2A 托卡马克装置上开展的实验结果的集成模拟分析与验证工作还比较少.

此外,集成模拟在新型托卡马克装置、未来聚变堆的放电运行模式设计中也扮演着极为重要的角色.在完成对装置的零维设计后,采用集成模拟对装置的运行模式进行设计,有利于校核设计参数,以及为放电实验提供重要的参考依据.目前已采用OMFIT 的集成模拟流程,完成了对ITER[16]和CFETR[17-19]等装置的运行模式设计与优化.此外,基于METIS和CRONOS 程序,在DEMO上完成了稳态运行模式的设计工作[20].

本文共分为五个部分.第一部分引言介绍OMFIT集成模拟对实验结果验证的研究概况.第二部分介绍本文所采用的OMFIT 集成模拟工作流程及主要程序模块的物理模型.第三部分介绍HL-2A 装置实验放电参数和初始等离子体温度和密度剖面,以及集成模拟计算结果.第四部分进一步分析芯部不稳定性,探讨了芯部约束改善的物理机制.第五部分为本文的总结.

2 OMFIT 集成模拟工作流及物理模型

2.1 OMFIT 集成模拟工作流

本文所采用的集成模拟工作流如图1 所示.该工作流按照任务目标分成两部分:芯部实验结果的集成模拟验证和芯部不稳定性分析.在实验结果验证的集成模拟流程中,包含了平衡程序EFIT[21]、电流演化及源项程序ONETWO[22]和剖面演化程序TGYRO[6,7].ONETWO 同时耦合了辅助加热计算模块,如计算中性束的NUBEAM[23]、计算低杂波的GENRAY[24]、计算ECRH 的TORAY[25]程序模块.本工作流程重点在于芯部实验结果的集成模拟与不稳定性分析,所以台基剖面将固定为实验剖面.

图1 OMFIT 芯部等离子体剖面集成模拟流程图Fig.1.The integrated simulation workflow of core plasma with OMFIT.

考虑到电流扩散时间尺度远大于输运时间尺度,为了数值上的准确性,集成模拟中采用输运时间尺度作为迭代的时间步长.一个完整的电流演化过程需要非常多的时间迭代步数才能实现,因此在单一程序内同时对电流分布和等离子体动理学剖面如电子密度ne、电子温度Te、离子温度Ti等进行高效地演化是不现实的.这里将计算电流演化和等离子体动理学剖面演化的过程分开计算,由ONETWO进行各种源项计算和电流演化,而等离子体动理学剖面的演化则由包含了芯部湍性和新经典输运的TGYRO 程序进行计算.

将电流和输运分开计算的基本假设是,在特定条件下的平衡态物理计算过程中,电流和输运可以分开计算多次迭代集成达到最终平衡态,等离子体可以在一个远小于电阻扩散时间尺度的时间内,输运达到热平衡状态.这样我们可以按照在电流演化的时间尺度上取相应的步长,电流每演化一个步长,直接计算新的电流分布下输运能够达到的稳定状态,从而避免了按时间演化输运所需要的大量计算成本,同时保证了所关心的平衡状态演化的计算精度.该流程无需对小尺度输运过多关注,完成对电流分布和动理学剖面演化后,即可在此基础上引入新模块对更为细致的物理问题进行研究.

具体的集成模拟工作流程如下:

1) 从HL-2A 放电平顶段的某一时刻t0出发,将该时刻处理后的实验数据导入OMFIT 集成平台,包括电子密度ne、电子温度Te、离子温度Ti、杂质离子密度nimp或有效电荷Zeff、辐射功率Qrad、环向旋转速度Ω的剖面分布,以及等离子体电流Ip、环电压V、环向磁场强度Bt和平衡位形等.

2) 保持等离子体动理学剖面不变,在一个很短的输运时间步长τ内使用ONETWO 计算电流分布,同时由ONETWO 调用辅助加热计算程序来计算热源项.由于本文所分析的实验只有中性束注入(neutral beam injection,NBI)加热,因此只需调用NUBEAM 程序来计算驱动电流以及快离子分布.在HL-2A 的实验中,由于能量约束时间通常在10 ms 的量级,为了能够根据剖面及时更新加热、电流驱动,因此将时间步长设置为1 ms.

3) 把ONETWO 计算得到的压强p和电流分布FF′剖面传递给EFIT,由EFIT 进行平衡计算,得到新的电流分布下的等离子体平衡.

4) 将2)得到的热源项和3)得到的等离子体平衡交给剖面演化TGYRO 程序,由TGYRO 计算芯部区域的等离子体动理学剖面,如电子密度ne、电子温度Te、离子温度Ti.

5) 将实验上的台基区等离子体剖面与TGYRO计算得到的芯部等离子体剖面相结合,便得到归一化小半径从0到1 上的完整等离子体剖面.

6) 基于更新后的剖面和平衡,重新循环进行第2)—5)步的操作,在迭代N个循环后,等离子体剖面计算达到收敛,最终得到一个稳态的解.

2.2 物理模型

EFIT 是托卡马克平衡重建最常用的程序,其原理是计算磁面坐标下的Grad-Shafranov 方程:

其中p是压强;F是电流通量;ψ是约化后的极向磁通.在给定初始值p(ψ)和F(ψ) 后,即可由EFIT程序计算得到二维的平衡位形.

ONETWO和TGYRO 都是计算输运的程序,但原理上不尽相同.ONETWO 采用GLF23[26]和一些简化的新经典输运模型分别计算湍流输运系数和新经典输运系数,TGYRO 则是采用TGLF和NEO 程序来分别计算湍流输运和新经典输运.在得到输运系数后,ONETWO 求解输运方程来获得随时间演化的温度、密度、电流剖面,TGYRO[6,7]则是通过通量匹配的方法来得到演化终态的温度、密度、环向旋转剖面.如2.1 节所述,采用TGYRO来模拟温度、密度剖面,用ONETWO 来模拟电流剖面.下面详细介绍TGYRO 相关的物理模型与求解方法.

TGLF 既可以作为线性程序[10]计算不稳定性波模的本征函数、频率和增长率,也可以作为准线性输运模型[27]计算湍性输运通量,包括粒子、能量和环向动量通量等.在集成模拟工作流程中,TGLF以准线性模型得到湍性输运通量,传递给TGYRO进行剖面演化计算.目前TGLF 湍性通量计算仅适用于捕获离子模(trapped ion mode,TIM)、捕获电子模(trapped electron mode,TEM)、离子温度梯度模(ion temperature gradient,ITG)、电子温度梯度模(electron temperature gradient,ETG)的漂移波湍流,而无法计算由动理学气球模(kinetic ballooning mode,KBM)、环向阿尔芬本征模(toroidal Alfvén eigenmode,TAE)以及高能量粒子模(energetic particle mode,EPM)等引起的湍性输运.TGLF 的饱和准则通过与GYRO 或CGYRO非线性动理学模拟结果校核而确定的,发展出SAT0[27]和SAT1[28]两种不同的饱和准则.通常认为SAT0更适用于低q95的情形(如q95≤6),后者适用于高q95的情形[29].本文研究的实验q95～4,因此选择使用SAT0.

新经典物理程序NEO[30-32]是直接求解漂移动理学方程的第一性原理程序,基于线性Fokker-Planck 碰撞算子,通过直接求解一系列由DKE 推导出来的方程来得到等离子体的分布方程,进而得到相应的输运系数.它可以比较精确地求解新经典效应主导的任何物理量,比如新经典输运通量、自举电流、等离子体旋转等.在计算新经典输运通量方面,NEO 与流体程序NCLASS和解析公式Chang-Hinton 都有很好的校核,但NEO 比后两者的适用范围都广.

TGYRO 可以快速地将输运计算到稳态,忽略中间的演化过程.TGYRO 通过调用TGLF和NEO分别计算湍性输运通量和新经典输运通量,实现总通量Qσ(σi,e)与源项计算得来的目标通量相匹配[32].

给定rr*处的边界条件Tσ(r*)后,即可积分反演得到整个剖面:

3 实验结果验证

3.1 初始实验等离子体温度和密度剖面

HL-2A 托卡马克是一个具有上、下封闭式偏滤器的中型聚变研究装置,其大半径为R1.65 m,小半径为a0.4 m .通常情况下,装置运行在下单零偏滤器位形.典型的等离子体放电参数为:等离子体电流IP180—200 kA,线平均电子密度ne(1.5-4)×1019m-3,环向磁场Bt1.2—1.6 T .此外,在HL-2A 装置上具有1.5 MW 的中性束注入加热系统(NBI),5 MW 的电子回旋加热系统(ECRH)和1 MW 的低杂波电流驱动系统(LHCD).基于HL-2A 的第37012 炮高比压高约束模放电进行模拟验证,加热方式只有NBI 加热,无射频波加热,诊断数据较为齐全,其放电参数如图2 所示.根据放电参数,选取稳态运行阶段的某一时刻作为模拟的初始时刻t0,这里t0选取的是进入H 模后的1020 ms,该时刻具体的参数见表1.第37012 炮是常规偏滤器下单零位形放电,采用两束NBI 进行加热,等离子体参数相对较高.

表1 第37012 炮在1020 ms 时的参数Table 1.The parameters of the shot #37012 at the 1020 ms.

图2 第37012 炮放电参数(a) 等离子体电流 Ip ;(b) 等离子体储能 WE;(c) 归一化比压 βN和极向比压 βp ;(d) 线平均电子密度 ;(e) NBI 加热功率;(f)DαFig.2.The dischargement parameters of the shot #37012:(a) Plasma current Ip;(b) stored energy WE ;(c) normal-izedbeta βNand poloidal beta βp;(d)line-averaged elec-tron density;(e) NBI heating power PNBI;(f)Dα.

实验诊断获得的动理学剖面主要有电子密度ne、电子温度Te、离子温度Ti、离子旋转速度ω.电子密度由FIR和微波反射两种诊断方式进行测量,电子温度由ECE 诊断进行测量,离子温度和旋转速度则由CXRS 系统给出.CXRS 测量 C6+的光谱,通过假定氘离子与碳离子的温度与旋转速度相等来获取氘离子的信息.由于芯部电子密度较高,ECE诊断中出现了截止,因此芯部电子温度也将结合离子温度以及等离子体总储能来近似构建.对第37012 炮1020 ms 时刻的离子温度和旋转速度数据处理如图3 所示.

图3 第37012 炮在1020 ms 时的(a)离子温度和(b)旋转速度剖面Fig.3.The ion temperature profile (a) and rotation profile(b) of the shot #37012 at the 1020 ms.

电子密度有两种诊断方式,在边界通过微波反射进行测量,精度较高,可以在排除异常数据点后直接使用,而靠近磁轴的电子密度是通过激光干涉进行测量,数据点较少,可以提供弦积分密度信息.对电子密度的处理则是以微波反射测量值为准,同时参考激光干涉弦积分后的芯部剖面,在实验误差允许的范围内,二次重建一个同时满足反磁测量得到的线平均密度和微波反射诊断数据的电子密度剖面,如图4 所示.

图4 第37012 炮在1020 ms 时的电子密度剖面处理Fig.4.The treatment of electron density profile of the shot#37012 at the 1020 ms.

在诊断数据不够齐全且不自洽时,采纳精度高、数据全的诊断数据,通过构造的方式得到与宏观参量(线平均密度、储能、βN)相匹配的且诊断未给出的电子温度剖面.本文参考离子温度剖面形状,经过计算,最终得到与其他测量结果自洽的电子温度剖面,如图5 所示.在3.2 节中,使用图3 拟合后的离子温度和旋转实验剖面、图4 二次优化后的电子密度实验剖面,以及图5 构造的电子温度剖面,一起作为OMFIT 集成模拟验证的初始输入.

图5 第37012 炮在1020 ms 时得到的电子温度剖面Fig.5.The profile of electron temperature of the shot#37012 at the 1020 ms.

3.2 集成模拟与实验结果对照与分析

按照2.1 节中介绍的计算流程,在HL-2A 装置上对第37012 炮进行了模拟验证.在ONETWO计算中,采用耦合的NUBEAM 程序对中性束粒子进行追踪计算.模拟时间段是1020—1030 ms,在这个时间段内等离子体处于平顶段稳定运行阶段,各个等离子体参数几乎没有变化.

由于台基区物理的复杂性以及诊断上的误差,目前对台基区域的计算和对照不是很可靠.在对第37012 炮的模拟中,芯部 (ρ＜0.8) 的湍流输运和新经典输运采用TGYRO 进行计算,而在边界处 (0.8＜ρ＜1) 固定剖面,与实验保持一致.在第37012 炮的实验诊断里,缺少对杂质密度或有效电荷的诊断,这里按照其他实验炮的诊断,将有效电荷设置为全空间恒定的常数2.6.

经过多轮迭代后,各项宏观参数及动理学剖面演化趋于收敛.从图6—图8 中可以看出,温度、密度以及压强略有下降,磁轴处温度最大误差也控制在10%以内,但整体上与实验上吻合得很好,可以证明以上的集成模拟流程可信.对各种成分的电流剖面计算结果如图7 所示,可以得到自举电流的份额为52%,NBI 驱动电流份额在30%左右,欧姆电流为18%.

图6 集成模拟计算中各物理量的多次迭代收敛性 (a1),(b1),(c1) 迭代前后对比;(a2),(b2),(c2) TGYRO 计算点的收敛过程Fig.6.The astringency of each physical quantity in the integrated simulation:(a1),(b1),(c1) The comparison between before and after the iteration;(a2),(b2),(c2) the convergence process of the TGYRO calculating points.

图7 第37012 炮在1020 ms 时刻的各成分电流剖面Fig.7.The current profiles of each composition of the shot#37012 at the 1020 ms.

图8 第37012 炮在1020 ms 时刻剖面的模拟结果与实验结果对照 (a) 压强剖面;(b) 电子密度剖面;(c) 离子温度剖面;(d) 安全因子 q 剖面Fig.8.The experiment and simulation profiles comparation of the shot #37012 at the 1020 ms:(a) Pressure;(b) electron density;(c) ion temperature;(d) safety factor q .

在第37012 炮中,由两束同向的NBI 对主等离子体进行加热,总功率为1.4 MW.经过NUBEAM程序计算,可以从图9 中看出,本炮放电采取的NBI 离轴注入,能量沉积位置主要在归一化小半径为0.2 处附近,且对离子成分的加热效果更加显著.

图9 NBI 能量密度沉积分布Fig.9.The distribution of NBI deposed energy density.

4 芯部输运分析

在计算达到稳定状态后,对TGYRO 最终结果进行分析,从图10 可以看出,TGYRO 计算的总通量与目标通量可以相匹配,满足(2)式—(4)式.电子能量输运通道是湍性输运(即TGLF 计算出的湍性通量)占据主导,而离子能量输运通道上,湍性和新经典输运(即TGLF 计算出的新经典通量)所占份额几乎相当,在ρ0.2-0.3 附近的ITB(internal transport barrier)区域氘离子能量输运由新经典输运占主导,湍性输运得到抑制.

图10 第37012 炮放电在1020 ms 时刻模拟后得到 (a) 离子能量通量;(b) 电子能量通量Fig.10.The ion energy flux (a) and electron energy flux(b) of the shot #37012 at the 1020 ms.

使用TGLF-scan 模块对上述结果进行静电漂移波不稳定性分析,图10 为0.2—0.9 区间内两支最不稳定的本征模式扫描结果,横坐标为归一化小半径,纵坐标为归一化波数kθρs,颜色深浅代表不稳定性增长率的高低,而红色代表的是电子抗磁漂移方向 (ω ＞0),蓝色代表的是离子抗磁漂移方向(ω＜0).通过图11 的波数分析进行模式识别,在0.2—0.3区域内,高k的不稳定性受到抑制,主要是ITG 模主导的湍流,在0.7 以外,有TEM和ITG 得到增长.

图11 0.2—0.8 区域内两支最不稳定的本征模式的频谱Fig.11.The spectrum of two most unstable eigenmode in the 0.2—0.8 region.

在经过TGLF 对不同位置(ρ0.3,0.5,0.8)的微观不稳定进行计算后,给出相应的增长率与波数的关系,如图12所示.图中横纵坐标分别为归一化的波数kθρs和γ(cs/a),其中cscs/Ωs,ΩseB/mic,a为装置小半径.按照线性分析的结果来看,高k模式在ρ0.3 附近受到抑制,湍流由红色小圆点代表的ITG 主导.在偏离ITB 区域(ρ0.5)内,ETG 逐渐增长,成为湍流的主要因素.而在更靠近边界处(ρ0.8),TEM和ITG 都得以增长,但主导输运的仍是ETG 不稳定性.

图12 ρ=0.3,0.5,0.8 处线性不稳定性的增长率与波数的关系(蓝色为电子抗磁漂移方向,红色为离子抗磁漂移方向)Fig.12.The relationship between the growth-rate and wavenumber of the linear instabilities in the ρ=0.3,0.5,0.8(the blue points represent the electron diamagnetic drift direction and the red points represent the ion diamagnetic drift direction).

5 结论

在OMFIT 平台上,通过集成EFIT,ONETWO/NUBEAM,TGYRO/TGLF/NEO 得到了可以实现对HL-2A 装置第37012 炮高比压放电实验结果的集成模拟验证流程.在该流程中,验证了HL-2A 装置的程序适用性和可靠性,最终芯部集成模拟结果整体上可以很好地与实验诊断数据相吻合,磁轴处温度最大误差仍然控制在10%以内.在第37012 炮的模拟结果基础上,进一步采用TGLF 模型进行线性不稳定性分析,结果显示NBI 沉积位置附近的高k微观ETG 不稳定性受到抑制,输运由低k的ITG 占据主导,输运水平降低至接近新经典的水平,从而实现了NBI 离轴加热导致的H模约束改善.

目前的结果尚未考虑台基和边界区物理,无法对边界杂质输运影响进行更多细致分析.同时,这里需要更为全面详细的诊断数据作为支撑,未来将在本文工作基础上进一步完善对HL-2A 实验的集成模拟验证流程.

感谢GA 方面提供的OMFIT 平台支持和GACODE,感谢Orso Meneghini、毛瑞、李航、吴木泉对本文工作的指导、建议和帮助.