任珍珍 申伟

1) (安徽大学物理与光电工程学院,合肥 230601)

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

负三角形变位形下的托卡马克具有更低的湍流输运和更好的能量约束,被认为是未来聚变堆一个更好的选择.为了探索负三角形变位形下高能量粒子激发不稳定性的特征,使用动理学-磁流体混合模型程序M3D-K 开展了此位形下高能量离子激发鱼骨模的线性不稳定性和非线性演化的模拟研究.基于类EAST 参数条件,模拟发现负三角形变解稳理想内扭曲模不稳定性,但会致稳鱼骨模不稳定性.非线性模拟发现在没有磁流体非线性效应时,负三角形变位形下的鱼骨模更不容易饱和,可能的解释是相比于正三角形变位形,在负三角形变位形下的高能量离子轨道更接近与芯部,因而更容易驱动鱼骨模不稳定性.这些结果表明考虑高能量粒子激发的鱼骨模不稳定性后,负三角形变位形相比于正三角形变位形并没有明显优势.

1 引言

托卡马克是最有前景的磁约束聚变装置之一,其中等离子体截面的形变对于其约束与稳定性具有显著的影响.因此,对于未来的聚变装置,形变参数是关键的设计参数之一[1−4].近年在TCV 托卡马克装置上的研究发现相比于传统的正三角形变,负三角形变可降低湍流进而改善等离子体约束[5,6].随后,DIII-D 托卡马克装置上也采用了负三角形变位形在L 模下放电,发现达到了较高的归一化比压值(βN=2.7)以及相当于H 模的约束性能(H98y2=1.2),并且没有边界局域模发生[7].另外,负三角形变位形可以降低偏滤器的热负载,这些优点使其被认为是未来聚变堆的一个更好的选择[8].

高能量粒子物理对于未来聚变堆中的燃烧等离子体来说是一个很重要的研究方向[9].这是由于在未来的磁约束核聚变堆中,实现自持燃烧的等离子体中的氘氚聚变反应会产生大量的快阿尔法粒子.此外,辅助加热下如中性束注入、离子回旋加热等也会产生大量的快粒子.这些快粒子可以激发出各种不稳定性,这些不稳定性反过来会引起很强的快粒子输运并导致快粒子损失到装置的第一壁,大量的快阿尔法粒子输运会严重降低其对背景等离子体的加热效率.另外,快粒子可能对磁流体不稳定性如新经典撕裂模、内扭曲模、电阻壁模等有显著的影响,而这些磁流体不稳定性有可能严重降低等离子体约束性能甚至会导致托卡马克破裂.

高能量粒子激发的不稳定性包括阿尔芬本征模不稳定性、快粒子模、鱼骨模等.其中,鱼骨模不稳定性最早在PDX 装置上有中性束注入的条件下被观测到[10],是环向模数和极向模数都为 1 的一种芯部不稳定性.鱼骨模不稳定性在实验中被观测到之后,有两种理论被用来解释鱼骨模不稳定性形成的物理机制.1984 年Chen等[11]认为鱼骨模是捕获高能量离子通过进动频率共振激发的.另外,1986 年,Coppi 和Porcelli[12]在考虑热离子逆磁效应后认为鱼骨模是一种类磁流体模,其频率与热离子逆磁漂移频率相当.随后,鱼骨模在多个装置上被观测到[13−20],因此鱼骨模是一种常见的高能量粒子激发的不稳定性.

然而,目前针对负三角形变位形的磁约束等离子体开展的高能量粒子物理相关的工作较少,DIII-D 装置上开展了负三角形变位形下的高能量粒子激发阿尔芬本征模不稳定性的实验,发现相比于正三角形变位形阿尔芬本征模引起的高能量粒子输运并没有明显减弱[21],但是相关的工作没有理论模拟方面的分析,主要是实验观测结果.因此,本工作针对负三角形变位形等离子体,使用动理学-磁流体混合模型程序M3D-K 开展了高能量粒子激发的鱼骨模不稳定性的模拟研究,线性模拟结果发现负三角形变位形解稳理想内扭曲模不稳定性,但会致稳鱼骨模不稳定性.然而,非线性模拟发现在没有磁流体非线性效应时,负三角形变位形下的鱼骨模更不容易饱和.这些结果表明考虑高能量粒子激发的鱼骨模不稳定性后,负三角形变位形相比于正三角形变位形并没有明显优势.

本文分为以下几个部分,第2 节简单介绍了M3D-K 程序所用的物理模型以及模拟所用的主要参数.第3 节展示了负三角形变下鱼骨模不稳定性的线性模拟结果.第4 节展示了负三角形变下鱼骨模不稳定性的非线性模拟结果.最后,第5 节给出了本工作的结论.

2 M3D-K 程序模型以及参数设置

2.1 M3D-K 程序模型

本文使用了动理学-磁流体混合模型程序M3D-K[22,23].该程序使用电阻磁流体模型描述热等离子体,使用漂移动理学或者回旋动理学模型描述高能量粒子.磁流体方程组在极向截面采用有限元方法求解,在大环方向采用4 阶有限差分或者拟谱法.由于M3D-K 代码采用柱坐标系,在磁轴处不存在奇异点,因而适用于等离子体芯部磁流体不稳定性的数值模拟.对于高能离子部分,为了降低粒子噪声以及模拟所需粒子数量,程序采用了δf方法.M3D-K 程序通过在动量方程中包含高能离子压强张量,将高能离子与背景磁流体自洽耦合.对于背景等离子体部分,程序还可以包含双流体效应,以及包括热离子逆磁漂移、电子压强等物理效应.M3D-K 程序已经成功被用来模拟各种磁流体不稳定性与高能离子的相互作用[24−35],是模拟高能离子不稳定性的有力工具.

M3D-K 程序所用的方程组如下,高能量粒子效应通过压强耦合进入动量方程:

连续性方程:

能量方程:

法拉第定律:

安培定律:

欧姆定律:

这里ρ,P,v,E,B,J,γ,µ0,κ,η,ν分别表示质量密度、压强、流体速度、电场、磁感应强度、电流密度、绝热系数、真空磁导率、热传导系数张量、电阻率、黏滞系数.d/dt=∂/∂t+v·∇.

高能量粒子压强张量在忽略非对角项的情况下使用Chew-Goldberger-Low (CGL)形式:

式中,I表示单位矩阵;b是磁场方向单位矢量;平行和垂直磁场方向的压强利用高能量粒子分布函数F在回旋中心坐标 (X,v//,µ) 中计算,

其中ρh=v⊥×b/Ω是回旋半径矢量,Ω ≡eB/m是回旋频率,e是电子电荷量,(x,v) 是高能量粒子相空间坐标,µ是磁矩,θ是回旋角.

高能量粒子回旋中心分布函数F可以表示为

其中δ是狄拉克函数.

高能量粒子的运动用回旋动理学或者漂移动理学来描述:

其中,E是总的电场,B是总的磁场,b0是沿平衡磁场方向的单位矢量,B0和δB分别是平衡和扰动的磁场.〈·〉 表示回旋平均,如果是漂移动理学的话直接用粒子导心处的电磁场.变量B＊和B＊＊分别表示为

2.2 初始参数与平衡剖面设置

主要的平衡剖面和参数根据EAST 托卡马克装置的#71320 炮设置[36],主要参数如下: 小半径a=0.45 m,大半径R0=1.86 m,椭圆形变κ=1.63,三角形 变δ=-0.436,磁轴处磁场B0=1.615 T,磁轴处电子密度ne0=5.52×1019m-3,阿尔芬速度vA=B0/(µ0ρ0)1/2=3.35×106m/s,阿尔芬时间τA=R0/vA=5.55×10-7s,磁轴处包括热等离子体和快离子的总压强βtotal,0=4.61%.这里µ0是真空磁导率,ρ0是磁轴处的热等离子体质量密度.平衡的安全因子剖面和压强剖面如图 1 中蓝线所示,其中Ψ是归一化极向磁通,其在磁轴处值为 0,在等离子体最外闭合磁面处为 1.q=1 面的位置如 图 1 中红色虚线所 示,在Ψ=0.171 处.在本工作中为了简化分析,热等离子体密度的径向分布取为常数.

图1 安全因子与总压强平衡剖面Fig.1.Equilibrium profiles of safety factor and total pressure.

考虑中性束注入产生的高能量离子,其入射能量为E0=58 keV.高能量离子在速度方向以及在实空间的径向分布表达式为

3 线性模拟结果

首先分析在负三角形变条件下的鱼骨模的线性不稳定性结果.图2 对比了正三角形变δ=0.436,无三角形变δ=0,负三角形变δ=-0.436位形下快离子激发内扭曲模以及鱼骨模的情况,图中Phot,0/Ptotal,0是快离子压强与总压强比值,Phot,0是磁轴处快离子压强值,Ptotal,0是磁轴处总压强值.可以看出,当没有快离子效应时,三种位形下的理想内扭曲模都是不稳定的,并且负三角形变对理想内扭曲模不稳定性起解稳作用,正三角形变对理想内扭曲模起致稳作用.当快离子压强比值Phot,0/Ptotal,0增大时,快离子的动理学效应会致稳内扭曲模[37−39].然而,当快离子压强增大到超过一定阈值时,鱼骨模不稳定性被快离子共振激发.另外,从图 2 可以看出,无三角形变或者负三角形变下鱼骨模很难激发(激发阈值值比较高),但一旦被激发则增长更快;正三角形变下激发鱼骨模的阈值比较低,但被激发后却增长得慢一些.对于三角形变影响内扭曲模的分析,之前Eriksson 和Wahlberg[40]的工作发现正的椭圆形变和正的三角形变对理想内扭曲模起稳定作用.Martynov等[41]发现三角形变对理想内扭曲模的势能 δW有重要的贡献,考虑到三角形变后解析结果与KINX 程序结果一致.另外,他们使用KINX 磁流体本征程序模拟发现一定参数条件下负三角形变对理想内扭曲模起解稳作用.这些理论分析结果与本文模拟的结果一致.

图2 不同三角形变下模频率和线性增长率与快离子压强比值 Phot,0/Ptotal,0 的关系Fig.2.Mode frequency and linear growth rate as a function of the fast ion pressure fraction Phot,0/Ptotal,0.

图3 给出了不同高能量离子压强比值的模结构.在负三角形变位形下,当快离子压强比值Phot,0/Ptotal,0=0 时,理想内扭曲模频率为 0,且模结构是上下对称的(图 3(a));当快离子压强比值Phot,0/Ptotal,0=0.15时,模结构变得扭曲并且有一定大小的模频率(图 3(b));当快离子压强比值Phot,0/Ptotal,0=0.4时,模式变为鱼骨模不稳定性,并且有更高的模频率(图 3(c)).作为对比,图 3(d)给出了正三角形变位形下的鱼骨模,可以看出模结构与负三角形变位形下的鱼骨模很相似.

图3 不同高能量离子压强比值 Phot,0/Ptotal,0 下的流函数U(a) δ=-0.436,Phot,0/Ptotal,0=0 ;(b) δ=-0.436,Phot,0/Ptotal,0=0.15 ;(c) δ=-0.436,Phot,0/Ptotal,0=0.4 ;(d) δ=0.436,Phot,0/Ptotal,0=0.4Fig.3.Velocity stream function U at different fast ion pressure fraction Phot,0/Ptotal,0 : (a) δ=-0.436,Phot,0/Ptotal,0=0 ;(b) δ=-0.436,Phot,0/Ptotal,0=0.15 ;(c) δ=-0.436,Phot,0/Ptotal,0=0.4 ;(d) δ=0.436,Phot,0/Ptotal,0=0.4.

图4 芯部压强剖面平坦下模频率和线性增长率与快离子压强比值 Phot,0/Ptotal,0 的关系Fig.4.Mode frequency and linear growth rate as a function of the fast ion pressure fraction Phot,0/Ptotal,0 with flat pressure profile.

另外,由于未来聚变堆运行时需要先稳定磁流体模式,选取了没有高能量粒子影响下理想内扭曲模稳定的算例开展研究.考虑到环效应对内扭曲模稳定性的影响时,当q=1 面之内的压强梯度足够小时,内扭曲模会趋于稳定[42].因此,选取了一个芯部压强剖面平坦并且磁轴处压强值保持不变的算例,重新计算了不同快粒子驱动时3 个三角形变参数下的内扭曲模的线性稳定性,压强剖面如图 1的红线所示.新压强剖面下的鱼骨模增长率以及模频率如图 4 所示.可以看出,在新压强剖面下鱼骨模的激发阈值更高,并且正三角形变下鱼骨模的激发阈值仍然低于负三角形变.

4 非线性模拟结果

本部分考虑没有磁流体非线性效应下的鱼骨模非线性演化,这里将热等离子体的磁流体响应通过只保留n=1 环向扰动而限制为线性的.由于线性部分 中Phot,0/Ptotal,0=0.4 条件下负三角形变位形的鱼骨模比正三角形变位形的鱼骨模增长率高.为了更好地对比,选取了一个更低总比压值βtotal,0=3.91%的负三角形变位形下的鱼骨模算例,并且快离子压强比值仍保持为Phot,0/Ptotal,0=0.4,如图 5 所示.可以看出,图5(a) 中正三角形变鱼骨模与负三角形变鱼骨模的线性增长率基本一样,增长率均为γτA=0.00906.但是正三角形变的鱼骨模可以非线性饱和,而负三角形变的鱼骨模不能饱和.图5(b) 给出了两种三角形变位形下鱼骨模向下扫频的特征,正三角形变位形下鱼骨模频率从ω=0.0413ωA向下扫频到ω=0.00707ωA,负三角形变位形下鱼骨模频率从ω=0.0336ωA向下扫频到 0.图6为Phot,0/Ptotal,0=0.4,βtotal,0=4.61%,δ=-0.436 算例的动能和磁能的n=1 分量随时间演化.可以看出,它们都随时间线性增长随后饱和,线性增长率一致.但是,由于本文的非线性模拟是考虑没有磁流体非线性效应下的演化,在演化中磁能的n=0 分量被限制为不变而只演化n=1 分量.

图5 没有磁流体非线性效应的鱼骨模非线性演化 (a) 动能的 n=1 分量的 演化;(b) 模频率演化Fig.5.Time evolution of the fishbone without fluid nonlinearity: (a) n=1 kinetic energy;(b) mode frequency.

图6 动能和磁能的 n=1 分量随时间演化Fig.6.Time evolutions of n=1 kinetic energy and magnetic energy.

为了进一步确认三角形变对于鱼骨模非线性演化的影响,分析了不同三角形变参数下鱼骨模的非线性演化特征.图7(a)给出了在Phot,0/Ptotal,0=0.4 以及βtotal,0=4.61% 条件下三角形变参数 从δ=0 增大到δ=0.436 时鱼骨模的非线性演化.可以看出,当没有三角形变时,鱼骨模非线性并不能饱和;当三角形变从δ=0 增大到δ=0.3 时,鱼骨模的线性增长率逐步降低,但是鱼骨模在非线性阶段仍然不能饱和;当三角形变增大到δ=0.436时,鱼骨模的增长率进一步减小,其非线性演化到饱和状态.作为对比,图7(b)给出了Phot,0/Ptotal,0=0.4 以及βtotal,0=3.91% 条件下三角形变参数从δ=0 减小到δ=-0.436 时鱼骨模的非线性演化.当三角形变δ=0 时,鱼骨模的非线性演化到饱和状态,与图7(a)中δ=0 的算例对比说明总比压值降低有利于鱼骨模非线性演化趋于饱和.当三角形变从δ=0 减小到δ=-0.2 时,鱼骨模的非线性演化仍然是饱和状态,当三角形变进一步减小至δ=-0.3以下时,鱼骨模的非线性演化趋于无法饱和.值得注意的是,三角形变参数从δ=-0.2 降低到δ=-0.436 时,鱼骨模的线性增长率增幅不大,但是非线性演化却变化非常明显.

图7 不同三角形变参数下的鱼骨模动能的 n=1 分量的演化(a) βtotal,0=4.61% ;(b)βtotal,0=3.91%Fig.7.The n=1 kinetic energy evolution of the fishbone with different triangularity: (a) βtotal,0=4.61% ;(b)βtotal,0=3.91%.

为了解释负三角形变参数下鱼骨模非线性演化不容易饱和的原因,分别在图 5 的两个算例中取了两组快离子对比,每组快离子的能量、抛射角以及所在的平均径向位置都基本相同,如图8所示.图中红色的粒子轨道对应粒子能量E=57.7 keV,Λ=1.05,Pφ=-0.27,其中Pφ是快离子的环向角动量,对应于快离子在径向的位置,Pφ=ZeΨ+mv//RBφ/B.另外,黑色的粒子轨道对应粒子能量E=49.8 keV,Λ=0.985,Pφ=-0.33.可以看出,两组高能量离子都是捕获粒子且位于等离子体芯部,与鱼骨模结构的位置基本相同.研究发现,在负三角形变位形下的高能量离子轨道更接近芯部,因此更有利于快离子与鱼骨模之间的能量交换,这可能是在负三角形变位形下鱼骨模更不容易饱和的原因.

图8 不同三角形变位形下的捕获高能量离子轨道 (a)δ=-0.436,βtotal,0=3.91% ;(b) δ=0.436,βtotal,0=4.61%Fig.8.Orbits of trapped fast ions with different triangularity: (a) δ=-0.436,βtotal,0=3.91% ;(b) δ=0.436,βtotal,0=4.61%.

5 结论

本文使用动理学-磁流体混合模型程序M3D-K研究了负三角形变位形下高能量离子激发鱼骨模的线性特征以及非线性演化.基于类EAST 参数条件,线性模拟结果发现负三角形变对内扭曲模起解稳作用,正三角形变对内扭曲模起致稳作用.当快离子压强增大到超过一定阈值时,不同三角形变位形下的鱼骨模不稳定性都会被快离子共振激发,并且正三角位形下鱼骨模的快离子压强激发阈值更低.另外,正三角形变会致稳鱼骨模不稳定性,而负三角形变对鱼骨模起轻微的致稳作用.非线性模拟研究结果发现降低总比压值有利于鱼骨模非线性演化趋于饱和.另外,发现负三角形变下鱼骨模非线性演化不容易饱和,其可能的原因在于负三角形变位形下的高能量离子轨道更接近芯部,因此更有利于高能量离子驱动鱼骨模时将能量传输给鱼骨模不稳定性.

本研究没有考虑热等离子的磁流体非线性,因而在鱼骨模的非线性演化中没有考虑模模耦合以及n=0 带状流,这些效应应该有助于鱼骨模的非线性演化趋于饱和,这将作为以后工作考虑的物理效应.另外,本文采用的是类EAST 参数条件,未来可以对ITER 装置下做进一步的负三角形变鱼骨模不稳定性的理论模拟分析.

感谢浙江大学傅国勇教授的讨论.本文的数值计算是在中国科学院等离子体物理研究所的“神马”大型高性能计算集群上开展的.