基于COMSOL辉光放电带电粒子特性模拟分析

Simulation and Analysis of the Characteristics of Charged Particles in Glow Discharge Based on COMSOL

摘要: 为了研究空芯光纤中氩气的放电特性，采用软件COMSOL的有限元方法计算相应的数学等式，基于等离子体放电理论，建立了圆柱型二维直流辉光放电流体模型。得到氩气放电在不同电压下反应腔内电子流密度、电子密度、激发态氩原子密度、电子温度等粒子密度的变化和在阴极区不同电压下粒子密度和电子温度的变化，最后对距离阴极区0.1 m的空间电荷密度做了分析。结果表明：随着电压的增大，阴极区粒子密度呈现出逐渐增加的趋势并且电子流密、电子密度和电子温度在100 V时变化明显。由于反应腔内的氩气量的影响，激发态氩原子密度和空间电荷密度随着电压的增大趋于稳定。

Abstract: In order to study the discharge characteristics of argon in hollow-core optical fiber, the finite ele-ment method of the software COMSOL was used to calculate the corresponding mathematical equa-tions. Based on the plasma discharge theory, a cylindrical two-dimensional DC glow discharge fluid model was established. The changes of particle density such as electron current density, electron density, excited state argon atom density, and electron temperature in the reaction chamber of ar-gon discharge under different voltages, and the changes of particle density and electron tempera-ture in the cathode area under different voltages were obtained. Finally, the space charge density at a distance of 0.1 m from the cathode area was analyzed. The results show that the particle density in the cathode region tends to increase gradually with the increase of the voltage, and the electron current density, electron density, and electron temperature change obviously at 100 V. Due to the influence of the amount of argon in the reaction chamber, the excited state argon atom density and space charge density tend to be stable with the increase of voltage.

1. 引言

空芯光子晶体光纤(photonic crystal fiber, HC-PCF)由Cregan等人基于光子带隙原理研制，实现了光在空芯中传播 [1] 。HC-PCF也利用其特性 [2] [3] ，扩展到很多传感器领域研究 [4] [5] [6] 。气体放电是我们经常看到的现象，如天空中的闪电和家庭中的日光灯，这是由于在电场作用下或其他激励下气体产生的电离。而惰性气体放电成为气体放电研究中的热点，国内的石鑫 [7] 和国外的Hong Y J，Lee S M，Kim G C等 [8] 都做出对应的研究。

辉光放电会产生大量的等离子体，它是物质的第4种状态，并且包含大量的带电粒子。等离子体的特性和辉光放电的特性紧密相连，如放电的电压、温度、气压等。国内的邵先军第一次使用有限元法建立了一维氩气直流辉光放电自洽流体模型，并对带电粒子的浓度、二次电子和总通量密度做出了研究 [9] 。高媛等以氢气为研究对象，对放电过程中的21种碰撞反应类型进行了模拟，起作用的有11种非弹性反应类型，并且这些是形成等离子体的主要反应 [10] 。2021年国外的Lisin Evgeny A.和Kononov Evgeny A.等首次对粒子悬浮区域的粒子间相互作用和地层电场进行了定量研究。在辉光放电中也首次观察到离子流对下游粒子的放电效应 [11] 。

1962年Ward等人在基于带电粒子的连续性方程和泊松方程的基础上对直流辉光放电等离子体的特性进行了研究，并得到了丰硕的成果 [12] 。近年来学者对等离子体粒子 [13] [14] 和电子温度 [15] [16] 做了大量的研究。由于在放电的过程很复杂，只是简单的实验不能很好的理解放电过程，数值模拟以其制作成本较低等优势并可以很好的解释放电原理，成为了现代研究辉光放电方式的热点 [17] [18] 。为了更加直观的观察到粒子的变化、进一步了解粒子的反应，本文在以上研究的基础上，笔者通过COMSOL有限元软件，建立二维的直流辉光放电模型，充入氩气惰性气体进行模拟仿真。研究不同电压下的电子流密度、离子流密度、空间电荷密度、电子温度等的变化情况，对氩气直流放电特性进行了分析。

2. 模型建立

2.1. 直流氩气放电几何模型

2.1.1. 直流放电实验装置

Figure 1. DC discharge experimental device

图1. 直流放电实验装置图

为了对实验整体有一个更好的理解并加深对实验过程的印象，现对直流放电装置进行构造。如图1所示，模型中左边的为阴极右边为阳极，两电极间的距离为0.568 m，气压稳定在3 Torr，外接一个直流电源和电压表，通过观察电压表能够比较直观的观察到电压的变化。为了研究充入氩气后的反应，实验开始前需要把反应腔里的气体抽出，以免对实验造成干扰。

2.1.2. 二维流体网格图

Figure 2. 2D fluid mesh diagram

图2. 二维流体网格图

图2为基于COMSOL建立的二维流体网格图，对应于图1所示，左边为阴极，右边为阳极。通过外电压，反应腔内由气体转化为导电的粒子，如电子、离子、激发态氩离子等。

2.2. 直流氩气放电流体模型

2.2.1. 数学模型

汤生放电是辉光放电的基础，它们之间最主要的差别是放电电流密度，后者比较大，空间电荷效应也起重要的作用。两者之间的关系可以表示为：

${\nabla }^{2}V=-\nabla \cdot E=\frac{e}{{\epsilon }_0}\left(n_e-n_i\right)$ (1)

其中，V是电势，E是场强， ${\epsilon }_0$ 为真空介电常数。

$\frac{\partial n_e}{\partial }+\nabla \cdot {\Gamma }_e=R_e-\left(u\cdot \nabla \right)n_e$ (2)

$\frac{\partial n_{\epsilon }}{\partial }+\nabla \cdot {\Gamma }_{\epsilon }\text{+}E\cdot {\Gamma }_e=R_{\epsilon }-\left(u\cdot \nabla \right)n_{\epsilon }$ (3)

其中 $n_e$ 是电子密度， ${\Gamma }_e\text{=}-n_e\left({\stackrel{\to }{\mu }}_e\cdot \stackrel{\to }{E}-{\stackrel{\to }{D}}_e\cdot \nabla n_e\right)$ 计算的是电子流密度，表示在E和浓度差作用下的电子流密度； ${\stackrel{\to }{\mu }}_e$ 代表电子的迁移率； ${\stackrel{\to }{D}}_e$ 代表电子的扩散率； $R_e$ 代表电子产生或消失的速度，可以由化学反应系数来表示； ${\Gamma }_{\epsilon }\text{=}-n_{\epsilon }\left({\stackrel{\to }{\mu }}_{\epsilon }\cdot \stackrel{\to }{E}-{\stackrel{\to }{D}}_{\epsilon }\cdot \nabla n_{\epsilon }\right)$ 计算的是电子能流密度， ${\stackrel{\to }{\mu }}_{\epsilon }\text{=}\left(5/3\right){\stackrel{\to }{\mu }}_e$ ；

考虑中性粒子的扩散作用：

$\frac{\partial n_n}{\partial t}-D_n{\nabla }^2{n}_n=S_n$ (4)

其中 $D_n$ 为中性粒子的扩散系数； $S_n$ 为源项即：

$S_n={\displaystyle {\sum }_G{k}_G}{n}_1{n}_2-{\displaystyle {\sum }_L{k}_L{n}_1}{n}_2$ (5)

其中 $n_n$ 为中性粒子密度，G为相应粒子的产生项，L为相应的粒子损失项。

电子源项为：

$R_e={\displaystyle {\sum }_{j=1}^N{n}_n}{k}_j{n}_e$ (6)

其中j表示反应腔中第j种反应； $k_j$ 为第j种反应的反应系数。

将弹性碰撞损失的能量和非弹性碰撞损失的能量加在一起可得电子能生成率为：

$R{}_{\epsilon }={\displaystyle {\sum }_{j=1}^N{n}_n{k}_j{n}_e\Delta {\epsilon }_j}$ (7)

其中 $\Delta {\epsilon }_j$ 为第j种反应能量的损耗。

$k_j={\displaystyle {\int }_0^{\propto }\epsilon \sigma \left(\epsilon \right)f\left(\epsilon \right)\text{d}\epsilon }$ (8)

其中 $\sigma \left(\epsilon \right)$ 为电子能的碰撞截面函数， $f\left(\epsilon \right)$ 为电子能量分布函数。

重物质输运(非电子物质)主要是求解混合物平均扩散模型：

$\rho \frac{\partial w_k}{{\partial }_t}+\rho \left(\stackrel{\to }{u}\cdot \nabla \right)w_k=\nabla \cdot {\stackrel{\to }{j}}_k+R_k$ (9)

$\rho =q\left({\displaystyle {\sum }_{k=1}^N{Z}_k{n}_k}-n_e\right)$ (10)

其中 $w_k$ 是第k种物质的质量分数； $j_k$ 是扩散通量向量；并且 $j_k=\rho w_k\cdot V_k$ ， $V_k$ 是第k种的扩散通量矢量； $R_k$ 是第k种物质的速率表达式，由COMSOL进行计算； $\rho$ 是空间电荷密度,目前国内外使用的都是混合模型。

2.2.2. 边界条件

辉光放电是自持放电，它由离子轰击阴极产生的二次电子发射来维持。电子的二次发射表示为：

$\hat{n}\cdot {\Gamma }_{e,s}={\displaystyle {\sum }_{k=1}^N{\gamma }_k{n}_n\left({\Gamma }_k\cdot \hat{n}\right)}$ (11)

其中 $\hat{n}$ 是垂直于阳极或阴极表面的矢量； ${\gamma }_k$ 是二次发射系数； ${\Gamma }_k$ 是表面的离子流。

$\frac{\partial {\rho }_s}{{\partial }_t}=\hat{n}\cdot {\stackrel{\to }{j}}_i+\hat{n}\cdot {\stackrel{\to }{j}}_e$ (12)

其中 $\hat{n}$ 是反应腔体表面的输出单位矢量； ${\stackrel{\to }{j}}_e$ 为反应腔表面的总电子流； ${\stackrel{\to }{j}}_i$ 为反应腔表面的总离子流。

${\sigma }_s=-\hat{n}\cdot \stackrel{\to }{D}$ (13)

其中方程(13)表示表面电荷密度； ${\sigma }_s$ 用作边界条件。

2.2.3. 氩气碰撞化学反应

表1是氩气的等离子体化学反应，比如弹性、激发和电离反应的不同横截面能量。

Table 1. Argon collosion reaction

表1. 氩气碰撞反应

表2是氩气表面反应，表中规定了粘附系数，描述了亚稳态原子与壁接触时恢复为基态原子的概率。

Table 2. Surface reaction

表2. 表面反应

3. 仿真结果与分析

使用COMSOL Multiphysics Plasma模块对氩气进行了仿真，计算得到了电子流密度、电流密度、激发态氩原子密度和电子温度的分布。通过改变电压在25 V~300 V，研究分析粒子在阴极区随电压的分布。最后研究了距离阴极0.1 m处电压对空间电荷密度分布的影响。

3.1. 粒子密度和电子温度分布

在两极设置电压为25 V、125 V、225 V，气压稳定在3 Torr。计算模拟粒子密度和电子温度的变化，其中横坐标左端为反应腔的阴极，横坐标是距离阴极的长度，模拟结果如下。

Figure 3. Spatial distribution of the electron current density

图3. 电子流密度空间分布

图3是不同电压下电子流密度的分布图。可以看出，电子流密度随着电压的增加而增加。当电压为25 V和125 V时，电子流密度变化不大。距离阴极0.1 m处，电压为225 V时，电子流达到0.003 A/m²的峰值，然后减小并趋于稳定。由于电压的增加使得初始电子和二次电子产生越来越高的能量轰击阴极，引起电子和激发气体粒子的雪崩。能量越高，与氩原子碰撞的概率就越高，也造成了大量的电离。没有发生弹性碰撞和非弹性碰撞的粒子在积累能量后，使正柱区的气体电离，在电场中大量电子被加速，产生流向阳极的电子流。

Figure 4. Spatial distribution of the electron density

图4. 电子密度空间分布

图4是不同电压下电子密度的分布图。电子密度随电压的增加而增加，在距阴极区0.1 m的距离处达到最大值，然后减小并趋于稳定。当电压为125 V和225 V时，两者的电子密度在10¹⁶ 1/m³附近变化大。当电压为25 V时，电子密度整体在10¹⁰ 1/m³附近变化，与电压为125 V和225 V时的电子密度有很大不同。当电压为125 V和225 V时，由于电压的增加达到自持放电的条件，电子数继续增加，二次电子数也相对增加，阴极区域的电子密度也不断增加。

Figure 5. Spatial distribution of the excited argon atom density

图5. 激发态氩原子密度空间分布

图5是是不同电压下激发态氩原子密度的分布图。当电压为125 V和150 V时，激发的氩原子密度在阴极达到最高值，然后有下降趋势。由于达到自持放电电压，二次电子的数量也增加，电子集中在阴极，之后电子与氩原子碰撞产生激发态氩原子，从而增加激发的氩原子的数量。当电压为50 V时，由于未达到自持放电的电压，阴极附近的电子为慢电子，容易与氩原子碰撞形成激发态氩原子。此外，与125 V和225 V时的激发态氩原子密度相比，25 V时的激发态氩原子密度低为10⁴ 1/m³，这是由于低电压放电管中的电子数少，电子速率低，氩原子碰撞几率低，导致二次电子数量少。

Figure 6. Spatial distribution of electron temperature

图6. 电子温度空间分布

图6是不同电压下电子温度的分布图。电压为125 V和225 V时阴极区温度较高，温度峰值随电压升高而增大，然后呈现下降趋势。分析得出的结论是，由于阴极区域的电场相对较强，电子获得较高的能量，并且在阴极表面发生加热，其中额外的电子以高漂移速度从阴极表面发射使电子温度在阴极区域达到峰值。当电压为25 V时，整体电子温度保持在1.6 eV左右，阴极区域呈现相对稳定的趋势。负辉区的电子温度变化非常缓慢，这表明电子温度受电场的影响很大。

3.2. 不同电压下阴极区粒子密度和电子温度的分析

为了进一步研究电压对粒子密度和电子温度的影响，在电极阴极区进行分析。

Figure 7. Spatial distribution of electron current density at different voltages

图7. 不同电压下电子流密度空间分布

图7是电子流密度的分布图。从图中可以看出，在25 V~100 V期间，电子流的密度在0 A/m²之间波动。在100 V~250 V时，电子流密度突然增加到1.4 A/m²，这是由于达到维持二次电子连续放电所需的电压引起的。当电压为275 V时，电子流密度趋于波动。这是由于电位降主要发生在阴极区域附近，导致电子能量损失，导致电子流密度降低，但总体趋势仍在增加。阴极的电子流密度主要由外部电压提供。

Figure 8. Spatial distribution of electron density at different voltages

图8. 不同电压下电子密度空间分布

图8是电子密度的分布图。电压为25 V~100 V时电子密度保持在10⁹ 1/m³，电压超过100 V时电子密度增至10¹² 1/m³，这主要是由于二次电子激发。当电压为225 V时，随着电压的增加，电子密度趋于稳定。

Figure 9. Spatial distribution of excited argon atomic density at different voltages

图9. 不同电压下激发态氩原子密度空间分布

图9是激发态氩原子密度的分布图。电压在25 V~100 V之间的整体变化稳定，并有轻微下降趋势，之后激发态氩原子密度开始急剧增加。虽然225 V后电压一直在增加，但整体变化是平缓的。由于反应腔中氩气的量是恒定的，所以随着电压的不断增大，激发态氩原子的量越来越小。

Figure 10. Spatial distribution of the electron temperature at different voltages

图10. 不同电压下电子温度空间分布

图10是电子温度的分布图。随着电压的不断增加，电子温度也不断提高。当电压为25 V~100 V时，电子温度保持在2 eV~2.2 eV之间。在阴极区域，100 V后电子温度明显升高，这是由于随着电压的增加，电场也在增加，阴极中的高能离子由于场强的增加轰击阴极，这使得电子温度不断升高。

3.3. 不同电压下阴极区粒子密度和电子温度的分析

为了进一步探究电压对粒子密度的影响，模拟测试了空间电荷密度下的分布情况。其中图11是25 V、125 V、225 V电压下的对比分析，图12是在距离阴极区0.1 m处，电压从25 V~300 V时空间电荷密度随电压变化的情况。具体研究如下。

Figure 11. Space charge density spatial distribution at different voltages

图11. 不同电压下空间电荷密度空间分布

图11是电压为25 V、125 V和225 V时空间电荷密度的分布图。阴极区的空间电荷密度随电压的增加而增大，阴极区的空间电荷密度较高，然后逐渐减小。在距阴极0.1 m的距离处，空间电荷密度缓慢但最终达到零。由于电子密度和离子密度在正柱区几乎相互抵消，而电子密度对外中性，所以这个区域的空间电荷密度几乎为零。在距阴极0.1 m的距离处，空间电荷密度在125 V和225 V时变化相似，稳定在10⁻¹⁰ C/m³，而在25 V时空间电荷密度低两个数量级。

Figure 12. Spatial charge density of 0.1 m from the cathode varies with voltage

图12. 距离阴极0.1 m空间电荷密度随电压变化

为了进一步观察电压对距阴极某一点的空间电荷的影响，如图12选择距阴极0.1 m的距离来测量和分析空间电荷密度的变化。空间电荷密度随电压的增加逐渐减小，当电压为150 V时趋于稳定，主要是由于电子持续增加和离子抵消。

4. 结论

本文基于COMSOL Multiphysics，采用二维流体模型模拟氩气直流放电，虽然相比于实际放电进行了一定简化，但也能客观反映不同粒子间的变化。本文对电子流密度、电子密度、激发态氩原子密度、电子温度、空间电荷密度等进行仿真分析，得出以下结论：

1) 电压在25 V~100 V区间时，阴极区的电子流密度在0 A/m²之间波动，当电压超过100 V时放电模式发生变化，产生大量的二次电子，导致电子流密度和电子密度以及阴极中的电子急剧增加，电流密度主要由施加的电压提供。

2) 当电压低于100 V时，阴极区激发态氩原子密度保持在10⁹ 1/m³以下。随着电压的增加，由于二次电子的激发，氩原子的密度迅速增加，最终趋于稳定。

3) 当电压为25 V~100 V时，阴极区域的电子温度保持在2 eV~2.2 eV之间，然后有大幅度的升高。这是由于电场的增加和阴极的离子轰击，表明电子温度受电场的影响很大。

4) 距阴极区0.1 m，当电压为25 V~100 V时，由于电子离子的连续偏移，空间电荷呈现下降趋势，最终趋于稳定。可以看出，反应室中的氩气对空间电荷有很大的影响。

基金项目

激光与物质相互作用国家重点实验室开放基础研究课题(SKLLIM2113)，浙江省“尖兵”“领雁”研发攻关计划项目资助(2022C01094)。

参考文献