EAST托卡马克杂质输运蒙特卡洛模拟

查学军, 钟德俊, 王福琼, 陈一平, 卢洪伟, 胡立群（.东华大学理学院应用物理系,上海 060；.中国科学院等离子体物理研究所,合肥 3003）

EAST托卡马克杂质输运蒙特卡洛模拟

查学军1, 钟德俊1, 王福琼2, 陈一平2, 卢洪伟1, 胡立群2
（1.东华大学理学院应用物理系,上海 201620；2.中国科学院等离子体物理研究所,合肥 230031）

介绍DIVIMP程序的基本思想和采用的物理模型,研究EAST托卡马克欧姆放电条件下碳杂质的产生与输运.模拟结果与实验测量具有较好的一致性,说明将该程序应用于EAST装置中杂质输运模拟的合理性.并利用DIVIMP程序预测EAST在8MW加热功率条件下钨杂质的输运特性.

托卡马克；杂质输运；模拟；蒙特卡洛方法；DIVIMP

0 引言

托卡马克[1]是一种利用磁约束来实现受控热核聚变的环形装置.托卡马克边界等离子体与第一壁材料发生相互作用,所产生的杂质能通过边界等离子体的输运进入主约束区等离子体,并对其性能造成极大影响.一方面,输运进入聚变反应堆芯部的杂质通过稀释燃料浓度来降低聚变功率,同时通过杂质辐射降低等离子体温度,缩小等离子体电流通道,从而引发不稳定性,影响等离子体约束；另一方面,存在于托卡马克装置偏滤器区域的杂质,其辐射过程有利于降低到达偏滤器靶板的能流,进而又会减小偏滤器靶板材料的腐蚀及杂质的产生.因此,杂质输运一直是磁约束等离子体最重要的研究课题之一[2].托卡马克装置中杂质产生及输运特性的研究以及杂质控制是磁约束受控热核聚变能源商用化进程中必不可少的环节.大科学工程装置EAST是由中国科学院等离子体物理研究所建立的世界上第一个具有非圆截面全超导托卡马克实验装置.EAST边界等离子体杂质输运研究的目的是为了合理设计并优化第一壁和偏滤器靶板部件,更好地控制等离子体中的杂质.

实验诊断和程序模拟是研究托卡马克装置中杂质特性的两种主要方法[3].杂质输运与边界等离子体输运强烈的相互耦合：等离子体施加给杂质粒子的力控制了杂质的输运过程,而杂质辐射又会影响等离子体参数的分布.杂质与等离子体之间的这种相互作用最终决定了杂质密度剖面.由于,大多数杂质相关实验诊断,如光谱仪和热辐射仪等,仅测量了沿着观测线的积分信号,测量过程掩盖了许多杂质分布信息.刮削层中的等离子体温度及流速场分布是决定杂质输运的重要因素,但是目前大多数托卡马克装置中几乎没有有关这些物理量在整个刮削层中完整分布的实验测量.杂质分布及等离子体参数剖面信息的缺失,必须通过模拟手段予以弥补.

当然,近年来杂质相关实验诊断也得到了较好的发展.测量等离子体流速的马赫探针[4－5]以及测量偏滤器区域等离子体温度、密度的汤姆逊散射系统[6]相继在各托卡马克装置中成功开展.在许多托卡马克中,已经通过多通道热辐射测量仪（Bolometer）信号[7]成功建立了总辐射的二维（如图1所示的二维极平面）剖面,并通过光谱诊断[8]建立了杂质辐射的二维剖面.此外,在DIII-D等托卡马克装置中,已经能够通过大量的偏滤器测量推断杂质密度在极向平面的二维分布[9].上述各二维剖面信息均将为我们提供理解杂质输运的有效信息.但是,一些诊断技术的使用受到物理或工程限制,如Bolometer所采用的断层分析算法具有很大的不确定性,且仅能针对一些特定的电荷态.此外,目前边界等离子体温度测量主要限于测量电子温度,对离子温度的直接测量结果很少,而离子温度梯度力通常又主导了杂质离子向芯部的输运过程.因此,为了解释实验信息,详细理解托卡马克边界物理活动,程序模拟包括对杂质输运特性的模拟仍然必不可少.

综上所述,程序模拟不仅可以克服杂质相关实验诊断信息和诊断区域有限的缺点,还是解释实验诊断杂质信息的有效工具.另外,实验与模拟结果的对比分析,可以帮助证实模拟中所采用模型的合理性.杂质输运程序主要有流体程序和蒙特卡洛程序两大类[10].蒙特卡洛杂质输运程序可以方便地模拟发生在托卡马克边界的复杂分子、原子过程,不需要对托卡马克位形作出除环向对称以外的任何假定和近似.同时,由于杂质粒子尤其是低价态杂质离子的特征电离时间往往短于粒子之间的平均碰撞自由时间,流体近似条件[11]可能得不到很好的满足.此外,对于一些核电荷数较高的杂质粒子（如74钨）的输运过程,如采用流体程序模拟,需要考虑的离子种类及需要求解的方程总数较多,计算时间往往较长,而抽样跟踪单粒子运动轨道及电荷态的蒙特卡洛杂质输运程序可克服该缺点.因此,蒙特卡洛程序在模拟杂质输运特性方面,有其优越性.

实验先进超导托卡马克（EAST）[12]于2006首次投入运行,并于2010年成功实现高约束运行模式.目前,该装置尚处于运行初期阶段,一些诊断手段,特别是边界过程相关诊断,还不完善,如边界汤姆逊系统、测量等离子体密度剖面的束发射光谱诊断等尚未成功在EAST装置中予以开展.同时,尽管已成功开展一些光谱诊断,如Da信号及CIII信号的诊断等,但杂质相关诊断信息仍非常有限,并且没有边界等离子体参数二维分布信息.目前,CIII信号观测通道仅涵盖了外偏滤器区域,且没有除CIII以外的其它C杂质线辐射信号（如CII辐射等）的测量.此外,现有的实验测量也难免会有一些误差.因此,针对EAST装置开展有关杂质输运特性的模拟研究,不仅可在实验条件尚不完善的条件下,理解该装置杂质输运基本特性,还可通过模拟结果与现有实验数据之间的比较,检验模拟所采用的模型,为模拟工作的进一步开展奠定基础.

本文运用DIVIMP程序[13]模拟研究EAST托卡马克在欧姆加热情况下碳杂质的输运特性,并得到其相应的碳杂质二维空间分布.基于碳杂质密度分布信息和背景等离子的温度密度信息,计算碳杂质辐射分布.实验测量的CIII信号与相关模拟结果之间有很好的一致性,说明将DIVIMP用于EAST装置上对杂质输运特性的解释和预测的合理性.由于钨材料具有熔点高、耐腐蚀性强等优点,国际热核聚变实验堆（ITER）、聚变示范堆（DEMO）等运行参数更高的下一代聚变装置的第一壁材料均考虑采用钨材料.EAST托卡马克偏滤器材料也将逐步升级改造为钨材料.于是,采用DIVIMP程序预测了EAST高功率加热条件下钨杂质的输运特性.模拟结果可为EAST第一壁和偏滤器的优化升级和实验的优化设计提供参考[14].

1 模型与方法

DIVIMP是通过抽样跟踪偏滤器位形下托卡马克边界等离子体中杂质离子状态和轨道的蒙特卡洛程序.与传统的流体程序相比,DIVIMP考虑的是单个粒子的输运,模拟它们在物质中运动的历史,记录其在运动中相关物理模拟量的贡献,在对单个粒子运动历史进行大量的重复模拟之后,再对物理模拟量进行统计平均,得到所需要的物理结果.它与传统的解法相比,程序结构简单,易于实现,误差较小.以下是对DIVIMP中杂质产生和输运计算模型的详细描述.

图1 单零偏滤器位形Fig.1 Skematic of a single null divertor configuration

1.1 杂质的产生

图1是单零偏滤器位形下,托卡马克磁约束等离子体的二维剖面.附加环向（图1中垂直纸面方向）导体中有与等离子体电流（IP）同向的电流（ID）,由ID和IP所产生的极向磁场零点称为X点,经过X点的磁面为分界面.在分界面外,磁力线开放的区域为刮削层（Scrape-off Layer,SOL）.刮削层中的开放磁力线与偏滤器靶板直接相交,如图2所示[15].当背景离子沿着磁力线运动到达偏滤器靶板表面时,靶板材料会因为溅射等过程向等离子体中发射杂质.同样,返回偏滤器表面的杂质离子也可能引发自溅射.

溅射过程可用溅射产额Y来描述,其定义为平均每入射一个粒子从靶板表面溅射出来的原子数,即

图2 偏滤器区域杂质产生和输运过程Fig.2 Production and transport processes of impurities in the divertor region

模拟中,溅射产额Y对入射粒子能量的依赖采用2007年由Eckstein等人所推荐的关系式[16].假定背景离子的入射能量Eimp为

其中,ZB为等离子体离子的电荷数,Te0和Ti0分别为偏滤器靶板处的电子和离子温度.由背景等离子体中离子的轰击所引起杂质粒子的总溅射通量为

1.2 杂质粒子的输运

1.2.1 中性杂质粒子的输运

等离子体与偏滤器靶板相互作用产生的杂质粒子均以中性粒子的形式进入等离子体中.与杂质离子不同,中性杂质粒子的运动不受电磁场约束,在DIVIMP模拟中分别由不同的程序模块对二者的运动状态进行单独跟踪.在t＝0时刻,中性杂质粒子以上述计算所得概率从偏滤器部件发射,并遵从汤普森能量分布和余弦角分布.在中性粒子被电离、或打击到固体表面之前,他们保持直线运动.模拟中假定,当中性粒子打击到偏滤器靶板表面,其被反射,而当中性粒子打击到除靶板以外的第一壁表面则被吸收.在时间步长Δtn时间内,中性粒子发生电离的概率

1.2.2 杂质离子的输运

杂质离子产生后,可能经历电离、复合、电荷交换等原子、分子过程.定义：

杂质离子沿着磁力线方向的运动采用经典理论描述.模拟中,假定离子沿磁力线方向受力情况为

其中,等式右边的五项分别为杂质压力梯度力、与背景等离子体之间的摩擦力、电场力、电子温度梯度力、离子温度梯度力等,s为沿磁力线的距离,p为杂质离子的压强,VB为背景等离子体的流速,v为杂质离子运动速度,Z为杂质离子的电荷数,αe和βi分别为电子和离子温度梯度系数.从方程（6）可见,杂质离子的输运与背景等离子体流速场和电场、电子温度梯度、离子温度梯度等物理量有关.模拟中,这些背景等离子体信息由SOLPS（B2－EIRENE）[17]计算得到.假定杂质离子在垂直磁力线方向上的运动为反常扩散.反常扩散系数一般不能由第一性原理导出.根据相关参考文献和国际上其它托卡马克中杂质输运特性模拟所采用的相关假定[18],模拟中规定杂质离子垂直磁力线反常扩散系数D＝1.0 m2·s－1,并且不考虑杂质离子的径向内箍缩效应.在经历每个时间步长Δt后,DIVIMP分别记录每个杂质离子的电荷态、位置、速度等信息.达到稳态以后,计算得到杂质密度和杂质辐射等信息.

1.3 杂质密度

图3 DIVIMP的计算网格Fig.3 Grid for DIVIMPmodeling

DIVIMP计算将边界等离子体区域划分为计算网格,如图3所示,对杂质、粒子运动状态的跟踪和杂质相关信息的计算基于网格进行.网格基于实际托卡马克放电时的平衡位形.可以看出,DIVIMP计算区域包括整个刮削层区域和分界面以内的芯部外围区域.

每一个网格由沿着磁力线的x坐标和垂直磁力线的y坐标来描述.t＝0时刻发射的N个杂质离子,经过ts后位于网格（x,y）内并具有ZI电荷态杂质离子归一化密度为

当发射的第i个离子进入（x,y）网格并处于ZI电荷态时,Ci（x,y, ZI,t）＝1（对于自溅射Ci（x,y,ZI,t）＜1）；否则,Ci（x,y,ZI,t）＝0.达到稳态时,（x,y）网格内的电荷态ZI杂质离子的归一化密度为

因此,（x,y）网格内ZI电荷态杂质离子的绝对密度为（x,y,Z1）R,考虑到部分中性粒子直接沉积到壁表面而对离子数没有贡献,R为电离产生离子的中性粒子份额.

1.4 杂质辐射功率

DIVIMP中一项重要的输出信息是,单位时间内由某一电荷态杂质离子在等离子体中给定位置附近单位体积内所发射的具有特定波长（频率）光子数ΦPLRP.根据ΦPLRP信息和ADAS数据库中有关一个特定电荷态杂质离子在其被电离到其它电荷态之前所能释放的光子数PE（光子发射效率）,可以推断该电荷态杂质离子的产生速率

1）电子碰撞引起的杂质离子从基态的跃迁.这是最基本的过程；

2）电子－离子复合,该过程也可能导致激发态的产生,从而会释放光子；

3）从上一电荷态的电离,形成激发态并伴随光子发射；

4）电荷交换复合,此时需要输入二维中性氢粒子分布相关信息（如来源于ERIENE的计算等）.

在实际计算中,可根据求解问题的情况仅仅考虑上述过程中的某些过程.如只考虑过程1）,则

其中η为发射一个光子所发生的平均电离事件次数,Iiz（x,y,z1）为单位时间发生电离的杂质离子数.

其中τiz为杂质离子的电离特征时间.

本文首先采用上述介绍的计算模型,模拟了EAST托卡马克第12303炮放电,欧姆加热条件下C杂质的输运特性.与模拟相关的主要参数如表1所示.表1中,安全因子定义为磁面上所躺磁力线在极向转一圈后,沿着环向所绕行的圈数,q95为95％磁通面上的安全因子,环向比压βT＝〈P〉/[/（2μ0],极向比压βP＝〈P〉/[/（2μ0）],BT和BP分别为环向磁场和极向磁场,〈P〉为等离子体体积平均热压力,线平均密度为等离子体密度在沿R方向过磁轴的直线上的平均值.其它等离子体参数的定义如图4所示.然后,本文采用DIVIMP程序预测了EAST装置在8MW高功率加热条件下钨杂质的输运特性.

表1 与模拟研究相关的主要参数Table 1 M ain parameters used in simulation

图4 托卡马克磁面位形极向剖面Fig.4 Poloidal cross-section ofmagnetic surfaces in a tokamak

2 结果与讨论

图5是模拟得到的碳杂质密度在边界区域的二维分布.由图5可见,杂质离子的二维分布沿着磁力线方向和垂直磁力线方向均存在着不均匀性.杂质密度垂直磁力线方向（径向）不均匀性可由Engelhard模型定性解释.而杂质密度沿着磁力线的不对称性主要由存在于边界区域中的许多复杂杂质过程,如电离、复合等所导致.另外,由图5可见,欧姆加热条件下由于加热功率非常低,即使对于较容易被腐蚀的碳杂质在芯部等离子体中的最高密度约为1015量级.

如1.4节所述.基于DIVIMP计算得到的杂质密度信息和SOLPS所提供的等离子体参数信息,DIVIMP可输出由C＋2离子辐射产生的CIII信号（波长465.0nm）的二维分布.图6是CIII辐射信号沿着经过偏滤器区域的18个观测通道的积分信号,横轴1－18分别对应于不同的观测通道.值得注意的是,图6中的实验测量值为未绝对标定的相对值,而计算值为实际值.由图6可见,模拟结果与实验结果有较好的符合,说明了将DIVIMP程序用于EAST装置杂质输运特性模拟的合理性.

图5 DIVIMP计算所得碳杂质密度在边界的二维分布Fig.5 Two-dimensional distribution of C impurities in the edge of EAST from DIVIMP calculation

图6 CIII信号的计算值和测量值对比Fig.6 Comparison against calculated CIIIsignals and measured ones

利用DIVIMP程序预测了未来EAST装置在8MW加热功率条件下钨杂质的输运特性,结果如图7所示.在8MW的高功率加热条件下,芯部等离子体中钨杂质的最高密度约为1015量级.尽管随着加热功率的升高,偏离靶板处的等离子体温度及离子轰击靶板时的入射能量也将升高,但具有最高束缚能的钨原子（8.68eV）,其溅射产额没有明显的增加.此外,由于钨杂质粒子的质量较大,其从偏滤器靶板溅射的出射速度较低,而其特征电离时间又非常短,也就是说,溅射产生的钨杂质原子,大多数将在偏滤器靶板附近电离,存在于偏滤器区域的高强粒子流的冲刷,以及靶板表面鞘层电压对钨离子的加速作用,使得这些钨离子重新返回偏滤器靶板[19].上述各原因使得,即使当加热功率大幅提高（～50倍）时,如将偏滤器靶板材料由碳材料改变为钨材料,等离子体中杂质水平不会有太明显的提高.

图7 DIVIMP计算所得钨杂质在边界层的密度分布Fig.7 Two-dimensional distribution ofW impurities in the edge of EAST from DIVIMP calculation

3 结论

用DIVIMP程序模拟了EAST欧姆加热条件下的边界层碳杂质的输运,得到了相对应的杂质分布,并将计算得到的杂质辐射信息与实验测量信息进行比较.结果表明：（1）欧姆加热条件下,EAST装置中碳杂质的最高水平约为1015量级.（2）模拟结果与实验结果之间符合较好,证明将蒙特卡洛程序DIVIMP用于研究EAST装置杂质输运特性模拟的合理性.采用DIVIMP程序模拟预测了未来EAST在8 MW高功率加热条件下钨杂质的输运特性.结果发现,高加热功率条件下,如将偏滤器靶板由碳材料改为钨材料,杂质含量不会有明显的提高.

[1] 石秉仁.磁约束聚变原理与实践[M].北京：原子能出版社,1999.

[2] 徐伟,万宝年,谢纪康.HT-6M托卡马克装置杂质输运[J].物理学报,2003,52（8）：1970－1909.

[3] Shimizu K,Takizuka T,Sakasai A.A review on impurity transport in divertors[J].Journal of Nuclear Materials,1997,241－243：167－181.

[4] Stansfield B L,Boucher C,Gunn J P,et al.Measured plasma parameters in TdeV's closed poloidal divertors[J].JNucl Mater,1995,220－222：1121－1125.

[5] LaBombard B,Jablonski D,Lipschultz B,McCracken G,Goetz J.Scaling of plasma parameters in the SOL and divertor for Alcator C-Mod[J].JNucl Mater,1995,220－222：976－981.

[6] Allen SL,Brooks N H,Campbell R B.Developmentof a radiative divertor for DIII-D[J].JNuclMater,1995,220－222：336 －341.

[7] Reichle R,Clement S,Gottardi N.Controlled fusion and plasma physics[C]∥Proc 22nd EPS Conf,The European Physical Society,Bournemouth,UK,July 3-July 7,1995,19C（Ill）：85.

[8] Allen S L,Bozek AS,Brooks NH.Recent DIII-D divertor research[J].Plasma Phys Control Fusion,1995,A191.

[9] Fenstermacher M E,Porter GD,Rensink M E,etal.UEDGE and DEGASmodeling of the DIII-D scrape-off layer plasma[J].JNucl Mater,1995,220－222：330－335.

[10] Stangeby PC.Dynamics of impurities in the scrape-off layer[J].Contrib Plasma Phys,1988 28：507.

[11] Wesson J.Tokamaks[M].International Series of Monographs on Physics,Oxford University Press,2011.

[12] Wang FQ,Chen Y P,Hu LQ,et al.Onion skinmodel（OSM）analysis of EAST SOL plasmas[J].Nucl Fusion,2014,54：093002..

[13] Stangeby PC,Elder JD.Calculation of observable quantities using a divertor impurity interpretive code,DIVIMP[J].JNucl Mater,1992,258：196－198.

[14] 郭后扬.JET中偏滤器几何结构和化学溅射对杂质行为和等离子体性能的影响[J].国外核聚变,2001,4：48－63.

[15] Stangeby PC.The plasma boundary ofmagnetic fusion devices[M].Bristol：Institute of Physics Publishing,2000.

[16] Behrisch R,Eckstein W.Sputtering by particle bombardment：Experiments and computer calculations from threshold to MeV energies[M].1st ed.Topics in Applied Physics.Springer,Softcover Reprint of Hardcover,2007.

[17] Schneider R,Bonnin X,Borrass K.Plasma edge physicswith B2-eirene[J].Contrib Plasma Phys,2006,46（1/2）：3－191.

[18] Schmid K,Krieger K,Kukushkin K,Loarte A.DIVIMPmodeling of tungsten impurity transport in ITER[J].JNucl Mater, 2007,363－365：674－679.

[19] Järvinen A,Giroud C,Groth M,Krieger K,Moulton D,Wiesen S,Brezinsek S,JET-EFDA Contributors.DIVIMP simulation ofW transport in the SOL of JET H-mode plasmas[J].Phys Scr,2011,T145：014013.

M onte-Carlo M odelling of Im purity Transport in EAST Tokamak

ZHA Xuejun1, ZHONG Dejun1, WANG Fuqiong2, CHEN Yiping2, LU Hongwei1, HU Liqun2
（1.Department ofApplied Physics,Donghua University,Shanghai 201620,China；
2.Institute of Plasma Physics,Chinese Academy ofSciences,Hefei 230031,China）

Basic idea and physicalmodels in DIVIMP are introduced and DIVIMPmodelling of C impurity production and transport properties in edge plasma of EAST machine with ohmic heating are provided.Simulation results show favorable consistency with experimentalmeasurements,demonstrating that themodels coupled in DIVIMP are suitable for interpretive and predictivemodelling of impurity transport in EASTmachine.Transport properties ofW impurity in EASTwith heating power as high as8MW are also predicted by DIVIMPmodelling.

tokamak；impurity transport；simulation；Monte Carlomethod,DIVIMP

1001-246X（2015）06-0715-07

TL61＋2

A

2014－11－09；

2015－04－09

国家自然科学基金（11175045）和科技部ITER计划专项课题（2013GB109001）资助项目

查学军（1971－）,男,博士,副研究员,主要从事等离子体物理理论和数值模拟研究,E-mail：xjzha＠dhu.edu.cn