Effects of magnetic field on electron power absorption in helicon fluid simulation

磁场对螺旋流体模拟中电子功率吸收的影响

摘要

本研究开发了一种名为 “北京大学螺旋放电”（PHD）的代码，可以模拟背景磁场大于 500 G 和压力小于 1 Pa 的螺旋放电过程。代码中使用了两个流体方程。PHD 模拟得出了两个重要发现：（1）等离子体密度随背景磁场的时间演化呈现出第二次快速增长（称为第二次密度跃迁），类似于螺旋子等离子体中的模式转换；（2）在存在磁场的情况下，电子功率吸收的峰值位置出现在中心轴附近，这与没有磁场的情况不同。这些结果可能有助于加深对放电机制的理解。

介绍

螺旋等离子体具有高等离子体密度、低压运行、稳态放电、高功率耦合效率以及等离子体区域无需电极等理想特性。除了用于等离子体研究以及可变比冲磁等离子体火箭、氦质子等离子体推进器和氦质子双层推进器等应用[1-7]，它们还因其多种优势而被用于材料加工[1, 8]。

关于螺旋子等离子体放电机制的统一理论尚未建立[9]。氦质子等离子体主要表现出电磁和静电两种耦合模式。电磁模式为弱阻尼，称为氦子波，而静电模式为强阻尼，称为 Trivelpiece-Gould (T-G) 波 [10-12]。陈指出，朗道阻尼机制在螺旋子放电中发挥了重要作用 [13]；这一点也得到了实验验证 [14]。进一步的实验研究表明，快速电子的数量不足以提供完全电离，这证明朗道阻尼并不重要[15, 16]。一些关于功率沉积的模拟，如基于 ANTENA2 和 HELIC 的模拟表明，T-G 波是电离机制的重要组成部分 [11，17]。此外，实验也验证了 T-G 波的存在 [18]。目前，人们普遍认为螺旋波和 T-G 波可能是导致电离的主要因素[1]。

表 1.一些螺旋仪及其典型运行参数。

数值模拟是研究谐振子等离子体的有效方法。特别是，通过模拟可以很容易地获得空间和时间上每一点的所有物理量；相比之下，这些物理量通常很难通过实验获得。完整的模拟包括功率沉积、传输过程和电离过程 [19]。Bose 等人建立了一个详细的二维轴对称流体模型，并结合麦克斯韦方程对波的传播、等离子体传输和电离过程进行了数值模拟。研究发现，随着磁场强度的增加，天线激发的波会向等离子体更深处传播[8]。Yang 等人提出了一种利用三维流体力学方程的螺旋子等离子体放电直接数值模拟模型。他们研究了不同天线电流下的放电特性以及螺旋子放电中低场峰的出现 [19，20]。Jaafarian 等人利用粒子在小室中的蒙特卡罗碰撞模拟技术研究了工作参数对螺旋子放电的影响[21]。CST [10]、HELIC [17]、MAXEB [11]、ANTENA2 [22]、SPIREs [23] 和 ADAMANT [24] 被用来研究电离机制。

如表 1 所示，许多螺旋子装置在低压（∼0.5 Pa）、强磁场（∼1kG）和高功率（∼kW）条件下工作。迄今为止，数值模拟研究只关注波的传播和功率沉积，而电离和传输过程通常被忽视[10, 17, 22-24]。据我们所知，在背景磁场超过 500 G 和压力小于 1Pa 的同一时间内，包括功率沉积、传输过程和电离过程的模拟尚未见报道，尽管这些参数可以通过实验实现。因此，我们开发了北京大学螺旋放电（PHD）代码来研究这些参数下的螺旋子放电过程。PHD 模型在平行磁场方向使用了伏极扩散近似，大大提高了计算效率。此外，我们还首次发现了螺旋子放电过程模拟中的二次密度跃迁。

名古屋 III 型天线、半螺旋天线、博斯韦尔天线和单环天线等多种类型的天线都可以产生太阳等离子体[1, 32, 33]。前三种天线主要激发 =m±1 极模式，而单环天线主要激发 =m 0 模式。我们结合北京大学等离子体试验（PPT）装置[31]的实验研究，对单环天线产生的 = m 0 谐子等离子体的放电过程进行了研究。图 1 是 PPT 实验装置的示意图。PPT 实验使用亥姆霍兹线圈产生背景磁场，而 PHD 模拟使用均匀背景磁场。我们将流体方程和泊松方程与麦克斯韦方程相结合，对螺旋子放电过程（包括电子功率吸收、电离过程和传输过程）进行了数值模拟。本文的其余部分安排如下。第 2 节介绍了模拟模型和放电配置。第 3 节讨论了基准工作，包括 PHD 与 COMSOL MultiphysicsTM 和 PPT 实验结果的比较。第 4 节详细介绍了在各种背景磁场下获得的主要模拟结果。第 5 节讨论了这些结果并提出了主要结论。

2.模拟模型

本研究建立了一个圆柱坐标下的双流体排放模型。我们首先讨论本研究中考虑的假设和近似值。物理量在极坐标方向上是均匀的，即所有变量的 ，矢量有三个分量。粒子流用漂移-扩散近似法描述。

电子能量分布函数（EEDF）是麦克斯韦式的。此外，中性原子和离子的温度被认为恒定在常温（300 K）。输入功率仅来自电子功率吸收和天线的焦耳热。只考虑了氩的基态电离反应，而忽略了其他化学反应，如潘宁电离、复合反应和发光。负离子和高价正离子也被忽略。接下来我们将考虑非麦克斯韦 EEDF，因为在实验中已经观察到 [34-37]。之前的实验表明，中性剖面以及电离导致的中性损耗会显著影响等离子体剖面或其动态行为 [38-42]。未来将考虑更复杂的中性动力学，如中性温度、压力平衡等。

值得注意的是，我们对平行磁场方向（平行方向）的粒子流采用了伏极扩散近似。在均匀背景磁场和周期性边界条件下，我们可以将平行方向输运模型从漂移扩散近似简化为伏极扩散近似。等离子体密度和温度相似，误差率约为 5%。我们的研究发现，在求解泊松方程时，低压和高磁场条件下会出现数值不稳定性。这是因为平行方向和垂直方向的扩散系数相差几个数量级。为了克服数值不稳定性，我们使用了伏极扩散近似。使用近似值的时间步长可比不使用近似值的时间步长大 100 多倍，即计算效率提高了 100 多倍。

等离子体模型

PHD 代码采用氩放电，模型中考虑的粒子是电子、离子和中性原子。等离子体模型的控制方程[43]如下。

n_α、D_α、Г_α (α=e,i,n)分别为颗粒密度、通量、扩散率和流动性。中性原子密度的计算公式为n_n = 2.41 x 10¹⁴ p_ncm^-3,其中p_n为气压，其单位为Pa。下标 e、i 和 n 分别表示电子、离子和中性原子。T_e 是电子温度，单位为 eV。此外，e、V 和ε₀ 分别是基本电荷、电势和自由空间的介电常数。与电荷分离 E_s 相对应的静电场由电势梯度给出，即 E_s = -▽V·H_i = 15.7eV 是第一电离能，基态电离的速率系数为。方程 (1) 中的通量 Г_α 是由修正的漂移-扩散近似法给出的：

其中，B₀ 为外部静磁场，μ_α 为迁移率。我们沿 z 方向设置了一个均匀的背景磁场，即B₀ = B₀e_z。电子迁移率μ_e 和电子扩散系数 D_e 的计算公式为和。其中，电子碰撞频率 ν_e 由 ν_e = ν_en+ν_ei [19] 得出。ν_en 和 ν_ei 分别是电子与中性原子的碰撞频率和电子与离子的碰撞频率。

其中，lnΛ 是库仑对数，δ_en 是电子与中性原子之间的碰撞截面。
方程 (5) 包括三元方程和三元未知数。我们可以用以下方程求解方程 (5)

采用平行通量的伏极扩散近似，方程 (10) 变为

其中是两极扩散系数。

中性原子的通量只考虑方程 (2) 中的扩散项

在公式 (3) 中，计算能量通量q_e时忽略了磁场对电子温度传输的影响。

p_tot 是电子吸收功率，计算公式如下：

等离子体电流 J_p 和感应电场 E_mf 由麦克斯韦方程计算得出。

麦克斯韦方程

高频电场 E_mf 和感应磁场 B_p 是通过麦克斯韦方程和等离子体电流演化方程计算得出的，计算公式如下：

其中，m_e、μ₀、ε_r、分别为电子质量、自由空间的磁导率和自由空间的相对介电常数。**J**_a 是天线电流，计算公式为

S_a是天线的横截面积，**I**_a 是电流振幅。总输入功率由等离子体吸收的电子功率和天线的焦耳热产生。电流幅值 I_a 根据电子吸收功率进行调整，以保持恒定的总功率。等离子体电流 J_p 由电子电流产生。因此，方程 (21)-(23) 也表示电子的运动。

初始条件

模拟的初始条件如下：

参考文献[43]和 PHD 模拟表明，稳态与一定范围的等离子体密度无关。然而，在低压和强磁场条件下，初始等离子体密度必须较大（约 10¹⁰cm^-3）。否则，等离子体密度不会增加。实验中采用的这种大等离子体密度值可归因于电容模式，而 PHD 并未考虑这种模式。另一个重要方面是电磁场初始条件的设定。我们计算了初始等离子体密度的 E_mf 和B_p的 “本征 “结构，并将其设置为初始条件。我们设定 E_mf=B_p=0，然后在放电过程开始前计算 20 个射频（RF）周期的麦克斯韦方程，即等离子体状态保持不变，包括n_e,n_i和T_e我们的评估结果表明，这种方法可在较短时间内收敛，而使用其他初始条件则会导致发散。

排放配置和边界条件

图2显示了放电配置。真空室的边界位于 r₀=0.075m处,铝屏蔽位于R₀=0.125m处。径向宽度为 10 mm、轴向宽度为 30 mm 的单环天线的位置分别为r_a=0.09 m和z_a=0.3m。这种放电结构与 PPT 的等离子体源区大致相同。

边界条件如下

我们只考虑了等离子体源区域的放电过程。为了减少边界的影响，我们对 z 方向上的所有变量都采用了周期性边界条件。

数值方法

对于流体方程 (1)-(3)，我们采用了空间交错网格，即标量占一个完整网格点，矢量占半个网格点。空间离散采用中心差分法。时间离散采用全隐式格式，并通过连续松弛（SOR）求解。

对于泊松方程 (4)，我们也采用了中心差空间离散法，并用 SOR 求解。为了克服时间步长的限制，我们采用了一种半隐式方法来求解电动势[45, 46]。这种方法的主要思想是利用当前时间步长的电场和等离子体密度分布来预测下一时间步长的电场。

方程 (15)-(23) 采用有限差分时域法。我们在空间离散化中使用了Yee 的交错网格[47]，在时间离散化中使用了两步交替隐式格式法。

PHD 与 COMSOL 和 PPT 的比较

本节分为两部分。在没有背景磁场的情况下，对 PHD 和 COMSOL Multiphysics 软件进行了简单的基准比较。此外，还将 PHD 代码的结果与 PPT 的实验结果进行了比较。除非另有说明，仿真参数和常数均取自表 2。

表 2.模拟参数和常数。

使用 COMSOL 进行基准测试

COMSOL 是一款商业软件，广泛用于许多科学模拟。考虑到其电感耦合等离子体模块的可靠性以及在背景磁场下获得的结果并不广为人知，我们使用该模块建立了一个基准。仿真参数如表 3 所示。为了使模拟条件更接近 COMSOL 的模拟条件，我们考虑了一系列与 COMSOL 一致的涉及可转移原子的化学反应[43, 48, 49]。轴向端板的边界条件与侧壁的边界条件相同，电子碰撞频率 ne 也与 COMSOL 中的相同。

表 3.模拟参数。

图 3 显示了电子密度随时间的演变。PHD 计算出的n_e 的时间演变与 COMSOL 得出的结果非常相似。电子密度稳态结构的形状也很相似；此外，如图 4（a）和（b）所示，当色标调整到相似范围时，它们的形状更加相似。最大等离子体密度的两个径向位置都在轴线上，轴向位置与天线位置相同。等离子体主要位于天线附近。

我们进一步比较了z=0.3 m处的电子密度径向分布，如图 4（c）所示。使用 COMSOL 得出的最大电子密度为 9.88x10¹¹ cm^-3,而使用 PHD 得出的最大电子密度为1.05x10¹²cm^-3。两者相差约 6%，这在放电过程中是可以容忍的。

PHD 模拟与 PPT 实验的比较

我们将本次 PHD 模拟的结果与之前 PPT 实验的结果进行了比较。需要注意的是，在 PPT 实验中，实验等离子体密度只能在扩散区进行评估，而 PHD 只模拟了源区，以节省计算资源。图 5 显示了最大离子密度与 (a) 功率 P₀,(b) 气体压力 p_n,和 (c) 磁场 B₀ 的关系。模拟值和实验值的大小相似，趋势一致。

压力曲线在 0.6 Pa 以下的一致性较好。原因可能是随着压力的增加，源区和扩散区之间的等离子体密度差增大。总的来说，模拟结果与实验结果在定性上是一致的，定量比较的结果并不十分相似。这可能与模拟配置的简化、建立模型时使用的若干假设、边界条件的选择等有关。

就磁场曲线而言，在强磁场条件下，PHD 与 PPT 之间的偏差似乎越来越大。可能的原因如下。(1) PPT 实验中颗粒的主要损失是侧壁和端板，而 PHD 模拟只考虑了侧壁。(2) 模拟中使用的是经典扩散，颗粒的传输可能比实验结果慢。因此，粒子的损失变小，等离子体密度变大。(3) 诊断位置离源头较远，因此会有一定的变化。此外，朗缪尔探针在强磁场下的诊断误差可能会变大。

总之，两次模拟的生长过程和电子密度曲线相似，数值差异很小。因此，PHD 代码在一定程度上是可靠的。

背景磁场变化时的放电过程

本节主要介绍无磁场和有磁场时的放电过程。放电情况在P₀=2000W、p_n=0.35Pa、0≤B₀≤1180G条件下进行。如图 6 所示，电子密度经历了线性增长和非线性饱和，最终达到稳态。可以看出，背景磁场越强，稳态等离子体密度越高。有趣的是，我们发现在磁场存在的情况下，密度在达到稳态之前迅速增加。我们称之为第二次密度跃迁，这与螺旋子等离子体中的模式转换类似。

在B₀ = 1180G的情况下，第二密度跃迁时间约为 t=3000 T_rf。B = 1180G时的电子密度约为无磁场时的 82 倍。这里有两种可能的解释。首先，可能存在一个共振等离子体密度，它是磁场 B₀ 和气体压力 p_n 的函数。在这个特定的等离子体密度下，电磁波和等离子体达到共振状态。这种共振会导致高电子功率吸收。其次，磁场抑制了等离子体的传输，从而减少了等离子体在壁上的损耗。因此，等离子体密度会增加。然而，第二次密度跃迁的确切起源仍不清楚。

图 7 显示了 n_e 的稳态结构。随着磁场的增大，轴向分布的均匀性增加，等离子体逐渐从局部区域扩散到整个真空区域。

图 8 显示了稳定状态下电子温度 Te 和电子吸收功率 Ptot 的周期平均分布。Ptot 的负值意味着电子失去能量，电子的能量转化为波能。在此，我们将简要讨论放电发生的方式。电子功率的吸收会提高电子温度，从而使电离截面变大，进而产生大量电离过程。如果源项与等离子体损耗相平衡，就会达到稳定状态。随着磁场的增加，电子温度降低。除B = 50G的情况外，最大 Te 的轴向位置与天线的轴向位置相同。此外，在有磁场的情况下，最大电子功率吸收发生在中心轴附近，而在无磁场的情况下，则出现在边界附近，如图 8（b）所示。在螺旋推进器配置中也观察到了天线附近或源壁附近的高温区[50-54]。PHD 模拟结果在一定程度上再现了这一现象。另一个有趣的现象是，当引入磁场时，电子功率吸收区域在 z 方向从一个峰值变为两个峰值。通过仔细观察，我们发现径向图案的数量会随着磁场的增加而增加。具体解释如下

通常认为，由于阻尼率较弱，螺旋波可以向中心轴传播。周期性边界条件很可能产生驻波解。在实验中也观察到了驻波螺旋波[55-57]。因此，我们推测这些系统中的主波是驻波螺旋子波。图 9 显示了不同时刻极性电场E_mfθ和感应磁场B_pz的振幅。在初始时刻（t = 40 T_rf），E_mfθ的振幅从天线位置开始衰减。当触发第二次跃迁（t = 3160 T_rf）和第三次跃迁（t = 3400 T_rf）时，它变成了驻波结构，等离子体中的振幅高于天线位置。然而，电场在稳态（t = 20000 T_rf）时并不是一个明显的驻波结构，这可能是在实验中不易测得高驻波比（用等离子体中最大E_mfθ/最小E_mfθ表示）结果的原因。驻波产生了电子功率吸收的两个轴向和中心峰值。磁场越大，密度越高。密度越高，径向波长越短。因此，径向图案的数量也会增加。

从图 6-8 可以看出，电子密度、电子温度和电子功率吸收之间存在合理的关系。模拟结果符合放电定律。

图 9.在 2000 W、0.35 Pa、1180 G 的情况下，不同时刻 (a) 极电场E_mfθ的振幅（单位 V m^-1）和 (b) 感应磁场 B 的振幅（单位 T），其中t = 40 T_rf、t = 3160 T_rf、t = 3400 T_rf、t = 20000 T_rf分别对应图 6 中的 A、B、C、D 点。

讨论和结论

本研究建立了一个包含电子功率吸收、电离过程和传输过程的双流体放电模型，并利用该模型获得了数值解。PHD 代码可用于模拟低压（0.5 Pa）和强背景磁场（1kG）下的放电过程。在无磁场条件下得到的结果与使用 COMSOL 得到的结果基本一致，在典型的螺旋子放电参数下得到的结果基本符合放电过程的规律。扫描参数的结果在定性上与实验结果一致，但在定量上并不一致。与无背景磁场的情况相比，强磁场下的等离子体密度要高出一个数量级，这与实验验证的放电定律是一致的。在有磁场的情况下，稳定状态下的最大电子功率吸收发生在中心轴附近。这可能是由于驻留的螺旋波造成的。在完整的放电过程中发现了第二次密度跃迁。这些结果有助于我们理解放电机制。特别是，随着功率的增加出现稳态密度跃迁是一个重要的结果。然而，由于应用了周期性边界条件，等离子体损耗不足，因此无法获得全面的结果。这是本研究的局限性。总的来说，本研究的模拟是初步的，在边界条件、放电配置、垂直温度传输、对称性破坏、天线优化、陨落原子反应等方面还需要进一步研究。