鹿 畅，夏广庆，孙 斌，韩亚杰，关思琦

(1. 大连理工大学航空航天学院，大连 116024； 2. 工业装备结构分析国家重点实验室，大连 116024；3. 辽宁省空天飞行器前沿技术重点实验室，大连 116024； 4. 大连理工大学白俄罗斯国立大学联合学院，大连 116024)

0 引 言

随着微小卫星的蓬勃发展，微推进系统[1-3]的应用也愈加广泛。微推进系统可为微小卫星提供精准的姿态调整和轨道控制，因而对微小卫星的在轨应用具有举足轻重的作用。由于微小卫星的质量通常小于10 kg[4-5]所以要求其推力器必须满足体积小、重量轻、功率低、比冲高、推力器范围宽等要求[6]。传统化学推进或冷气推进无法很好的满足上述需求。因此具有高比冲、长寿命、高控制精度等特点的微型电推进系统成为最有竞争力的选项。研发高性能的微型电推进系统对微小卫星的进一步发展和应用至关重要。微推进系统作为微纳卫星的核心部件，可以为微纳卫星提供精准的姿态调整和轨道控制。碳纳米管阵列推力器[7-9]是新近提出的一种基于碳纳米管(Carbon nanotubes, CNT)场电离的微型离子推力器，可将离子推力器的尺寸缩小至1 cm。

碳纳米管阵列推力器的工作原理示意图及实物图如图1所示。碳纳米管阵列推力器利用碳纳米管的场电离效应在纳米管尖端高效的电离工质气体，然后利用栅极引出并加速离子。相较于其他类型的微电推进系统，碳纳米管阵列推力器具备易小型化、高效率、高比冲等特点。根据碳纳米管电极的极性还可以切换发射离子和电子，不需要中和器，所以也具有较长的使用寿命。此外，碳纳米管阵列推力器还可以结合冷气推进，作为冷气推进尾端的二次加速装置，并避免场发射推力器和胶体推力器中的浸润污染问题[7]。

图1 碳纳米管阵列推力器示意图及实物图[8]Fig.1 Schematic diagram of carbon nanotube array thruster[8]

目前有关碳纳米管阵列推力器的研究主要集中在实验验证。Hicks等[8]首次制造和验证了碳纳米管阵列推力器的可行性。同时，Hicks在实验中发现CNT场离子发射电流对工质流率和阴阳极电压非常敏感，阴阳极电压过大会导致场电离过程不稳定。2017年，Tajmar等[7]基于碳纳米管阵列推力器和冷气推力器为欧空局LISA任务研制了1 cm2MEMS离子推力器。该推力器在0.06 sccm流率的氩气下展示出了良好的性能，其极限电压为740 V，最大束流1.2 mA。电离效率约28%，比冲1700 s，总推力接近30 uN。

评价碳纳米管阵列推力器性能的一个重要指标是CNT的场增强能力。场增强能力越高，同等电离效率下电极所需施加的电压差越低。当电极施加的电压差过高时，CNT尖端的放电模式会进入其他放电模式，如电晕放电等，从而进入不稳定放电过程，严重时会发生击穿短路。这些情况不仅会影响离子引出及加速过程，还会大大降低推力器的使用寿命。因此，对于碳纳米管阵列推力器，CNT的场增强能力越高，则表明其性能越好。然而，CNT场增强能力与推进工质流率、放电电压、碳纳米管几何特性等参数均密切关系。然而上述文献均尚未针对这些问题展开研究。因此，为进一步优化提高碳纳米管阵列推力器的性能，亟需开展其放电特性研究。

鉴于数值模拟方法高效、经济的特点，本文针对碳纳米管阵列推力器的放电特性开展了数值模拟研究。本文首次基于直接蒙特卡洛(DSMC)算法[10]、浸入式有限元算法(IFE)[11-15]、粒子云及蒙特卡洛碰撞(PIC-MCC)算法[11-15]以及Ammosov-Delone-Krainov(ADK)[16]隧穿电离模型建立了碳纳米管阵列推力器工质气体和CNT场电离过程的三维仿真模型。由于采用了IFE算法，因此该模型可以较好的处理CNT的复杂几何结构。本文首先利用文献中的实验数据验证了模型的正确性，然后针对碳纳米管阵列推力器的相关放电特性进行了研究。

1 碳纳米管阵列推力器放电室三维仿真模型

本文仿真模型主要包括三个子模型：工质原子仿真模型、CNT场电离模型和离子运动碰撞模型。其中工质原子仿真模型采用DSMC算法进行求解。由于碳纳米管阵列推力器的电离效率较低(10%～30%)[7]，因此工质原子的求解过程中忽略了电离过程对工质原子分布的影响。CNT场电离模型采用了ADK隧穿电离模型。离子的运动及碰撞采用了IFE-PIC-MCC算法。下面分别对上述模型进行介绍。

1.1 基于DSMC算法的工质原子仿真模型

在电推进中，推进工质通常采用Xe、Ar、Kr等惰性单原子气体，因此本文拟基于可变硬球假设(VHS)建立DSMC仿真模型，且模型中只考虑两体碰撞，不考虑多体碰撞。下面对DSMC算法的主要内容进行介绍。

粒子碰撞频率及自由程：

在DSMC算法中需要首先给定粒子的初始状态、时间步长、网格长度等初始输入参数。其中，时间步长和网格长度是决定计算精度两个关键参数之一。这两个参数是根据粒子的碰撞频率及平均自由程确定的。其中，时间步长要求小于粒子最大碰撞频率的倒数，而网格长度要求小于粒子的平均自由程，这样才能保证所有碰撞可以被有效捕捉。粒子的平均碰撞频率如式(1)所示：

(1)

粒子的平均自由程λ如式(2)所示：

(2)

(3)

粒子碰撞采样：

DSMC算法的第二个主要计算内容是根据粒子的碰撞频率来随机抽取碰撞对。本项目拟采用非时间采样碰撞方法(NTC method)。假定某任意网格内有N个粒子，则NTC方法的基本过程为：首先令Δt表示时间步长，则Δt内，某粒子在该网格内发生碰撞的概率PD可以表示为：

(4)

式中：c表示该粒子的运动速度，VC表示网格体积，S表示碰撞截面。

那么在该网格内发生碰撞的最大概率(PD)max可表示为：

(5)

式中：下标max表示取最大值。

因此，在该网格内所需选择的最大碰撞粒子对的数量NNTC为：

(6)

最后，在该网格内随机抽取NNTC对粒子进行碰撞处理。

粒子碰撞处理：

由于本文拟采用VHS模型且仅考虑粒子间的两体碰撞，所以粒子间碰撞仅考虑两体弹性碰撞，在碰撞过程中无能量损失。则粒子碰撞过程可根据两个粒子体系的动量守恒和能量守恒进行描述。

最后，DSMC方法的计算流程如图2所示。

图2 DSMC算法流程图Fig.2 The flowchart of DSMC method

1.2 基于ADK隧穿电离模型的CNT场电离模型

Bauer和Mulsner[16]总结了数个用于描述场感应电离过程的理论推导模型，并通过与描述强外电场下的单电子库伦势时域薛定谔方程(TDSE)的数值解对比发现，当外加电场不大于临界电场[17]时，氢的ADK模型与TDSE数值解吻合最好。当外电场大于临界电场时，CNT电离机制由隧穿电离变为势垒抑制电离(BSI)。此时，ADK模型仍然可以与BSI电离模型吻合良好。

Bruhwiler等对ADK公式进行了简化[17]。简化后的ADK模型电离几率预测公式如式(7)所示：

(7)

式中：ξi为电子基态能量，n*为有效主量子数，Γ为标准伽马函数，E为电场强度。

1.3 基于IFE-PIC-MCC算法的离子运动碰撞模型

近几年发展出的将浸入式有限元(IFE)方法[18-20]与PIC-MCC方法相结合形成的IFE-PIC-MCC[21-23]方法是一种可以在结构化的网格中对等离子体进行仿真的较好的数值模拟方法。IFE方法在被界面分割的单元中使用了分片基函数的技巧，从而使得网格的划分不依赖于界面。而其他未被界面分割的单元则使用传统的有限元基函数。因此，IFE方法可以在结构化网格中有效的求解具有复杂边界条件的电磁场。而且，IFE方法相较于传统的有限元方法改动较小，因而易于编程实现。当计算区内无物体时，IFE算法还可以自动退化为普通FE算法。因此，本文采用IFE算法求解CNT电势及电场。

PIC算法的计算流程如下：首先，确定离子的初速度。假定场电离后，离子继承原子的速度，因此其初始速度符合麦克斯韦分布，如式(8)所示。

(8)

电子假定为全部被CNT吸收，因此在仿真区域内不考虑电子。因此，电势分布的泊松方程如式(9)所示：

(9)

上述方程由IFE算法求解。离子在电场中的运动由式(10)描述：

(10)

式(11)用PIC算法求解。离子与原子弹性碰撞及电荷交换碰撞采用MCC算法求解。以Ar气为例，Ar+与Ar之间的弹性碰撞截面如下所示：

σel(εi)=k1+k2lnεi

(11)

式中：k1,k2为弹性碰撞系数，其中k1=1.189×10-18m2，k2=-9.815×10-20m2 [24]；εi为Ar+离子动能；σel为弹性碰撞截面。

CEX碰撞截面如下所示：

σCEX(εi)=k3+k4lnεi

(12)

式中：k3,k4为CEX碰撞系数，其中k3=5.921×10-19m2，k4=-4.611×10-20m2 [24]；σCEX为碰撞截面。

单位时间内，Ar+与Ar间碰撞频率可表示为：

μi=σnn(xi)vi

(13)

式中：vi为Ar+离子相对Ar原子速度，nn为原子数密度，xi为Ar+离子所在位置，σ为总碰撞截面。

则碰撞概率可表示为：

Pi=1-exp(-μiΔt)

(14)

式中：Δt为时间步长，Pi为碰撞概率。IFE-PIC-MCC方法的具体计算流程可参考文献[11-15]。

2 仿真校验

本节通过对比文献[8]中的数据对本文提出的模型进行校验。文献[8]中采用的工质气体为Ar气。仿真所需的基本输入参数如表1所示。另外，本文仿真结果均匀无量纲值，无量纲参考量如表2所示。ADK模型所需的输入参数如表3所示。

表1 基本输入参数Table 1 The basic input parameters

表3 ADK模型输入参数Table 3 The parameters for ADK model

本节选取两根CNT阵列进行仿真，计算区域几何形状如图3所示。然后，基于两根CNT的仿真结果估算碳纳米管阵列推力器的总发射电流。最后通过对比不同放电电压下的结果对本文模型的正确性进行验证。

图3 仿真区域示意图Fig.3 The simulation domain

本节分别选取阳极(CNT)电压为350 V、400 V、450 V、500 V、550 V的五个工况进行仿真。仿真空间步长的无量纲值设为1。由于文献[8]中未给出CNT的密度，本文根据文献[25]中结果假定CNT间距为5 μm，因此整个计算区域大小设为20×5×50。为了满足CFL条件，每个工况下的无量纲时间步长Δt满足：

(15)

仿真区域的场边界条件设置为：Z=0及Z=Zmax面设置为Dirichlet边界，且Z=0面电势等于阳极电压，Z=Zmax面电势等于0 V。其余边界设置为Neumann边界。仿真区域的粒子边界条件设置为：Z=0及Z=Zmax面设置为吸收边界，其余边界设置为对称边界。

中性原子仿真结果如图4所示。可以看出中性原子在电离腔内基本呈均匀分布。

图4 Ar原子密度分布Fig.4 The density distribution of Ar atoms

CNT电压为550 V时的电势、电场强度、离子密度分布如图5和图6所示。另外计算得到CNT尖端的无量纲场强约为2600。可以看出，CNT尖端处电场强度得到极大地增强，Ar原子可在尖端处被电离。通过统计Z=Zmax端面离子电流得到该工况下两根CNT产生的离子电流约为4.60 nA。其中，Z=Zmax端面电流统计公式如式(16)所示：

图5 电势和电场强度分布Fig.5 The distribution of potential and electrical field

(16)

式中：Ii为统计的离子电流,M为统计的迭代步数,N为M步内流过Z=Zmax端面的仿真离子总个数,w为仿真离子权重。

其他工况下的离子电流如表4所示。

表4 不同工况下的CNT离子电流Table 4 CNT ion current under different working conditions

根据CNT的密度及推力器总面积可估算出总离子发射电流。根据仿真结果估算的离子电流与实验测量的离子电流对比结果如图7所示。可以看出，仿真结果小于实验结果，但两者趋势基本一致。导致这一结果的原因主要有两个方面：首先，CNT的密度基于假设给定，实际工况中CNT的密度可能远大于仿真工况。其次，CNT的直径只有2 nm，精确的计算CNT尖端场强较为困难，仿真所得电场强度仍然是实际情况的近似。

图7 实验结果及仿真结果对比Fig.7 Comparison of experimental results and simulation results

3 碳纳米管阵列推力器放电特性分析

本节首先对碳纳米管阵列推力器放电性能的评估方法进行了概述，然后重点分析了离子与原子间的碰撞对碳纳米管阵列推力器放电性能的影响。

3.1 碳纳米管阵列推力器放电性能的评估方法

碳纳米管阵列推力器的放电性能主要通过CNT的场增强能力评估。CNT的场增强能力通常通过其场增强因子γ表示。场增强因子γ定义为CNT几何增强因子β与电极间距λ的乘积：

γ=βλ

(17)

在不考虑电离的情况下，CNT几何增强因子等于其长径比，长径比越大CNT的几何增强因子越大。但是，电离过程产生的离子存在空间电荷效应会影响几何增强因子的大小。因此，稳态放电情况下的CNT几何增强因子需要通过计算Fowler-Nordheim场发射曲线的斜率得到。Hicks等将Fowler-Nordheim场发射曲线变换为如下方程[8]：

(18)

式中：x为电极电压倒数；y为x的函数，y=ln(Ix2)，I为总离子电流；φ为CNT表面功函数，约为5 eV；r为电极半径；A,B分别为两个系数，A=1.541×10-6A·eV/V2，B=6.831×109V/(eV3/2·m)。

式(18)为y=kx+b形式的线性方程，所以若已知电极电压及离子电流，即可求出方程(18)的斜率，进而反推CNT的几何增强因子，并最终求得CNT的场增强因子。

基于第2节中表3的结果可以估算得到CNT的几何增强因子约为2.6×107，场增强因子约为1308。

3.2 碰撞过程对CNT场增强能力的影响

由于电推进通常采用Xe(氙气)作为推进工质，因此本小节以Xe作为推进工质，分析碰撞过程对CNT场增强能力的影响。表5所列为Xe的输入参数，其余未列出的参数参考表1和表2。

表5 Xe原子输入参数Table 5 The parameters of Xe atoms

首先对比考虑碰撞与不考虑碰撞情况下的电势、电场强度和离子密度分布，如图8所示为550 V工况下的仿真结果对比。可以看出，考虑碰撞时，电场小于不考虑情况下，而离子密度却大于不考虑情况下。另外，表5对比了不同情况下的离子电流。根据表5中结果可以计算考虑及不考虑碰撞情况下的场增强因子分别为γ≈1665和γ≈1670。可见，考虑离子与原子间碰撞时，CNT的场增强能力有所降低。

图8 考虑碰撞(左侧)与不考虑碰撞(右侧)的仿真结果对比Fig.8 Comparison of simulation results considering collision (left side) and not considering collision (right side)

表5 不同情况下的离子电流对比Table 5 Ion current under different conditions

导致上述结果的原因为：碰撞过程阻碍了离子往栅极的运动，增加了离子在电离腔体内的密度，从而增大了离子的空间电荷效应，进而减弱了CNT尖端的电场，最终导致其场增强能力减弱。

4 结 论

本文基于DSMC、IFE-PIC-MCC算法以及ADK隧穿电离模型建立了碳纳米管阵列推力器工质气体和CNT场电离过程的三维仿真模型。通过与文献中的实验数据对比验证了本文仿真模型的有效性。针对碳纳米管阵列推力器在的放电特性进行了研究，并重点分析了离子与原子碰撞过程对碳纳米管阵列推力器放电特性的影响。结果表明：

1) 本文仿真模型所得仿真结果与实验结果趋势基本一致，但小于实验结果。分析认为主要由以下两个原因导致：首先，CNT的密度基于假设给定，实际工况中CNT的密度可能远大于仿真工况。其次，CNT的直径远小于网格尺寸，在仿真中只能按网格点赋值。所以模拟中无法捕捉到CNT实际几何特性对电场强度的影响，仿真所得电场强度仍然是实际情况的近似。

2) CNT可以极大地增强局部电场，稳态时场增强因子约为1308。

3) 碰撞过程对CNT附近电场及场增强因子均有影响，考虑碰撞时，电场及场增强因子均小于不考虑情况下，而离子密度却大于不考虑情况下。

4) 碰撞过程影响CNT放电的原因为：碰撞过程阻碍了离子往栅极的运动，增加了离子在电离腔体内的密度，从而增大了离子的空间电荷效应，进而减弱了CNT尖端的电场，最终导致其场增强能力减弱。