托卡马克中高能量粒子与内扭曲模相互作用的数值模拟研究 申伟 2015年8月28日

提纲 背景介绍 M3D/M3D-K代码模型 模拟结果: 1. 托卡马克中的锯齿振荡及其导致的高能量粒子输运的模拟研究 2. DIII-D装置中n = 1内扭曲模的模拟研究 总结

内扭曲模 内扭曲模不稳定性（internal kink instability）是托卡马克中常见的磁流体不稳定性，其极向模数为m = 1，环向模数为n = 1。内扭曲模主要发生在安全因子q = 1 有理面之内，其位移函数为： 电阻内扭曲模（resistive internal kink mode）：理想临界不稳定时考虑电阻对内扭曲模的影响，其增长率： γ ~S-1/3~η1/3 Lundquist数S≡τη/τA τη ≡ μ0a2/ η, τA ≡R0/vA John Wesson. Tokamaks, volume 149. Oxford University Press, 2011.

锯齿振荡Ⅰ 锯齿振荡（sawtooth oscillation）是发生在托卡马克等离子体近轴区域的周期性内部磁场重联过程，其模数主要为m = n = 1。 锯齿混合半径（sawtooth mixing radius）一般在q = 1有理面的外部，而锯齿反转半径（sawtooth inversion radius）基本上对应于q = 1有理面。 等离子体温度与安全因子剖面在锯齿崩塌（sawtooth crash）期间快速变化，这是由于锯齿崩塌对应于每个锯齿振荡最后的快速磁场重联过程，可以使等离子体中心区域与周边区域迅速混合。 Plasma Control ITER Physics Expert Group on Disruptions, MHD, ITER Physics Basis Editors. Chapter 3: Mhd stability, operational limits and disruptions. Nuclear Fusion, 1999. 39(12):2251.

锯齿振荡Ⅱ 锯齿振荡的名称来源于由软x射线（soft-x-ray，SXR）信号观测到的锯齿形中心等离子体温度变化。 在锯齿反转半径之外的温度与安全因子的变化方向与此半径之内的相反，说明了在锯齿崩塌发生时，q = 1之内的能量与电流被快速地输运到锯齿混合半径之内的区域中q > 1的部分。 Plasma Control ITER Physics Expert Group on Disruptions, MHD, ITER Physics Basis Editors. Chapter 3: Mhd stability, operational limits and disruptions. Nuclear Fusion, 1999. 39(12):2251.

高能量粒子 高能量粒子（energetic particles，简称为EP）是指其能量远高于热等离子体温度（5 − 10 KeV ），包括聚变反应产生的阿尔法粒子（氦离子，能量为3.52 MeV ）、辅助加热（中性束注入、离子回旋加热、低杂波加热或者电子回旋加热）产生的高能量离子或电子（几十到几百KeV）。 高能量离子通过与背景热电子和热离子碰撞后形成所谓的慢化分布（slowing down distribution）： Robert J Goldston, Paul Harding Rutherford. Introduction to plasma physics. CRC Press, 1995.

M3D/M3D-K代码模型 高能量粒子压强张量: 高能量粒子回旋中心分布函数： 利用回旋动理学（gyrokinetic）或者漂移动理学（drift-kinetic）来描述高能量粒子的运动： W. Park, S. E. Parker, M. Chance et al., Phys. Fluids B 4, 2033 (1992). G. Y. Fu, et al., Phys. Plasmas 13, 052517 (2006).

提纲 背景介绍 M3D/M3D-K代码简介 模拟结果: 1. 托卡马克中的锯齿振荡及其导致的高能量粒子输运的模拟研究 2. DIII-D装置中n = 1内扭曲模的模拟研究 总结

研究背景 锯齿振荡可以引起托卡马克中高能量粒子的显著输运。至今为止，锯齿模引起的高能量粒子输运已经在很多实验装置中被观测到。 对于高能量粒子的输运研究，之前有两种方法被用于描述背景锯齿模的演化。 第一种方法利用了理论模型来描述背景锯齿振荡。Kolesnichenko等人[1,2] 用Kadomtsev模型描述锯齿振荡。 第二种方法则基于实验测量来描述背景电磁场的演化。Zhao和White[3]利用实验参数构造了锯齿模结构，并用来模拟高能量粒子在Tokamak Fusion Test Reactor (TFTR)中的输运。Farengo等人[4,5]利用实验数据重建了锯齿模的空间结构和随时演化。 到目前为止，通过锯齿模的自洽模拟来研究高能量粒子输运的工作几乎没有。 我们利用M3D/M3D-K代码来模拟锯齿振荡，进而研究其所导致的高能量粒子的输运。其中背景锯齿振荡是用电阻磁流体模型模拟，而高能量粒子用测试粒子（test particle）方法推动。 [1] Y. I. Kolesnichenko and Y. V. Yakovenko, Nucl. Fusion 36, 159 (1996). [2] Y. I. Kolesnichenko, V. V. Lutsenko, et al., Phys. Plasmas 4, 2544 (1997). [3] Y. Zhao and R. B. White, Phys. Plasmas 4, 1103 (1997). [4] R. Farengo, H. E. Ferrari, et al., Plasma Phys. Controlled Fusion 54, 025007 (2012). [5] R. Farengo, H. E. Ferrari, et al., Nucl. Fusion 53, 043012 (2013).

锯齿振荡模拟结果：平衡剖面以及初始参数 环径比R0/a=4的圆截面托卡马克。 电阻分布为Spitzer形式： η(T)= η0(T/T0)-3/2，其中η0=10-5，相应地 Lundquist数S=1/η0=105。 热等离子体中心比压值为β0 =4.38%。 初始平衡是由平衡代码VMEC得出。 安全因子：q0=0.7, q1=3.6 ψ是归一化的极向磁通（在磁轴处ψ =0，在边界处ψ =1）。 高能量粒子分布： 其中Λ ≡ μB0/E是抛射角参数，⟨Ψ⟩是粒子轨道平均后的ψ值。

线性磁流体模拟结果 模数为n=m=1的线性不稳定的内扭曲模。 U 线性磁流体模拟结果 模数为n=m=1的线性不稳定的内扭曲模。 该模主要位于q=1共振面之内，线性增长率为γ τA=0.0283，其中τA是阿尔芬时间，τA ≡R0/vA，vA≡B0/(μ0ρ0)1/2。 线性增长率正比于S-1/3。根据定标律此模是电阻内扭曲模。

非线性磁流体模拟结果 在非线性模拟之前，电流和压强源项需加入。 图 (a)显示了总动能与中心温度T(0)的随时变化，说明周期性的锯齿振荡被成功地模拟出来。 如图(b)所示，保留了从n=0到n=6的环向模数，不同的模耦合在一起，其中n=1的模是最重要的。

一个锯齿周期内的磁面随时演化

高能量粒子输运模拟结果: 理论背景 根据Kolesnichenko等人[1,2]给出的物理图像，假设△ rb <<rmix，其中△ rb代表高能量粒子轨道宽度，而rmix代表锯齿的混合半径。两个特征时间决定了高能量粒子在锯齿模中的输运行为： （1）环向进动时间τpr，它与粒子环向漂移运动相关； （2）锯齿崩塌时间τcr 。 如果锯齿崩塌足够快的话（τcr << τpr ），高能量粒子被冻结在磁力线上并随着变化的磁面一起运动，但是高能量粒子的输运会被其环向进动漂移减弱，这是由于粒子的环向进动会导致粒子运动和锯齿模之间解耦，而环向进动的作用会被粒子纵向运动减弱。 τpr = τcr给出了捕获粒子的临界能量: Ecrit = 2πMκsrsR0ωB/τcr 其中M是离子质量，κ是椭圆率，ωB是回旋频率，“s”代表在q = 1有理面上的值。 [1] Y. I. Kolesnichenko and Y. V. Yakovenko, Nucl. Fusion 36, 159 (1996). [2] Y. I. Kolesnichenko, V. V. Lutsenko, et al., Phys. Plasmas 4, 2544 (1997).

高能量粒子输运模拟结果总述 高能量粒子再分布程度的定义是： 其中 代表了粒子的径向位置，能量由Ecrit归一。 Ecrit = 2πMκsrsR0ωB/τcr对应于：vh/vA=0.7275, ρh/a =0.04656 下图显示了高能量粒子的再分布程度依赖于粒子的抛射角和能量。对于通行粒子，Λ=0.0；对于捕获粒子， Λ=1.0。

捕获粒子结果分析Ⅰ averaged ψ≡(ψmax+ψmin)/2，其中ψmax和ψmin分别为粒子在一个反弹周期（bounce period）内所经过的极向磁通的最大值和最小值。 当捕获粒子的能量高于0.8 Ecrit时，在q = 1面内的粒子有很强的环向进动漂移，因此它们几乎不被输运。 捕获粒子的进动频率。虚线（黑线）代表q = 1面的位置，而实线（白线）代表ωpr =2π/τcr。实线和虚线相交于0.8 Ecrit附近。

捕获粒子结果分析Ⅱ Ecrit=2πMksrsRoωB/ τcr 图(a)中Ecrit用其对应的不同锯齿模崩塌时间计算。模拟中临界能量 Esim用|∂g/∂E| < 0.02的第一个点定义。如图(b)所示， Esim大致反比于锯齿崩塌时间。 如下图所示，捕获粒子的临界能量和粒子质量无关。

通行粒子结果分析Ⅰ 当通行粒子的轨道宽度很小时（Δrb ≪ r，其中r是径向坐标），并且等离子体横截面是圆形的，通行粒子的进动频率为[1]： 其中： ωB0 = eB0/M，B0 = B(r = 0)，ρ = v/ωB0，Δ是Shafranov位移, Δ′ = dΔ/dr，Δ′′ =d2Δ/dr2，q′ = dq/dr， 是磁剪切，σv ≡ sign(v∥)， 以及 左图比较了由以上方程计算出的通行粒子进动频率与我们模拟中的对应频率，我们发现在q = 1有理面之内两者符合的很好。 [1] Y. I. Kolesnichenko, R. B. White, and Y. V. Yakovenko, Phys. Plasmas 10, 1449 (2003).

通行粒子结果分析Ⅱ 在锯齿模区域内的通行粒子进动频率如图（a）和（b）所示。 对于同向通行粒子，在此区域内的大部分位置ωpr < 2π/τcr，所以同向通行粒子被很强的再分布。 相反地，反向通行粒子的进动频率比同向通行粒子大得多。当反向通行粒子的能量增加时，ωpr的值将会高于2π/τcr，相应地反向通行粒子的再分布会减弱。 实线（白线）代表了ωpr = 2π/τcr。

提纲 背景介绍 M3D/M3D-K代码简介 模拟结果: 1. 托卡马克中的锯齿振荡及其导致的高能量粒子输运的模拟研究 2. DIII-D装置中n = 1内扭曲模的模拟研究 总结

DIII-D实验结果简介 在中性束加热的DIII-D等离子体放电中，如图所示，在q = 1面附近的局部电子回旋辐射（electron cyclotron emission，简写为ECE）诊断所得的频谱包含了连续锯齿崩塌（竖实线）之间的多个m = 1, n = 1振荡，其中包括频率大致从15.5 kHz变为14 kHz的类鱼骨模（A）,与频率固定在14 kHz附近的振荡（B）。 在t = 3000 ms时，除了在轴中性束注入，离轴中性束开始注入。由于束粒子分布随着时间越变越宽，（A）逐步演化为（B）。 我们模拟所选的等离子体平衡剖面是在t = 3110 ms (竖虚线)处。

初始剖面及相关参数设定 DIII-D放电的主要参数和剖面如下： 大半径(major radius)：R0=1.63 m 小半径(minor radius)：a=0.616 m 拉长比(elongation)：κ=1.65 三角形变(triangularity)：δ=0.036 磁轴处环向磁场：B0=1.90 T 中心密度：n0=5.6×1019 m-3 热等离子体中心压强：βthermal,0=3.40% 束离子中心压强：βhot,0=1.67% 安全因子： q(0)=0.911, q(1)=4.533, q(0.191)=1 束离子分布函数： 束离子同时包含了同向通行（co-passing）粒子和捕获（trapped）粒子，并且归一化的粒子速度和回旋半径分别为：

U 磁流体模拟：线性结果 束离子用磁流体模型描述（即忽略束离子动理学效应）。 线性增长率为γτA = 0.0141的理想磁流体不稳定的扭曲模，此模主要位于q = 1有理面之内，并且模数主要为n = m = 1。

磁流体模拟：非线性结果 Lundquist数S = 1/η0 = 105。粘滞和垂直热导在径向均匀分布，分别为: , 左图：中心压强P(0)和不同环向模的动能指出此模在初始锯齿崩塌之后演化为三维准稳态，并且n = 1是最主要的模数。 下图比较了在有和没有高环向模数的情况下n = 1模的动能演化。在两种情况下n = 1模都在锯齿崩塌后饱和，并且高n模数增强了n = 1模动能的饱和程度。

磁面与压强剖面 非线性饱和阶段的磁面几乎随时不变，磁面结构如左图所示。 同一个环向角横截面上的压强剖面如右图所示，压强在q = 1面之内变平，并且其剖面与磁面结构相符。

磁流体非线性模拟结果： 改变电阻与垂直热导 电阻，粘滞以及垂直热导比值固定（即ν/η0 = 3，χ⊥/η0 = 30）。如图(a)所示，当电阻值在η0 = 3.33 ×10−6 ∼ 9 × 10−5范围内时，所有的算例结果都是三维准稳态。 图(b)显示了固定电阻和粘滞值，并取不同χ⊥ 值的算例的动能随时演化。当垂直热导值下降到比临界值低后，扭曲模的非线性演化由准稳态转变为锯齿振荡，与Halpern等人的结果相似[1]。 [1] F. D. Halpern, et al. Plasma Phys. Controlled Fusion 53, 015011 (2011).

包含束离子动理学效应模拟： 线性模结构对比 如图(a)所示，在磁流体极限下模结构为上下对称，且模频率为0。 考虑束离子动理学效应：当离子束压强足够大时，n = 1模具有一定的模频率以及扭曲的模结构，如图(b)和(c)所示。 本征模的速度流函数（Velocity stream function）U。(a) 没有高能量粒子。(b) 有高能量粒子且参数为ΔΨ = 0.39，Pbeam,0/Ptotal,0 = 0.329。(c) 有高能量粒子且参数为ΔΨ = 0.25，Pbeam,0/Ptotal,0 = 0.418。

包含束离子动理学效应的模拟:非线性结果 在相同的热压强Pthermal以及径向积分后的离子束压强 的条件下，我们选择束离子径向分布不同的两个算例来进行非线性模拟。 左图显示了在两个算例中，内扭曲模都在首个锯齿崩塌后达到饱和的三维准稳态，并且n = 1的模是最主要的。 (a) ΔΨ =0.39，Pbeam,0/Ptotal,0 = 0.329。 (b) ΔΨ =0.25，Pbeam,0/Ptotal,0 = 0.418。

动能比较 左图比较了在有和没有高能量束离子的情况下n = 1和n = 0模的动能。 如图(a)所示，束离子在线性阶段对n = 1模起去稳作用，但是在非线性阶段降低了n = 1模的饱和幅度。 如图(b)所示，n = 0模的饱和幅度在有束离子的情况下更高。 束离子更窄的径向分布会导致更大的线性增长率和更高的n = 0模的饱和幅度。

模频率随时演化 下图显示了两种束离子分布下模频率的随时演化。 对于较宽的束离子分布，在非线性阶段的模频率比最初的线性频率低一点。 然而，对于更窄的束离子分布，模频率在非线性阶段下降的更加明显。

内扭曲模的饱和是由于磁流体非线性 只保留n = 1的环向扰动，并且将本底等离子体的磁流体响应限制为线性的。 没有磁流体非线性条件下n = 1模动能的随时演化。实线（红线）对应于ΔΨ =0.39和Pbeam,0/Ptotal,0 = 0.329的算例，而虚线（蓝线）对应于ΔΨ = 0.25 和Pbeam,0/Ptotal,0 =0.418的算例。

非线性阶段的模旋转由极向带状流引起 图(a)和(b)分别显示了非线性饱和阶段环向和极向速度的n = 0分量。极向的n=0速度远比环向的大，因此n=0速度主要在极向。 从图(b)和(c)可以看出，模位置的半径为rmode ∼ 0.2a，在模位置的平均极向带状速度为vp,ZF ∼ 2 × 10−3 (avA/R0)，因此极向带状速度引起的模频率大致为ωZF = vp,ZF /rmode ∼ 0.01 ωA。非线性阶段模频率在0.006 ωA附近，这与极向带状流引起的频率值相符。 在非线性阶段的模旋转是由流体非线性产生的极向带状流引起的。 ΔΨ = 0.39和Pbeam,0/Ptotal,0 = 0.329算例在t = 2000 τA的等高线剖面。(a) n = 0环向速度。(b) n = 0极向速度。(c) ϕ = 0 极向截面内的速度流函数U。

内扭曲模对束离子分布的影响 左图显示由类鱼骨模（即ΔΨ = 0.25和Pbeam,0/Ptotal,0 = 0.418 对应的算例）引起的速度为v = 0.705 vA的束离子再分布。 首先，在初始锯齿崩塌（t = 500 τA）之后，同向通行粒子和捕获粒子都在锯齿模范围内被很强的再分布。 之后，在扭曲模的非线性饱和阶段（从t = 1500 τA到t = 2500 τA），同向通行粒子和捕获粒子的分布都在q = 1面内变得更平。 由ΔΨ = 0.25和Pbeam,0/Ptotal,0 = 0.418对应的类鱼骨模引起的束离子再分布。(a) Λ = 0.2和v/vA = 0.705的同向通行粒子。(b) Λ = 1.0 和v/vA = 0.705的捕获粒子。

总结 使用磁流体-动理学混合模型代码M3D-K对高能量粒子与内扭曲模相互作用进行了数值模拟。 托卡马克中的锯齿振荡及其导致的高能量粒子输运： 锯齿振荡引起高能量粒子在等离子体中心区域径向再分布，并且再分布的程度依赖于其抛射角与能量。 捕获粒子和通行粒子的进动频率起着重要的作用。 DIII-D装置中n = 1内扭曲模： 内扭曲模在锯齿崩塌之后演化为饱和扭曲模。 磁流体：垂直热导降低时，由准稳态转变为锯齿振荡。 考虑束离子动理学效应：n = 1模非线性阶段饱和与旋转的最主要机制是磁流体非线性。