CFETR参数下α粒子慢化过程的数值模拟*

吴相凤 王丰 林展宏 陈罗玉 于召客 吴凯邦 王正汹

(大连理工大学物理学院,大连 116024)

氘氚聚变产生的高能量α 粒子是维持未来托卡马克反应堆等离子体高温的主要加热源,良好的α 粒子约束对于维持稳态燃烧等离子体至关重要.在持续发生聚变反应的系统中,α 粒子远离热平衡,呈现非麦克斯韦分布.如果忽略轨道效应,基于局域库仑碰撞的假设可以得到α 粒子的经典慢化分布,然而由于α 粒子存在较大的漂移轨道宽度,空间输运不容忽视,为得到更为准确的α 粒子分布函数,需要开展相关的数值计算.本文使用模拟程序PTC (particle tracer code)在中国聚变工程试验堆(CFETR)不同的放电模式下,采用粒子轨道跟踪和蒙特卡罗碰撞方法,对α 粒子慢化过程进行了数值模拟,获得了更为真实的α 粒子分布函数,并将其与经典慢化分布进行了对比.结果显示分布函数在中等能量附近和经典慢化分布存在较大差异.进一步的分析表明,这是由于中等能量下α 粒子的较强的径向输运引起的.本文的研究结果对准确评估α 粒子加热背景等离子体的能力具有重要参考价值.

1 引言

氘氚聚变反应是目前最有希望实现商业化应用的聚变反应之一[1,2].随着国际热核聚变实验堆(ITER)工程的持续推进,人们开始关注在ITER 之后的聚变装置的设计和建造问题,以便尽早地实现聚变能源商业化[3,4],中国聚变工程试验堆(CFETR)将作为ITER 和聚变电站商用示范堆(DEMO)之间的关键衔接点,以推进中国聚变能源发展,推动聚变能源商业化进程,以及实现可持续地清洁能源方面发挥重要作用[5,6].

在CFETR 燃烧等离子体中,因聚变生成的携带有 200 MW 聚变功率的α 粒子是加热等离子体的主要方式之一,然而,这些高能量的α 粒子只有被有效慢化才能起到维持等离子体自持燃烧的作用[7−10].在氘氚等离子体中关于3.5 MeV 的α 粒子在等离子体中的慢化分布或能量沉积的研究已经进行了半个多世纪[11−15].α 粒子慢化过程中,在非均匀磁场中漂移并受到电子和离子的摩擦以及与背景离子的角度散射,基于局域库仑碰撞假设求解福克-普朗克方程的稳态解可以得到经典慢化分布[16,17],它被广泛应用于包括α 粒子、中性束注入粒子等在内的快离子相关物理研究,如加热、电流驱动、快离子与Alfvén 本征模的相互作用及快离子输运等[18−21].然而这种慢化分布具有一定局限性,在真实情况下,轨道效应和空间输运是不可忽略的,同时,粒子漂移轨道宽度也会影响α 粒子在等离子体中的输运[20−22].这种简化模型无法准确描述α 粒子的行为,有必要采用更丰富的物理模型描述非局域的α 粒子慢化过程[23,24].

为了探索α 粒子的慢化过程,本文回顾了经典慢化分布理论模型并采用数值模拟方法进行对比分析.在数值模拟方面使用粒子轨道追踪耦合蒙特卡罗模拟程序PTC(particle tracer code)[25]在CFETR 稳态运行模式和混杂运行模式下对α 粒子的慢化过程进行了计算,模型考虑了粒子的有限轨道效应以及新经典输运,采用蒙特卡罗方法计算α 粒子与背景等离子体之间的库仑碰撞,得到了等离子体平衡条件下α 粒子的慢化分布函数.同时,统计了α 粒子热化和损失的比例,计算了α粒子对背景等离子体的加热功率,并得到了α 粒子的总体粒子约束时间和能量约束时间.本研究通过在不同运行模式下对α 粒子整个慢化过程中的α 粒子的产生、约束、对背景等离子体的加热以及整个分布演化的数值模拟分析,还有理论和程序模拟的对比分析,为CFETR 的物理设计提供了有价值的参考依据.

本文第2 节描述了经典慢化分布理论模型、CFETR 稳态和混杂运行模式下的等离子体参数、PTC 程序物理模型.第3 节是两种运行模式下α 粒子慢化过程的模拟结果,包括α 粒子的粒子数变化、对背景等离子体的加热功率、α 粒子的能量变化、α 粒子的慢化分布函数等.第4 节是讨论分析,包括两种运行模式下慢化过程的分析、约束时间的分析、理论与程序模拟的慢化分布函数的分析.第5 节是总结.

2 经典慢化分布与PTC 模型

2.1 经典慢化分布

经典慢化分布是指高能粒子与热背景等离子体碰撞导致高能粒子的经典减速形成的速度或能量分布,它可以作为描述α 粒子慢化过程的一种理论模型.经典慢化分布的推导基于忽略速度扩散项和轨道效应的假设,根据这个假设,可以简化福克-普朗克方程,得到如下式子[16]:

其中fα为α 粒子速度分布函数,Sα为α 粒子源项.此外,假设vi≪vα ≪ve,其中vi和ve为背景离子和背景电子的热速度,vα为典型的α 粒子速度,并假设背景电子温度与背景离子温度相等,即Ti=Te,有高能α 粒子的慢化碰撞算子:

其中聚变反应速率S0=nDnT〈σv〉,nD为氘的粒子数密度,nT为氚的粒子数密度,σ为反应截面,v为相对速度大小.对于能量为3.5 MeV 的α 粒子,初始速度vb=1.3×107m/s.将(2)式和(3)式代入(1)式并进行计算,得到稳态时三维的α 粒子速度慢化分布函数:

其中H(vb-v) 为阶跃函数,将(4)式乘以 4πv2,得到一维的α 粒子速度慢化分布:

根据能量E与速度v之间的变换关系以及阶跃函数的定义,得到一维的α粒子能量慢化分布:

其中,Eb为α 粒子初始能量,Ec为临界能量,根据(6)式,可得到CFETR 稳态运行模式下的局域稳态能量慢化分布如图1 所示.

图1 电子温度分别为27.78,14.4 和6.7 keV,对应电子密度分别为1.14×1020,9.34×1019 和7.47×1019 m−3 参数下得到的经典能量慢化分布f1, f2, f3Fig.1.Classical energy slowing-down distributions f1,f2 and f3 obtained for the electron temperatures of 27.78,14.4 and 6.7 keV,and their corresponding electron densities of 1.14×1020,9.34×1019 and 7.47×1019 m-3.

接下来,利用各向同性的α 粒子慢化分布计算出α 粒子的密度n和总能量nE,对分布函数(4)在三维空间积分,得到

根据粒子平衡方程和能量平衡方程:

接下来的模拟中,将使用PTC 程序对α 粒子的慢化过程进行模拟,并与理论模型进行比较分析.

2.2 CFETR 参数与PTC 物理模型

基于CFETR 稳态运行模式和混杂运行模式(v201903)参数下进行了α 粒子慢化过程的模拟,具体等离子体参数如图2 所示.由图2 可看到,在稳态运行参数下,芯部区域和边界区域等离子体参数相对比较均匀,混杂运行模式下的温度略高于稳态模式,密度量级与稳态参数一致.

图2 CFETR 中的背景等离子体参数 (a) 稳态运行模式下的密度、温度和安全因子剖面;(b) 混杂运行模式下的密度、温度和安全因子剖面Fig.2.Background plasma profiles in CFETR: (a) Density,temperature,and safety factor profiles in steady-state scenario;(b) density,temperature and safety factor profiles in hybrid scenario.

CFETR 稳态运行模式是一种基于等离子体稳定性和安全性的非感应电流驱动的运行模式,混杂运行模式是感应和非感应电流混杂驱动的脉冲式的运行模式.这两种运行模式不是独立的,它们可以相互转换和组合,以实现更优秀的聚变反应效果和能量输出.

PTC 程序是一个拥有自主知识产权的粒子追踪程序,可以模拟聚变反应中产生的α 粒子在等离子体中的慢化过程[25].PTC 程序可以利用全轨道模型和漂移轨道模型来追踪粒子运动,在α 粒子慢化过程模拟中,由于全轨道计算较为耗时,而且根据α 粒子的能量和位置估算出α 粒子的回旋半径在0.5—5 cm,回旋半径对慢化过程的影响是可以忽略不计的,因此,可以采用漂移轨道的导心运动方程追踪α 粒子.

在模拟过程中,α 粒子源由单位体积的聚变速率来定义[25],粒子运动方程采用经典的四阶龙格库塔方法求解,粒子与背景等离子的碰撞通过蒙特卡罗方法求解.PTC 程序将慢化过程看作一系列弹性散射,在每个弹性散射中,PTC 先根据α 粒子的初始速度和位置,利用含时演化的轨道方程计算α 粒子的运动轨迹和所处的空间位置;然后,根据该位置区域内的背景电子和离子的密度、温度,计算出散射截面;接下来,使用概率分布函数来模拟碰撞过程中的能量和动量转移.同时,在每个小时间段内会随机采样一定数量的α 粒子进行统计,记录α 粒子的位置、能量等信息,通过大量的随机采样和弹性碰撞模拟,就可以得到α 粒子慢化过程中的分布函数.相较于经典慢化分布理论模型,PTC 程序在轨道运动的基础上,考虑了粒子的新经典输运过程,包括考虑径向扩散和粒子漂移轨道宽度等因素,从而能够更全面、准确地模拟α 粒子的慢化过程.

3 模拟结果

3.1 α 粒子慢化过程

首先,在PTC 程序的输入模块中输入CFETR稳态运行模式下的平衡磁场和热等离子体剖面、扰动电磁场以及α 粒子初始分布;然后在极向截面划分三角形网格,将背景等离子体信息插值到网格,载入α 粒子源;接着,使用漂移轨道模型追踪粒子,使用蒙特卡罗方法求解粒子与背景等离子体的弹性碰撞;在模拟过程中,程序会统计α 粒子数量、计算每个α 粒子的能量损失率、根据碰撞截面计算α 粒子对背景的加热功率.同样利用CFETR混杂运行模式下的平衡磁场、等离子体剖面等参数模拟α 粒子的慢化过程,得到两种模式下慢化过程中的物理量变化对比如图3 所示.

图3 CFETR 稳态运行模式(实线)和混杂运行模式(虚线)下的各个物理量随时间的变化 (a) α 粒子数量;(b) α 粒子损失率;(c) α 粒子对背景等离子体的加热功率;(d) α 粒子平均能量Fig.3.Time evolution of various physical quantities in CFETR steady-state scenario (solid lines) and hybrid scenario (dashed lines): (a) Number of α particles;(b) loss rate of α particles;(c) heating power of α particles to the background plasma;(d) average energy of α particles.

由图3 可以看到,在两种运行模式下,被统计的在等离子体中的α 粒子数量、热化的等离子体数量、损失的α 粒子数量变化趋势几乎一致,分别在0.68 s 和0.73 s 时,稳态运行模式和混杂运行模式下在等离子体中的α 粒子数量趋于平稳,达到稳态时的α 粒子数量在稳态运行模式和混杂运行模式下分别为 1.795×1020和1.709×1020,稳态运行模式下的数量略高于混杂运行模式.稳态运行模式下α 粒子损失率最后维持在 0.41%,混杂运行模式下α 粒子损失率最后维持在 0.38%.稳态运行模式和混杂运行模式下对背景等离子体的加热功率分别在0.57 s 和0.59 s 时开始保持不变,稳态运行模式下α 粒子对背景电子的加热功率达到110 MW,对背景氘离子和氚离子的加热功率分别为45.7 MW 和31.1 MW,整体略高于混杂运行模式结果.稳态运行模式和混杂运行模式下α 粒子的平均能量大约经过0.64 s 和0.72 s 由初始的3.5 MW 趋于稳定的1.61 MW,即α 粒子经过慢化过程达到稳态.

此外,也在归一化极向磁通ψ空间对α 粒子的加热功率密度进行了模拟计算,结果如图4 所示,根据图4 发现,在两种运行模式下α 粒子对背景电子和背景离子的加热功率密度量级一致,在不同的ψ位置有所不同.根据粒子追踪结果也得到了在两种运行模式下达到稳态时α 粒子密度分布和在能量空间、极向磁通空间的分布,如图5 和图6所示.

图4 ψ 空间的加热功率密度Fig.4.Heating power density in the ψ space.

图5 稳态时α 粒子的密度分布 (a) CFETR 稳态运行模式;(b) CFETR 混杂运行模式Fig.5.The α particle density in steady-state: (a) CFETR steady-state scenario;(b) CFETR Hybrid scenario.

图6 α 粒子分布函数 (a) 能量空间;(b) 归一化极向磁通空间Fig.6.The α particle distribution function: (a) Energy space;(b) normalized poloidal magnetic flux space.

3.2 约束时间的计算

根据模拟结果,由粒子平衡方程和能量平衡方程可以得到程序模拟得到总体的α 粒子约束时间和能量约束时间,α 粒子源项每秒钟产生的总粒子数在稳态运行模式和混杂运行模式下分别为3.5×1020和3.24×1020,稳态时在等离子体中的粒子总数在稳态运行模式和混杂运行模式下分别为 1.795×1020和1.709×1020,计算得到α 粒子总体粒子约束时间分别为 0.51 s和0.52 s.同样,根据α 粒子初始总能量和稳态时总能量计算得到α 粒子的总体能量约束时间在稳态运行模式和混杂运行模式下分别为 0.24 s和0.22 s.

4 讨论

4.1 慢化分布对比分析

当α 粒子在非均匀磁场中运动时,粒子在轨道上的漂移会导致粒子的输运方向发生变化,对于高能量α 粒子,由轨道效应引起的输运是可以忽略的,而对于中等能量下的α 粒子,轨道效应引起的α 粒子速度扩散较强,输运不可忽略.如图7 所示,在中等能量E=0.2—2 MeV 附近,α 粒子输运较强,α 粒子在径向的扩散较为明显,PTC 模拟结果较为平缓.在高能量尾部(＞3.5 MeV),在模拟中考虑了库仑碰撞的热化效应,存在一定的能量扩散[26].在不同的空间位置,由于磁场的不同,轨道效应对速度扩散的影响也不同,导致在不同位置处的局部分布函数与理论分布函数符合程度不一致,在本文中,不同空间位置由不同的归一化极向磁通(ψ)表示,ψa和ψb为稳态运行参数下的两个不同位置,ψc和ψd为混杂运行参数下的两个不同位置.从图7 可以看到,PTC 模拟结果与经典慢化分布在稳态运行模式ψb=0.5—0.6 的符合程度高于ψa=0.1—0.2,混杂运行模式下ψd=0.4—0.5 的符合程度高于ψc=0—0.1.

图7 PTC 程序得到的能量慢化分布与理论能量慢化分布的对比 (a) 稳态运行模式下 ψa=0.1—0.2和ψb=0.5—0.6 ;(b) 混杂运行模式下 ψc=0—0.1和ψd=0.4—0.5Fig.7.Comparison between the energy slowing-down distribution obtained by PTCcode and the classical energy slowing-down distribution:(a) In steady-state scenario at ψa=0.1 −0.2 and ψb=0.5 −0.6 ;(b)in hybrid scenario at ψc=0 −0.1 and ψd=0.4−0.5.

为了进一步验证分布函数的不同是由输运效应引起的,分析并对比了Wilkie[27]使用回旋动理学方法考虑存在径向输运得到的修正慢化分布函数.其中,利用高能量下以径向扩散项占主导的扩散系数Drr,得到的动理学方程为

经过修正的三维的稳态慢化分布函数为

将其乘以 4πv2并根据能量E与速度v之间的变换关系以及阶跃函数的定义,得到修正的一维的α 粒子能量慢化分布:

在其他一些研究中,也发现中等能量下的径向输运是不可忽略的.Hauff等[22]的研究发现,经典慢化分布中的局域慢化假设在中等能量下失效,需要考虑轨道效应引起的径向输运,通过垂直去相关机制的研究发现,漂移轨道和拉莫尔轨道导致粒子通量是个关于能量E的函数,粒子扩散系数随着能量增大以1/E减小,中等能量下,粒子扩散系数很大,径向输运很强.Wilkie等[21]也发现经典慢化分布中的通量面局域碰撞导致慢化的假设不成立,碰撞和输运相关的时间尺度会相互竞争,在中等能量下,粒子径向输运时间小于碰撞特征时间,因此经典慢化分布不再适用.同样,在JET 的D-T 实验中[28],也观测到中等能量下的径向通量较大,输运较强.因此准确的α 粒子慢化分布函数,需要考虑径向扩散,更真实的模拟结果有助于更好地预测α 粒子的输运行为和能量沉积分布.

4.2 约束时间的讨论分析

根据经典慢化分布,可得到理论的α 粒子的粒子约束时间和能量约束时间的表达式.选取稳态运行模式下ψ=0.1—0.2 内的平均电子温度和密度,得到理论α 粒子慢化时间为 0.93 s,粒子约束时间和能量约束时间分别为 0.76 s和0.31 s.ψ=0.5—0.6 内α 粒子慢化时间为 0.42 s,粒子约束时间和能量约束时间分别为 0.48 s和0.17 s.对比PTC稳态运行模式下的模拟结果,α 粒子总体慢化时间约为 0.68 s,总体粒子约束时间为 0.51 s,能量约束时间为 0.24 s,理论计算结果与程序模拟结果量级一致.

同样局部选取混杂运行模式下ψ=0.1—0.2 内的平均电子温度和密度,计算得到理论α 粒子慢化时间为 1.01 s,粒子约束时间和能量约束时间分别为 0.79 s和0.33 s,ψ=0.5—0.6 内α 粒子理论慢化时间为 0.43 s,粒子约束时间和能量约束时间分别为 0.48 s和0.17 s.稳态运行模式下芯部的电子温度和密度较为平缓,混杂运行模式下芯部的电子温度和密度变化较快,使得芯部区域的约束性能有所不同,此外,稳态运行模式下的聚变功率略大于混杂运行模式,单位时间内产生的α 粒子更多,慢化时间也有所不同.

5 结论

本研究在CFETR 参数下用PTC 程序模拟研究了聚变产物α 粒子的慢化过程,对比分析了α粒子在慢化过程中不同运行模式下对背景等离子体的加热功率、粒子数变化、约束时间、能量变化的区别,并得到了考虑轨道效应和输运之后更准确的α 粒子慢化分布函数.经过PTC 模拟结果与经典慢化分布的对比发现,在中等能量下二者分布的峰值和能量变化存在不一致,径向输运会抹平中等能量下的能量慢化分布,使慢化分布峰值降低.本文进一步对比了考虑输运之后的修正慢化分布,发现考虑径向扩散之后,修正慢化分布只在中等能量下与慢化分布存在差异,且修正结果与PTC 模拟结果符合程度较好.研究表明,经典慢化分布在中等能量下具有较大误差,PTC 程序考虑了轨道效应与输运过程,可以得到更为准确的α 粒子分布函数,尤其在中等能量下的分布更符合真实的α 粒子的行为.本文的研究结果对于准确预测α 粒子沉积分布和评估α 粒子加热背景等离子体的能力具有重要参考价值.