空心阴极羽流形态与热特性的变化机理

2022-11-21于博徐亚男康小录

中国空间科学技术 2022年5期

于博，徐亚男，康小录

上海空间推进研究所，上海 201112

1 引言

空心阴极是一种利用热电子发射机制为电推力器提供点火原初电子以及羽流中和电子的装置。电推进空心阴极的工作性能与其放电模式息息相关，在传统空心阴极的多次试验[1-4]中，空心阴极存在3种放电模式，而不同的放电模式对应的阴极性能及放电状态有所不同。

对于任何工质气体的阴极来说，都存在不同放电模式[3-4]：点状模式的弧光主要集中在触持极出口位置，呈现亮斑状形态；在羽状模式下，放电弧光呈现分散形态；弥散模式是指放电弥散到整个舱体内的状态。一般地，电推进应用中的空心阴极放电模式通常指点、羽两种模式。阳极电压及放电振荡的变化通常视作放电模式转变的信号，根据相关报道，在点状模式和羽状模式的过渡过程中会出现阴极的最佳工作状态，在该状态下，发射体温度和放电电压会回归到极小值[5-6]，是阴极工作的最高能效状态。上述现象与发射体材料、阴极结构尺寸不相关。无论是六硼化镧阴极[7]、钡钨阴极[2]、C12A7阴极[8]或是阴极推力器[9]，都至少存在点、羽两种放电模式。在阴极的不同放电模式中，阳极电压、阴极温度分布均有不同程度的变化[10-11]。对于一般钡钨空心阴极而言，点状模式的阳极电压在20 V以下，但也有研究人员设计一种新结构的六硼化镧空心阴极，可将阳极电压降至15 V以下[12]。此外，羽状模式还伴随着更高的电磁干扰[13]，并且，羽状模式下的高能离子状态参数[14]以及电极的腐蚀率的增加[15-16]也是阴极性能衰减的共性问题。在2015年后，研究人员逐渐开始关注高比冲电推进技术，将推进工质的研究重点转移到Kr[17]，对于20 A级空心阴极来说，Xe的点-羽模式转换点约为10.5 ml/min，对应阳极电压12.8 V，而Kr的转换点在30.8 ml/min，并且在任何相同的阳极电压下，Kr所需要的气体流率都会高于Xe。进一步地，有学者针对Kr工质阴极的工作点进行了详细研究[18]，工作参数覆盖阳极电流0.1～1 A，气体质量流率1.1～2.6 ml/min，获得关于Kr阴极的工作模式矩阵：点状模式主要分布于阳极电流高于0.6 A及气体流率高于2.2 ml/min的范围内。

目前为止，关于空心阴极放电模式的研究多集中于阴极各类放电特性和参数规律的归纳，而对于阴极结构设计而言，放电模式转变机理的研究才是重点：Myron M提出点状模式可以提供足量的正离子以减少等离子体对触持极流入电流，而羽状模式则需要触持极提供电子加速鞘层来维持放电[3]。 Salhi A认为当阴极表面的返流电荷达到较高数值时，阳极电压可能达到最小值，并且这种低能耗状态是等离子体流场的供能较高的结果[19-20]。以上两方面机理研究都能够从试验与数值模拟方面得到相关验证，但迄今为止，研究人员对点-羽两种模式的转变机理研究依然不够深入，鲜有针对阴极结构优化设计方面的启迪。为此，本文将对上述问题进行进一步研究。

空心阴极放电模式的转变过程与内部微观过程有密切关系，具体涉及粒子碰撞过程以及等离子体参数分布，本文采用比较适合捕捉阴极内部粒子流动与碰撞过程的混合模型——单元粒子模型/化学动力学碰撞模型(PIC/PCD)[21-24]来实现阴极内部等离子体流动过程模拟。接着，通过开展相同结构、工况的空心阴极放电试验，以羽流光色的计算图像与试验照片进行比对来验证数值模型，即从放电模式捕捉精度的角度来验证模型。在此基础上，对不同气体流率、不同阴极顶孔径下的阴极内部过程进行计算(计算工况覆盖点、羽两种模式)，获得等离子体参数的空间分布以及各类能耗数据的变化规律，以揭示阴极放电模式转变的内在机理。

2 数值模型

根据空心阴极的几何结构特点，本文采用2D计算域来表征阴极内部流场，所涉及等离子体流场主要包含3个物理过程：1)粒子在电场中的输运过程；2)粒子之间的碰撞过程；3)粒子与壁面的碰撞过程。具体计算域网格划分如图1所示，电子由发射体表面发射进入计算域，而电子、原子在计算域中发生各类碰撞(电离、激发和弹性)，由电离碰撞而生成离子，离子将参与和相关壁面的碰撞过程。

图1 空心阴极的2D计算域Fig.1 A 2D computational domain of the hollow cathode

给出阴极内部数值模型的一些假设：1)中性原子处理为背景气体，即利用第三方商用软件(例如稀薄流体相关计算软件)来单独计算原子数密度分布，并将该计算结果导入计算程序，参与电子-原子的碰撞计算；2)忽略重粒子之间的碰撞过程，原因为该类碰撞过程的平均自由程远大于电子-原子碰撞，或者该类碰撞对整体粒子能耗的计算精度影响较小；3)仅考虑离子-壁面碰撞，而忽略电子-壁面碰撞，原因为电子-壁面碰撞对整体电子能耗的影响较小。采用上述假设的目的在于简化模型，提高计算速度。

数值模型由两部分计算组成：等离子体流场计算和阴极组件热分布计算。流场计算为热计算提供输入条件，即通过流场计算获得表面1、2和3的热流密度分布，以该参数分布作为热计算的第二类边界条件，再通过热计算所求解的温度来更新流场计算中的发射体温度，具体流程可参考图2。

图2 数值模型的流-热迭代计算流程Fig.2 A flow chart of the flow-heat iterative calculation of the simulation model

2.1 粒子输运过程

粒子输运过程主要涉及电子和离子两种粒子，电子是从计算初期就进入计算域的粒子，而离子则是在电子输运过程中与原子发生电离碰撞所产生的，属于计算域中的再生成粒子，因此，离子的初始位置是不确定的，但初始速度均为0。

针对电子、离子两种粒子的模拟过程，本文采用单元粒子(PIC)模型来处理粒子追踪以及参数节点分配的问题。关于PIC的算法细节可参考文献[23]，故这里仅介绍本文的设定：1)Xe原子作为背景气体处理，其数密度分布参与各类碰撞计算，并统计当地网格内的原子消耗。2)采用双粒子PIC设定(非粒子加速法)，即电子、离子均参与输运、碰撞计算，且具有不同时间步长。离子时间步长2×10-9s，电子时间步长1×10-11s，这样的设定可以提高对阴极鞘层电势的模拟精度。但是，这也意味着计算量增多以及计算时间变长，因此需要选择一种计算速度快的碰撞模型与之配合。3)计算粒子权系数(1个计算粒子代表多少个真实粒子)为1×109。

发射体的电子入口表面的电子发射数量可根据理查森发射电流密度公式[25]计算：

(1)

式中：发射系数A0=1.204×106A/(m2·K2)；透射系数D=0.98；Temit为发射体工作温度，可通过图2中的迭代计算结果确定，初始值设为1 200 K；φemit为发射体逸出功，对于钡钨411发射体，取2.08 eV；e为元电荷；E为鞘层电场强度；ε0为真空介电常数。

鞘层电势可以通过Child提出的鞘层模型[27]计算：

(2)

式中：Te为电子温度；Ti为离子温度；mi为离子质量；me为电子质量。

2.2 粒子间的碰撞过程

前文有述，粒子间所涉及的碰撞过程主要为电子-原子碰撞，本文提出一种求解电子-原子碰撞过程的数值模型——PCD(plasma and chemical dynamics)模型，常见于磁流体模型中。PCD与PIC模型的兼容需要进行一些处理。PCD模型可捕捉碰撞反应速率，该参数可转化为对网格内发生某类碰撞的次数，这样即可与PIC模型形成连接。由于PCD模型针对整个网格内的电子进行碰撞判断，而不是针对某一个电子，因而PCD模型比MCC模型具有更高的计算速度，并且计算结果更加光滑且规律性明显。

在一个网格内，发生第j种(j参考表1中的编号)碰撞反应的次数为：

Nreaction,j=Sreaction,jVmesh

(3)

式中：Sreaction,j为第j种碰撞反应的反应密度速率；Vmesh为面网格所对应的环形体网格体积，为对应面网格面积与环周长的乘积。

根据表1所提供的化学反应系数ξj，可以获得Sreaction,j的计算方法：

Sreaction，j=nennξj

(4)

式中：ne为当地网格的电子数密度；nn为原子数密度。

表1给出了本文所考虑的电子-原子的碰撞反应类型，包括Xe原子典型的8条能级结构、一阶和二阶电离能级(小功率霍尔推力器的放电电压较低，三阶及以上的电离、重粒子引发的7 s、8 s、5 d、6 d之间的跃迁可以忽略，但二阶电离以及亚稳态电离不可忽略)以及相应的电子能量阈值和化学反应系数，而忽略了库伦碰撞、离子复合等放电室内部自由程较高的碰撞类型。当整个网格内的所有电子进行各类碰撞反应数量得到计算后，每一类碰撞发生的具体电子则可以在超过对应反应阈值的电子中进行随机选取，以确定每一个电子的碰撞数据。

表1 e-Xe碰撞反应的相关数据

对碰撞后的粒子运动状态进行设定：1)电子发生弹性碰撞后，原子运动状态视为不变，电子能量不变，但电子运动方向发生散射变化；2)电子发生激发碰撞后，原子运动状态不变，在激发态寿命末期释放相应能量的光子(光子仅在计算羽流光色中使用，不参与粒子计算过程)，电子能量减少激发时的跃迁能级，运动方向发生散射变化；3)电子发生电离碰撞后，原电子、新生电子与新生离子的能量之和等于原电子发生电离前的能量与电离能级的差，新生离子能量与中性气体背景气体的温度一致，方向随机分布，而原电子和新生电子的运动方向依然遵循电子散射变化规律。

在发生碰撞后，电子与碰撞前运动角度的夹角称为散射角，可描述为：

(5)

式中：Ek，e为碰撞前原电子的能量；R为随机数。假设碰撞前原电子运动方向与Z轴夹角为θ，那么，碰撞后的电子运动方向与Z轴的夹角θ′描述为：

cosθ′=cosθcosχ+sinθsinχcosφ

(6)

式中：φ为回转方位角，在360°内随机选取，即φ∈(0，2πR)。

2.3 粒子与壁面的碰撞过程

在阴极放电的等离子体流场中，各种粒子与各个壁面的碰撞十分频繁。若考虑所有的壁面碰撞类型，计算效率较低，因此，对粒子与壁面碰撞进行简化处理。

首先，电子与任何壁面的碰撞在趋势上均服从能量守恒过程，可视作镜反射。其次，除负电表面外，离子与其他各壁面的碰撞视作镜反射过程。

离子与负电表面的碰撞是最为复杂的。由于鞘层的存在，并不是所有正离子都会与发射体表面发生碰撞[25-26]：当待进入鞘层的离子速度低于波姆速度时，认为离子无法进入鞘层，以镜反射来处理碰撞过程；当离子速度不低于波姆速度时，认为离子可以进入鞘层。在鞘层的加速下，离子会产生较高的能量，并将能量全部沉积在发射体表面。发射体下游前端到阴极孔顶末端的壁面热流沉积密度分布由式(7)计算：

(7)

式中：Ekinetic为离子轰击壁面前的动能；Ti为离子温度；Eionization为Xe电离能；φ为鞘层电势，可参考式(2)求解；S为离子轰击发射体表面的(i，j)面网格的数量；A为(i，j)面网格的面积；Δt为时间步长。

2.4 热分布计算模型

以发射体表面和阴极顶表面的热流密度分布为第二类边界条件，环境温度300 K为第三类边界条件，可对整个空心阴极进行稳态的热分布计算。网格间的热传导采用傅里叶导热模型，而空间内的热辐射换热并不采用传统的四次方定律(本文工况中的表面数量较多，单个表面的温度分布也有差异，利用角系数计算的热流密度误差较高)，取而代之，采用辐射传递系数来表征空间内网格之间的辐射传递角度和比例，该系数定义为由表面A出发的、经过多次反射达到表面B而被表面B吸收的能量，与表面A所辐射的总能量的比例。具体可通过追踪光路的算法来实现，例如，网格m对网格n的辐射传递系数Rm→n为：

(8)

式中：Nm为从网格m辐射出的全角度光线数量；Nm→n为网格n接收到的来自网格m的光线数量(这些光线包括直接到达的，也包括经过反射到达的)。

网格m所吸收的净辐射热流密度qrad,m为：

(9)

式中：εm为网格m的表面发射率；Am为网格m的面积；Tm为网格m的温度；QT，m为网格m所吸收的总辐射热流，可表示为：

(10)

将热辐射所求解的网格内的热流密度作为第二类边界条件代入导热计算中，可求解出空心阴极组件内的温度分布，进而可输出发射体的平均工作温度Temit。

2.5 验证试验

为验证和修正数值模型，在直径为1.2 m、长2 m的真空舱内开展空心阴极的放电试验。由于试验测量方法的限制，关于阴极内部等离子体参数的验证难度较大，因而本文仅针对与阴极放电模式相关的宏观物理参数进行对比验证：阴极表面测点温度(定量验证)和阴极羽流光色(定性验证)。

整个试验平台布置见图3(a)，抽空系统包括三级抽空机组，极限真空度可达到6×10-6Pa以下，在空心阴极工作时，真空度可以达到2×10-3Pa以下，能够满足试验需求。空心阴极采用额定放电电流为4.5 A的原理试验件(LYJ-4.5)，如图3(b)所示。发射体采用钡钨411材料，表面逸出功为2.08 eV，阳极板材料为Mo，距离阴极触持极顶2 cm。测温采用S型铂热电偶，将4个热电偶固定在触持极外壳表面，该热电偶含铑10%，含铂90%，长期(高于10 s)测量的最高允许温度为1 300 ℃，短期测量的最高允许温度为1 600 ℃。此外，以光学温度计对阴极顶表面进行测温，具体测量图像见图3(c)。

图3 数值模型验证试验系统Fig.3 The systematic diagram of verification test for the simulation model

(1)温度验证

图3(a)中，4个热电偶对应测点编号分别为1～4，阴极顶表面的测点为编号5(见图3(c))。表2给出5个测点对应的计算值与试验值对比以及相应的误差。

通过阴极5个测点的温度计算值与试验值对比，可知试验值普遍低于计算值，这与本文数值模型的接触热阻设定、导热系数与真实值差异有关，但考虑误差在5%以内，对后文的机理分析影响可以忽略，因此认为该计算模型的精度能够满足研究需要。

表2 LYJ-4.5阴极的温度计算值与试验值对比

需要说明的是，温度参数是表征阴极工作状态的重要参数，是在流场与热场迭代计算中同时涉及的物理参数，可以在一定程度上体现整个数值模型的精度。然而，在现阶段，仅能从间接角度通过温度来验证流场计算精度，而在未来测量方法得到进一步发展后，有望直接对等离子体流场的计算精度进行验证。

(2)光色图验证

通过本文数值模型针对阴极流场的原子激发情况进行计算，在网格内统计发生激发碰撞的次数，见图4～图7。

图4 光色图和试验照片的对比(2.5 ml/min，羽状模式)Fig.4 Comparison of color picture and test photo(2.5 ml/min，plume mode)

图5 光色图和试验照片的对比(3.0 ml/min，羽状模式)Fig.5 Comparison of color picture and test photo(3.0 ml/min，plume mode)

图6 光色图和试验照片的对比(3.5 ml/min，最佳工作点)Fig.6 Comparison of color picture and test photo(3.5 ml/min，best running point)

图7 光色图和试验照片的对比(4.0 ml/min，点状模式)Fig.7 Comparison of color picture and test photo(4.0 ml/min，spot mode)

羽状模式下的羽流光色更偏向品红色(波长380～420 nm)，而点状模式下的羽流光色则更偏向蓝紫色(波长440～460 nm)，在这一方面，计算结果与试验照片在定性上几乎一致。这说明，数值模型在捕捉激发能级的计算精度方面拥有较低误差，可以较精确地模拟阴极羽流状态，这为后文分析点-羽模式相关机理的有效性提供了验证。

3 结果与分析

根据相关文献所述，空心阴极放电模式的转变现象主要为：点状模式中，阳极电压较低，阴极顶温度较低，并伴随亮斑状羽流；羽状模式中，阳极电压较高，阴极顶温度较高，羽流较为发散。本文将从现象入手来分析这种现象的内在机理。计算所涉及工况参数见表3，其中，尺寸a为阴极顶孔半径，尺寸b为阴极顶与触持极间距，尺寸c为触持极孔半径。在每个工况后的括号中还给出由计算数据推断的阴极放电模式。

首先，讨论气体流率对阴极温度的影响规律。图8共给出6组工况的温度分布结果：由于阴极放电电流均为4.5 A，故所有工况的发射体温度几乎相同，根据肖特基效应下的零场发射电流公式反推计算，发射体温度在1 450 K左右，但空心阴极发射体始终处于“空间电荷场”状态(共3种状态：拒斥场、空间电荷场和加速场)，由此，3个工况的发射体温度在数值上接近，但不会完全相同。值得注意的是，6组工况的阴极顶温度关于气体流率有所不同，例如，case 1的阴极顶温度高出case 3近80 K，而case 4的阴极顶温度比case 6高出110 K，这说明阴极热流输入表面的热流密度特性有所变化。

其次，讨论阴极结构变化对温度变化的影响规律。通过对比相同流率下、不同结构阴极的温度分布，可见，随着阴极顶孔径的增加，阴极顶大于发射体的温度差值会逐渐上升。进一步地，阴极组件温度分布与相关表面的热流输入数据相关，因此，本文将继续考察阴极内部离子数密度的分布情况，以从中获得温度分布变化规律的诱因。

表3 LYJ-4.5阴极各工况结构及气体流率

图8 阴极温度分布云图(case 1～case 6)Fig.8 Contour of cathode temperature distribution(case 1～case 6)

图9给出离子数密度分布，同时标注了流线和各关键表面的离子碰撞热功率(括号内为各表面流入热功率的百分比)。由计算结果可知，随着气体流率的减小，每种结构阴极表面热功率都在逐渐从阴极顶向发射体过渡。这意味着，低气体流率工况的离子-壁面碰撞热量更多地贡献给了阴极顶，而高流率工况的离子-壁面碰撞则更多地贡献给发射体。该机制恰好为图8中6个工况阴极顶温度存在差异阐明了诱导原因。而产生离子这种变化规律的原因在于低流量工况中的离子数密度分布在阴极顶外部(出口的下游)，而高流量工况的离子数密度多分布于阴极顶孔内，并且非常集中，这可以极大程度向发射体腔内输送离子。为揭示离子这种分布特性的内在机理，需要考察离子产生位置(电离核心区)的相关数据。电离核心区不仅与离子数密度峰值区相关，还主要取决于电子数密度峰值区和空间电势的骤降区。

图9 离子数密度分布云图(case 1～case 6)Fig.9 Contour of ion number density distribution(case 1～case 6)

续图9Fig.9 Continued

为分析离子数密度在孔区的分布特性内在机制，给出等离子体流场的空间电势分布计算结果，见图10。首先，case 1、case 2、case 4和case 5的放电电压较高，均在20 V以上，而case 3和case 6的放电电压较低，上述结果与点-羽工作模式下的阳极电压特性相符。此外，case 1、case 2、case 4、case 5与case 3和case 6相比较，前四者的空间电势在阴极顶孔区外部的数值较case 3和case 6升高，这会导致电子温度在孔区外将达到较高数值，可在孔区外触发较多的电离碰撞，那么case 1、case 2、case 4和case 5的电离核心区主要分布在阴极顶孔区外部，而case 3和case 6的电离核心区应该位于阴极顶孔内部。上述结果可以解释点-羽模式下电离核心区发生变化的原因。

图10 阴极空间电势分布云图(case 1～case 6)Fig.10 Contour of cathode temperature distribution(case 1～case 6)

为深入剖析点-羽模式中电离核心区从顶孔内部向外部推移的原因，需要继续考察两方面因素：1)电子数密度，该参数是指征电离核心区位置的参数；2)电子温度和中性粒子数密度，两者是决定电离碰撞概率的因素，可以进一步指向电离核心区。

图11为电子数密度分布，case 1、case 2、case 4和case 5的电子数密度峰值区均位于孔外，而case 3和case 6的电子数密度峰值区位于孔内部，这一点与电离核心区在点-羽模式过渡中的变化特性相符。

图11 阴极电子数密度分布云图(case 1～case 6)Fig.11 Contour of electron number density distribution(case 1～case 6)

图12为电子温度分布计算结果，图13为背景气体压强分布结果。结合上述两个计算结果，case 1、case 2、case 4和case 5的电子温度骤降位置在阴极顶孔出口位置，分布较发散，该位置可指向这些工况的电离核心区位形(推断可能为羽状模式，羽流光色发散不集中)，而case 3和case 6的电子温度骤降位置位于阴极顶孔内，而分布较集中，该位置可指向这两个工况的电离核心区位形(推断可能为点状模式，羽流光色集中)。

图12 阴极电子温度分布云图(case 1～case 6)Fig.12 Contour of cathode temperature distribution(case 1～case 6)

无论是气体流率变化或是阴极顶孔半径变化，压强分布可以解释电离区位置变化的原因：case 1、case 2、case 4和case 5的阴极顶孔区内部压强较低(电离核心区的数密度均低于9.0×1022m-3，该数值是4.5 A放电电流下对应的原子数密度临界值，低于该数值时，电离碰撞概率将下降显著)，而case 3和case 6的阴极顶孔区内部压强较高，孔区内距离出口0.1 mm位置的原子数密度在9.5×1022m-3以上。对于羽状工况而言，阴极顶孔区数密度不足，导致电离率不够，难以维持孔区内部的高电离放电，对发射体直接轰击的离子数量就会降低，遂放电电压提高，令电子温度升高以提高更多的电离概率，使电离核心区从孔区内部推移到孔区外部。

图13 阴极原子数密度分布云图(case 1～case 6)Fig.13 Contour of atom number density distribution(case 1～case 6)

图14给出阴极内部原子发生激发碰撞的次数分布，其中激发碰撞最集中位置(白区)可定义为激发核心区，可见激发核心区与电离核心区的位置几乎一致(对比图14和图11)，并且点状模式(case 3和case 6)的激发核心区都位于阴极顶孔内，而羽状模式(case 1、case 2、case 4和case 5)的激发核心区大部分都位于阴极顶孔外。这个特性与电离核心区几乎一致，原因为激发与电离碰撞概率均是电子温度的函数且单调性几乎一致，两者有同增同降的相似规律。因此，在羽状模式下，电离核心区向阴极顶孔外推移就指征着激发核心区同样向孔外推移，导致退激位置的外推和发散，令羽流表现出发散的羽状形态。同理，点状模式下的阴极羽流会表现出较为集中的点状形态。至此，阴极放电模式的工作特性与内在物理机制之间就建立了联系。