PIC仿真增大真空介电常数对等离子体周向波的影响

The Influence of Increasing Vacuum Dielectric Constant on Plasma Azimuthal Wave in PIC Simulation

摘要: 在霍尔推力器全粒子PIC (Particle in cell)模拟方法中，常常采用增大真空介电常数的方法来增大等离子体德拜长度，进而在划分网格时可以将网格尺寸放大，减小网格数量从而达到减小计算量的目的。但由等离子体形成的周向振荡对真空介电常数变化反应比较敏感，相应等离子体参数随真空介电常数的变化是否符合离子声波色散规律需要得到验证。本文以2.5D-3V的霍尔推力器周向–径向模型为基础，通过改变真空介电常数放大倍数对等离子体分布情况、空间电势分布情况及模拟所得离子声波参数与理论值对比发现：在研究霍尔推力器等离子体周向振荡相关问题时，真空介电常数扩大倍数越大，所得周向波参数与理论估算值偏差越大，原则上不宜采用真空介电常数扩大的方法进行模拟，若所研究区域范围较大，可适当扩大真空介电常数至4倍以内。

Abstract: In the PIC (particle in cell) of Hall thruster, increasing the dielectric constant of vacuum is often used to increase the Debye length of plasma, and then the mesh size can be enlarged and the num-ber of meshes can be reduced to reduce the calculation amount. However, the azimuthal oscillation formed by the plasma is sensitive to the change of the vacuum dielectric constant. Whether the change of the corresponding plasma parameters with the vacuum dielectric constant conforms to the law of ion acoustic dispersion needs to be verified. This paper is based on the circumferential and radial model of hall thruster of 2.5D-3V. By changing the amplification factor of the vacuum di-electric constant, the plasma distribution, the space potential distribution and the simulated ion acoustic wave parameters were compared with the theoretical values, and it was found that: the larger the amplification of vacuum dielectric constant is, the larger the azimuthal wave parameters are. In principle, it is not suitable to use the method of expanding the vacuum dielectric constant to simulate the larger deviation from the theoretical estimate. If the research area is large, it can be appropriately expanded to within 4 times.

1. 引言

霍尔推力器工作时，当稳定的霍尔漂移运动遭到破坏时，其导向中心发生偏移，在电场作用下，电子会产生向阳极的迁移运动，即电子传导 [1]。相关研究表明 [2] [3] [4]：在霍尔推力器放电过程中，等离子体振荡会引起电子异常传导，进而产生电子传导电流。近年来，有很多针对等离子体振荡对电子传导电流的影响的相关研究 [5]，都是通过霍尔推力器放电通道全粒子PIC模型进行的。在粒子模拟过程中，为减小计算量，加快计算速度，常常采用Szabo [6] 提出的增大真空介电常数的方法。但这种方法会导致模型中的物理过程发生变化，尤其是会削弱等离子体之间的自洽电场强度。而等离子体的振荡频率与真空介电常数的大小息息相关。因此，必须考虑此方法对所得模拟结果的影响。近年来，有研究考虑了PIC模拟中真空介电常数的改变对电子运动的影响 [7]。本文通过建立以霍尔推力器放电通道周向–径向(θ-r)为研究平面的2.5D-3V PIC数值计算模型，研究放大真空介电常数不同倍数，对空间等离子体分布情况、电势分布情况及周向振荡情况所产生的影响，来判断能否通过真空介电常数增大来进行模型加速。

2. PIC模型

本文建立以霍尔推力器放电通道周向–径向(θ-r)为研究平面的2.5D-3V PIC数值计算模型。模拟霍尔推力器放电通道加速区为一环形平面区域，如图1所示。

将此环形区域沿着半径R方向切开。由于电子漂移不稳定性波长在mm量级 [8] [9] [10]，而数值研究 [9] 表明，曲率效应对于霍尔推力器(通道半径通常为几十厘米的量级)的影响并不大。因此对模型进行简化处理，这里忽略曲率对模型的影响，认为R方向为无限大。在这种情况下仿真区域近似为一个矩形。对应于柱面坐标系的 $\left(\theta ,r,z\right)$ 建立直角坐标系 $\left(x,y,z\right)$，如图2所示。

Figure 1. Simulate the region diagram

图1. 模型计算区域图

Figure 2. PIC model

图2. PIC模型

霍尔推力器放电通道中，离子或电子数密度在10¹⁷~10¹⁸在模型中，将磁场沿着与推力器通道环形壁面垂直(沿Oy负方向)设置，初始磁场强度设为 $B_0=200\text{G}$，电场沿推力器轴向(沿Oz正方向)设置为固定电场，初始电场强度为 $E_0=1\times {10}^4\text{V}/\text{m}$。PIC计算平台对于边界条件的处理主要包括确定边界电势以及对仿真粒子与边界相互作用的设置。在此模型中，主要涉及了三个方向的边界条件：x方向(周向)边界，采取周期性边界条件进行处理；y方向(径向)边界，为霍尔推力器双壁面结构；如下图3所示：

Figure 3. Model boundary

图3. 模型边界情况

同时，为保证模型稳定，对z方向边界做如下处理：当离子碰撞到z方向上边界 $z=z_{\max }$，则离子消失，在z方向下边界 $z=z_0$ 处以离子初始能量 $T_i=0.1\text{eV}$ 给入新离子。当电子碰撞到z方向下边界 $z=z_0$，则电子消失，在z方向上边界 $z=z_{\max }$ 处以电子初始能量 $T_e=5\text{eV}$ 给入新电子。以这种形式进入模拟区域的带电粒子，其位置在x、z方向上随机分布，而y方向为了保证鞘层形成过程不受影响，粒子位置按其打出模拟区域的y方向位置重新分布。对于粒子碰撞z方向边界的其他几种情况，认为粒子以初始能量大小( $T_e=5\text{eV}$， $T_i=0.1\text{eV}$ )反打回模拟区域。模型的各项参数设置如下表1：

Table 1. PIC model parameter

表1. PIC模型模拟参数

模型验证

本模型模拟中电场是通过等离子体运动不断自洽求解的，模型计算中存在电场扰动，但磁场是沿径向固定施加的，并不存在磁场扰动，所以模拟中所产生的波均为静电波。根据文献 [5] 中相关的模拟结果，本文中所模拟的周向波从参数上来说符合离子声波的色散规律，接下来对模拟区域进行试算说明，并对离子声波进行验证。

采用表1中参数以及在此模型其他参数不变的情况下，增大初始离子温度为20 eV，使离子温度值高于电子温度进行对比模拟。图4分别是两种工况下电子密度在模拟区域的分布情况。在初始离子温度较低的情况下，放电区域内有明显的周向波形成现象，但当离子温度升高至高于电子温度时，由于离子朗道阻尼的出现，对离子声波的形成产生了抑制，这里可以判定，模拟中所产生的周向波均为离子声波。

因此，接下来的研究中，本文将根据离子声波色散关系对模拟结果进行讨论和分析，根据动理学理论 [11]，离子声波波长：

$\lambda =2\pi {\lambda }_D\sqrt{2}=2\sqrt{2}\pi \sqrt{\frac{{\epsilon }_0{k}_b{T}_e}{e^2{n}_e}}$

频率：

$\omega =k\cdot v_{i,0}\pm \frac{k{c}_s}{\sqrt{1+k^2{\lambda }_D^2}}=\frac{k{c}_s}{2\pi \sqrt{1+k^2{\lambda }_D^2}}$

其中， $c_s$ 为离子声速；k为离子声波沿波矢方向波数。

相速度：

$v_{ph}=\omega \lambda =c_s\sqrt{\frac{2}{3}}$

(a) (b)

Figure 4. The distribution of electron density in the simulated discharge region under different ion temperature conditions. (a) Electron density distribution at ion temperature of 0.1 eV; (b) Electron density distribution at ion temperature of 20 eV

图4. 不同离子温度放电区域内电子密度分布对比。(a) 离子温度为0.1 eV电子密度分布；(b) 离子温度为20 eV电子密度分布

3. 模拟结果分析

3.1. 增大真空介电常数对等离子体分布的影响

本文分别以：原真空介电常数以及将真空介电常数扩大4倍、9倍、49倍、100倍和400倍的情况进行模拟评估。

图5通过调整不同真空介电常数得到模拟区域电子密度的分布云图。从四种结果来看，在真实的真空介电常数 ${\epsilon }_0=8.85\times {10}^{-12}\text{F}/\text{m}$ 情况下，模拟放电区域中心点电子密度的最高值约为3.6 × 10¹⁷/m³，在四种工况中最高；当真空介电常数在此基础上扩大9倍后，模拟区域中心电子数密度的最高值约为3.0 × 10¹⁷/m³，相比原工况有所下降，当真空介电常数在原基础上扩大100倍后，模拟放电区域中心点电子数密度的最高值约为2.6 × 10¹⁷/m³，此时模拟放电区域中心的电子数密度最大值与原工况相比降低了1 × 10¹⁷/m³。以上工况说明：在真空介电常数合理增大的前提下，真空介电常数增大的越大，所对应模拟区域中心电子密度越小。介电常数是综合反映介质内电极化行为的物理量，介电常数越大，介质产生的感应电荷削弱原外加电场的能力越大，即原外加电场减小的越多，介质的极化特性越强。这里将真空介电常数增大，会使得模拟区域中粒子电荷更趋于平均化，造成了模拟区域中心电子密度下降的现象。但在真空介电常数扩大400倍时，真空介电常数扩大过大，已经影响了波的形成，介质内部粒子所形成的反向电场大大影响了原本施加的外电场，与真实真空介电常数情况下相比，真空介电常数扩大400倍后对电荷的平均化造成了波形被破坏，已经无法准确模拟推力器内相关物理现象。

(a) (b) (c) (d)

Figure 5. The distribution of electron density in the simulated discharge region under different ion temperature conditions. (a) Electron density distribution ( $\epsilon =4{\epsilon }_0$ ); (b) Electron density distribution ( $\epsilon =9{\epsilon }_0$ ); (c) Electron density distribution ( $\epsilon =100{\epsilon }_0$ ); (d) Electron density distribution ( $\epsilon =400{\epsilon }_0$ )

图5. 不同离子温度放电区域内电子密度分布对比。(a) 电子数密度分布( $\epsilon =4{\epsilon }_0$ )；(b) 电子数密度分布( $\epsilon =9{\epsilon }_0$ )；(c) 电子数密度分布( $\epsilon =100{\epsilon }_0$ )；(d) 电子数密度分布( $\epsilon =400{\epsilon }_0$ )

3.2. 增大真空介电常数对空间电势分布的影响

图6是不同真空介电常数，模拟放电区域空间电势在径向及周向的分布变化情况。随着真空介电常数增加，模拟放电区域内电势下降，且根据图6(a)中两侧壁面附近电势变化的斜率上看，较大的真空介电常数会造成鞘层明显增厚。图6(b)中从周向电势变化波形上看，其趋势与电子密度分布中大致相似，即在5 mm的模拟周向长度内形成约3个周向波，波长并不随着真空介电常数发生变化。

图7是选取原真空介电常数为 ${\epsilon }_0=8.854\times {10}^{-12}\text{F}/\text{m}$，及对其分别扩大9倍、49倍、100倍、400倍情况下，模拟区域内最高点电势值和壁面电势的变化情况。从图中可以发现，随着真空介电常数扩大倍数的提高，模拟区域电势降下降，而壁面负电势增大，说明增大真空介电常数，导体壁面负电荷聚集，鞘层电势降更大。在原真空介电常数 ${\epsilon }_0=8.854\times {10}^{-12}\text{F}/\text{m}$ 条件下，导体壁面附近电势为−0.42 V。

(a) (b)

Figure 6. Comparison of radial and azimuthal changes of potential at different vacuum permittivity. (a) Comparison of radial changes of potential at different vacuum dielectric constants; (b) Comparison of azimuthal changes of potential at different vacuum dielectric constants

图6. 不同真空介电常数情况下电势沿径向及周向变化对比。(a) 不同真空介电常数电势沿径向变化对比；(b) 不同真空介电常数电势沿周向变化对比

Figure 7. The variation of center point potential and wall potential in the simulated dischar-

ge area at different vacuum dielectric constants

图7. 不同真空介电常数模拟放电区域中心点电势及壁面电势变化情况

对不同真空介电常数情况下模拟得到的传导电流进行统计得到图8。传导电流随真空介电常数增大而下降。由经典传导估算该模型参数下传导电流为：

$j_{ez}\approx \left(\frac{{\upsilon }_{ea}}{{\omega }_{ce}^2}\right)\frac{e^2{n}_e}{m_e}{E}_z=m_e{n}_e{\upsilon }_{ea}\frac{E_z}{B_r^2}\approx 0.056\text{A}$

而模拟中，在原真空介电常数情况下所得到的传导电流约为0.1 A。这是由于除了经典传导之外还有其他电子传导形式存在。但随着真空介电常数增加，模拟所得到的传导电流大幅度下降，已经远远偏离计算的经典传导电流和原真空介电常数条件下模拟的到的传导电流。

Figure 8. The electron conduction current varies with the vacuum dielectric constants

图8. 电子传导电流随真空介电常数变化

3.3. 周向振荡参数分析

为了更准确地讨论增大真空介电常数对周向波频率产生的影响，图9截取了原真空介电常数在模拟区域处及 $x=L_x/2$ 处，周向及径向电势随时间变化情况，其中图9(a)中显现出的波的斜率为周向波的相速度，分别统计不同真空介电常数情况下，波的相速度变化。具体方法是：分别取模型稳定后，与下图等时间间隔的不同时段各工况的周向电势随时间变化的图像三组，计算图中相速度，每种离子质量取三组相速度的平均值统计结果，并与不同真空介电常数情况下计算的理论相速度进行对比，结果如表2。表2中同时对模拟所得周向波的波长、频率以及周向波形成时间与估算理论值进行了对比。

(a) (b)

Figure 9. The azimuthal and radial potential varies with time. (a) The azimuthal potential varies with time; (b) The radial potential varies with time

图9. 周向及径向电势随时间变化情况。(a) 周向电势随时间变化情况；(b) 径向电势随时间变化情况

Table 2. Theoretical and simulation value of the azimuthal wave parameters at different vacuum permittivity

表2. 不同真空介电常数周向波各项参数与理论值对比

4. 总结

从周向波形成的角度：适度扩大真空介电常数，不影响霍尔漂移运动形成，在真空介电常数扩大100倍以内，模拟区域中均能形成稳定周向波，但较大的真空介电常数会影响电子漂移运动，当将真实的真空介电常数扩大400倍时，周向波已经无法形成。而从模拟的周向波各项参数来看：无论真空介电常数扩大多少倍，均会使周向波的各项参数偏离理论值，这说明此时的模拟结果已经不再符合离子声波的色散规律，真空介电常数改变对电子的漂移运动产生了不明确的新的运动机制，故，通过增大离子质量来减小此PIC模拟程序计算量的方法不再可靠。因此本文认为：为保证此模型计算准确性不宜采取扩大真空介电常数方法进行处理，若模拟区域较大，为满足计算需要可以适当扩大真空介电常数4倍以内。