同轴微阴极电弧推力器的粒子网格法数值模拟

，，，,*

1.北京航空航天大学 宇航学院，北京 100191 2.北京控制工程研究所，北京 100094 3.空间物理重点实验室，北京 100076

同轴微阴极电弧推力器的粒子网格法数值模拟

熊子昌1，耿金越2，陈轩3，王海兴1,*

1.北京航空航天大学 宇航学院，北京 100191 2.北京控制工程研究所，北京 100094 3.空间物理重点实验室，北京 100076

采用粒子网格(Particle-In-Cell,PIC)方法对同轴微阴极电弧推力器(μCAT)工作过程进行了模拟研究，并应用自相似方法对模型进行简化，获得了推力器羽流区的电子数密度分布、离子数密度分布、电势分布及离子轴向平均速度，通过改变磁感应强度和位形分析磁场对推力器内等离子体运动特性及推力器性能的影响。计算结果表明，电子被外加磁场捕获约束在磁力线附近，低速离子与高速电子形成的双极扩散电场加速离子喷出；在相同流量情况下，磁感应强度0.02 T时，离子返流严重，磁感应强度0.05～0.30 T时，磁感应强度变化对速度影响较小；磁场位形对离子运动和推力器性能有较大影响，磁力线与轴线夹角较小时离子速度下降明显，夹角较大时离子返流严重。

同轴微阴极电弧推力器；粒子网格法；自相似法；磁感应强度；磁场位形

微阴极电弧推力器是美国乔治华盛顿大学近年来研制的一种新型推力器，主要有同轴和环形两种结构[1]，同轴微阴极电弧推力器结构示意如图1所示。其主要工作原理是电弧烧蚀阴极材料产生较高电离度的等离子体，在电场和磁场的共同作用下高速喷出而产生推力，外加磁场能够提高等离子体轴向速度,减小羽流扩散角,提高推力器比冲。该推力器结构简单，不需要复杂的推进剂供给系统，比冲较高，因而受到广泛关注，在微小卫星推进系统应用方面有广泛的应用前景[2-5]。

图1 同轴微阴极电弧推力器结构示意Fig.1 Schematic drawing of co-axial micro- cathode arc thruster

从目前已有的研究来看，推力器结构、磁感应强度及分布、脉冲放电能量的大小和持续时间、阴极材料的选取等都会对推力器的推力、比冲等性能参数产生较大影响[6-9]。但对该影响的物理机制和物理过程的深入研究尚无报道[10]。为了深入了解推力器结构与工作参数的改变对推力器性能影响，为推力器的优化设计提供参考，缩短推力器研制周期，有必要对微阴极电弧推力器工作过程进行系统的数值模拟研究。

本文采用粒子网格(Particle-In-Cell,PIC)方法对同轴微阴极电弧推力器进行模拟，应用自相似方法对模型进行简化，获得电子数密度分布、离子数密度分布、电势分布、离子轴向平均速度分布等结果，对微阴极电弧推力器的工作过程进行分析，估计推力器比冲；并且对微阴极电弧推力器影响最大的因素——磁场进行进一步分析，研究了磁感应强度和磁场位形对推力器性能的影响。

1 模拟方法

1.1 PIC法

PIC法[11]：通过向计算域内放置一些宏粒子来模拟真实的等离子体粒子的运动，每个宏粒子可以看作一起运动的等离子体团，其荷质比与真实粒子相等。

PIC方法的主要流程为：1)将网格内粒子的电荷分配到网格节点上获得电荷分布；2)在网格节点上求解泊松方程计算电场分布；3)在电磁场分布的基础上，按照离子所受电场和磁场作用力来求解宏粒子的运动状况。

因为微阴极电弧推力器中等离子体电离度很高，可以看作完全电离等离子体，所以，本文不考虑激发、电离等碰撞反应。微阴极电弧推力器中的外加磁感应强度远远大于放电所产生的自感磁场磁感应强度，所以认为推力器中磁场不变。在每个计算周期求解运动方程(式(1)(2))与泊松方程(式(3))以推动粒子运动和求解电场。

式中：m为粒子质量；v为粒子速度；E,B分别为电场强度和磁感应强度;x为粒子位置；φ为空间电势；ρ为空间电荷密度；ε为介电常数。

为了使得数值算法精确可靠，时间步长的选取依赖于等离子体德拜长度，空间步长的选取依赖于等离子体频率[7]。

式中：Δx为空间步长；Δt为时间步长；λD为等离子体德拜长度；ωpe为电子振荡频率；ωce为电子回旋频率。

1.2 自相似方法

在同轴微阴极电弧推力器二维模拟中，德拜长度为10 μm量级，推力器羽流特征长度为厘米量级，对于全尺寸的二维PIC模拟，模型计算域需要划分的网格节点数为百万量级，为保证模拟精度，所需要模拟的粒子数将为上亿个，这样庞大的计算量对于计算机的存储和运算速度是一个巨大的挑战。因此，模型的简化是很有必要的。

传统的模型简化方法有减少重粒子质量法和增大介电常数法，这两种方法可以有效地减少收敛所需要的循环步数和求解所需要的网格数，但是会破坏等离子体的物理过程，使结果有较大失真。对此，Taccogna等提出了自相似法，并验证了其可行性[12-14]。自相似法可以通过保证一些重要的物理、性能参数不变，以保证缩比后的模型物理过程相似，很大程度上保留了缩放后物理特性，故本文选择自相似法对同轴微阴极电弧推力器模型进行简化。但由于霍尔推力器的PIC模拟需要考虑电离过程，所以努森数为Taccogna自相似方法中的重要缩放参数。而本文中，因为等离子体电离度很高，假设为完全电离等离子体，不考虑电离过程，所以并不需要维持努森数不变。

(1)离子速度不变

(2)电流密度不变

故离子数密度也不变，推力器的粒子电流I=jA，A为面积。所以I～L2。

(3)带电粒子回旋半径和长度的比值不变

所以有B～1/L。

表1 各物理量对应的相似准则数Table 1 Physical quantities and their self-similarity factors

1.3 模型建立

微阴极电弧推力器放电过程中，外加磁场磁感应强度大小远大于放电过程所产生的自感磁场磁感应强度大小，因此认为放电过程中磁感应强度不变。

(1)磁场计算

磁场计算在商业软件COMSOL中的“磁场”模块下进行，根据推力器线圈的位置、结构和材料对推力器通道内及羽流区磁场进行有限元计算，得到磁场分布如图2所示。

图2 同轴微阴极电弧推力器磁场分布Fig.2 Distribution of magnetic field of the co-axial micro-cathode arc thruster

(2)计算域选取

为了尽可能保留推力器的羽流变化特性，本文选取0.06 m×0.03 m的矩形计算域，下游出口边界为零电势边界，轴线为轴对称边界，其余边界为开放边界。并在此基础上对模型采用自相似方法进行简化，缩放比为0.05。计算域网格数为300×150，计算时间步长为1×10-12s。

平均自由程为离子在两次碰撞之间平均运动的距离，平均自由程越大，则粒子间发生碰撞的概率越小[15]。假设种类1的粒子与作为背景的种类2粒子发生碰撞，则两种粒子碰撞的平均自由程为：

式中：n2为背景粒子的数密度；Q为粒子间的碰撞截面。

推力器电子温度为几电子伏特，电子-离子碰撞截面为10-18m2量级，推力器等离子体数密度为1018m-3量级，故电子-离子碰撞自由程为1 m量级，大于推力器特征长度毫米量级，因而本文中忽略电子-离子碰撞导致的复合。

由于推力器羽流电离度很高，羽流等离子体99%被电离[16]，近似完全电离，所以模型中，假设等离子体只由Ti二价离子与电子组成。额定工况下，磁感应强度为0.3 T，Ti离子的质量流量为3.8×10-6mg/120 μs，阴阳极之间的电势差为50 V，文献[17]表明，即使磁场会影响推力器的最终速度，但是它不会改变Ti离子在电极间的初始速度，这个速度为2×104m/s，本文将此速度作为模拟中离子的入口速度。

2 数值模拟结果

2.1模型验证

由于同轴微阴极电弧推力器尚无详细的试验数据报道，本文首先对不同磁场下的环形微阴极电弧推力器进行模拟，得到不同磁场及对应流量下离子平均速度，并与文献报道的试验测量结果[18]进行比较，如图3所示。模拟结果与文献报道试验结果符合较好，验证了模型的正确性。

图3 离子平均速度对比Fig.3 Comparison of average ion velocity

2.2 磁感应强度为0.3 T额定工况模拟结果

图4和图5分别为电子数密度分布和离子数密度分布。图中左下角黄色区域为推力器位置。电子质量很小，在磁场中做回旋运动的Larmor半径很小，因此被外加磁场捕获约束在磁力线附近，形成电子聚集区。而质量大、速度慢的离子，则由于双极扩散电场的作用，聚集于电子聚集区附近并加速喷出。图5可以观察到明显的沿电子运动区域的离子主流区。

图4 电子数密度分布Fig.4 Distribution of electron number density

图5 离子数密度分布Fig.5 Distribution of ion number density

图6为电势分布。由图6可知，磁力线捕获约束电子，并在此区域形成了狭长的低电势通道，电子和大部分离子均沿此通道运动。在图6电势分布中，沿轴线方向存在高电势区，这是因为阴极表面高电势与电子低电势区存在径向电场，使得阴极表面发出的离子在电子低电势区附近震荡，使得离子数密度在轴线较高，而此处并没有电子与之中和，从而产生局部高电势区。

图6 电势分布Fig.6 Distribution of electric potential

图7给出了离子轴向平均速度沿轴线变化。离子轴向平均速度通过统计各轴向位置所在径向平面内离子的轴向平均速度而得到。由图7中0.3 T条件下离子轴向平均速度可知，离子从阴极表面喷出，经过阴极表面附近的高电势区域速度急剧下降，之后进入低电势通道后速度逐渐增加至25.3 km/s左右,对应的比冲约为25 800 m/s。

图7 不同磁感应强度离子轴向平均速度Fig.7 Axial average ion velocity in different magnetic flux density

由以上仿真结果可知，磁场在微阴极电弧推力器中有至关重要的作用，将直接影响等离子体的运动过程及推力器性能。下面将分别改变外加磁场磁感应强度大小及位形，研究其影响规律。

2.3 磁感应强度大小影响

为了研究磁感应强度大小对推力器内等离子体运动过程及推力器性能的影响，本文还采用图2相同的磁场位形，仅同比减小磁感应强度，使阴极斑点附近磁感应强度分别为0.15 T、0.05 T和0.02 T，并分别对其进行模拟。结果如图7～图9所示。

图8和图9分别为电子数密度分布和离子数密度分布。电子的分布情况将通过双极扩散电场影响离子分布，由图8和图9可知，随着磁感应强度的减弱，磁场对电子的束缚能力逐渐减弱，导致电子聚集区的“宽度”逐渐增加；在磁感应强度为0.02 T时电子聚集区宽度进一步增加，甚至部分电子开始脱离磁力线的束缚，同时离子也出现了明显的返流现象。

图7为不同磁感应强度下离子轴向平均速度沿轴向变化。由图7可知，当磁感应强度为0.05～0.30 T时，随着磁感应强度变化，离子轴向平均速度变化趋势基本一致，且数值相当。这与图4中平均速度随磁感应强度增大而增大不同，主要是因为图4中，不同磁感应强度采用离子流量不同，造成了各磁场条件下离子速度的差异。而本文中为了分析磁感应强度单一条件对等离子体流动的影响，采用了相同流量。磁感应强度为0.02 T时离子轴向平均速度与0.05～0.30 T有较大差别，主要是部分离子出现返流引起的。

图8 不同磁感应强度下的电子数密度分布Fig.8 Distribution of electron number density in different magnetic flux density

总的来说，当磁感应强度为0.05～0.30 T时，随着磁感应强度变化，离子轴向平均速度变化趋势基本一致，随着磁感应强度的增大而小幅增大；当磁感应强度为0.02 T时，离子数密度降低严重，甚至出现返流的现象。

图9 不同磁感应强度下的离子数密度分布Fig.9 Distribution of ion number density in different magnetic flux density

2.4 磁场位形影响

除了磁感应强度大小，磁场位形也对等离子体运动及推力器性能有较大影响。

通过改变电磁线圈与推力器阴阳极之间的相对位置以改变磁场位形，即分别沿z轴正负方向移动3 mm，使得穿过阴极端面与绝缘介质端面交界点的磁力线发生变化，磁力线与z轴的夹角由额定工况下的11°分别变化为22°、2°，磁场与推力器相对位置示意如图10所示，图中蓝色、灰色、红色部分为推力器移动后的位置，带箭头的曲线表示磁力线。分别模拟这两种磁场条件下额定工况的结果，模拟中改变计算域中磁场位置，而推力器在计算域中的位置不变。

图10 穿过阴极端面与绝缘介质端面交界点的磁力线与轴线夹角示意Fig.10 Schematic drawing of 2 degrees between axis and the magnetic field line passing through the junction point of cathode end face and insulator end face

图11为磁力线与轴线夹角为22°和2°情况下的电子数密度分布。由图11可知，由于电子被磁化，电子运动轨迹与穿过阴极端面与绝缘介质端面交界点的磁力线一致。磁力线与轴线夹角度数越大表明磁力线越靠近z轴，导致电子运动也越靠近z轴。

图12为磁力线与轴线夹角为22°和2°情况下的离子数密度分布。由图12可知，离子主流区随电子运动轨迹改变而改变，即随着磁力线的改变而改变。当夹角很大时(22°)，磁力线与阳极很靠近，有一部分离子会与阳极发生碰撞，导致碰撞后的离子很有可能形成返流。当夹角很小时(2°)，磁力线与阳极距离很远，不但不会形成返流，还会减轻阳极附近的离子震荡幅度，形成稳定束流。但是因为磁力线位置远离z轴，造成离子主流区出口速度方向与z轴夹角变大。

图11 不同磁场位形电子数密度分布Fig.11 Distribution of electron number density in different magnetic field topology

图12 不同磁场位形离子数密度分布Fig.12 Distribution of ion number density in different magnetic field topology

图13为不同磁场位形情况下离子轴向平均速度沿轴向变化。由图13可知，当夹角很大时，虽然出口速度有明显提升，但是会造成离子返流现象。当夹角为2°时，离子主流区出口速度与z轴夹角变大，导致离子轴向平均速度减小明显，对推力贡献不利。而这3种工况之间相对位置仅仅相差3 mm。

图13 不同磁场位形离子轴向平均速度Fig.13 Axial average ion velocity in different magnetic field topology

3 结束语

本文采用2D轴对称粒子网格法对同轴微阴极电弧推力器等离子体运动过程进行数值模拟，为了加速算法的收敛时间，将一种可以缩小发动机几何尺寸，同时保证模拟时重要的参数不变的自相似法应用其中，以保证缩放前后结果一致。

计算结果发现：

1)推力器产生的电子被外加磁场捕获约束在磁力线附近，低速的离子由于双极扩散电场的作用，聚集于电子聚集区附近并加速喷出。

2)当磁感应强度为0.05～0.30 T时，随着磁感应强度的减弱，离子轴向出口速度逐渐减小，但总体趋势基本一致，当磁感应强度为0.02 T时，磁场对电子的束缚减弱明显，离子返流严重，出口流量降低明显。

3)磁场位形对推力器性能影响明显，穿过阴极端面与绝缘介质端面交界点的磁力线与轴线夹角较小导致轴向速度分量下降明显，而夹角较大则会导致返流严重。

为了进一步了解同轴微阴极电弧推力器系统的工作过程，今后的研究需要建立更加精细的推力器模型例如考虑壁面效应、阳极效应等。

References)

[1] KEIDAR M,ZHUANG T,SHASHURIN A,et al.Electric propulsion for small satellites[J].Plasma Physics & Controlled Fusion,2015,57(1):014005.

[2] ZHUANG T,SHASHURIN A,DENZ T,et al.Micro-vacuum arc thruster with extended lifetime[C].45th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit,Denver,August 2-5,2009.

[3] ZHUANG T,SHASHURIN A,HAQUE S,et al.Performance characterization of the micro-cathode arc thruster and propulsion system for space applications[C].46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit,Nashville,July 25-28,2010.

[4] HAQUE S,KEIDAR M,LEE T.Low-thrust orbital maneuver analysis for cubesat spacecraft with a micro-cathode arc thruster subsystem[C].33rd International Electric Propulsion Conference,The George Washington University,Washington D.C.,October 6-10,2013.

[5] SHASHURIN A,KEIDAR M.Comparative analysis of micro-cathode arc thruster performance[C].33rd International Electric Propulsion Conference,The George Washington University,Washington,D.C.,October 6-10,2013.

[6] ZHUANG T,SHASHURIN A,KEIDAR M.Micro cathode thruster (μCT) plume characterization[J].IEEE Transactions on Plasma Science,2011,39(11): 2936-2937.

[7] VAIL P,PANCOTTI A.Performance characterization of micro-cathode arc thruster (μCAT)[C].47th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit,San Diego,July 31-August 3,2011.

[8] HURLEY S,KEIDAR M.Linear actuated micro-cathode arc thruster system [C].52nd AIAA/SAE/ASEE Joint Propulsion Conference,Salt Lake City,July 25-27,2016.

[9] BRIEDA L,ZHUANG T,KEIDAR M.Towards near plume modeling of a micro cathode arc thruster[C].49th AIAA/ASME/SAE/ASEE Joint Propulsion Conference,San Jose,July 14-17,2013.

[10] ZHUANG T,SHASHURIN A.Circular periodic motion of plasma produced by a small-scale vacuum arc[J].Plasma Sources Science and Technology,2011,20:015009.

[11] VERBONCOEUR J P.Particle simulation of plasmas: review and advances[J].Plasma Physics & Controlled Fusion,2005,47(47):A231.

[12] TACCOGNA F,LONGO S,CAPITELLI M,et al.Plasma flow in a Hall thruster[J].Physics of Plasmas,2005,12(4):78.

[13] TACCOGNA F,LONGO S,CAPITELLI M,et al.Self-similarity in Hall plasma discharges: applications to particle models[J].Physics of Plasmas,2005,12(5):78.

[14] TACCOGNA F,SCHNEIDER R,LONGO S,et al.Kinetic simulations of a plasma thruster[J].Plasma Sources Science and Technology,2008,17(2):1-16.

[15] SZABO J.Fully kinetic numerical modeling of a plasma thruster [D].Cambridge: Massachusetts Institute of Technology,2001:92-106.

[16] KEIDAR M,HAQUE S,ZHUANG T,et al.Micro-cathode arc thruster for PhoneSat propulsion [C].27th Annual AIAA/USU Conference on Small Satellites,Logan,2013: 1-9.

[17] ZHUANG T,SHASHURIN A,BEILIS I,et al.Ion velocities in a micro-cathode arc thruster[J].Physics of Plasmas,2012,19(6):661.

[18] ZHUANG T.Micro-cathode arc thruster system for cube satellite[D].Washington,D.C.: George Washington University,2012:83-85.

(编辑：高珍)

Particle-in-cellsimulationofaco-axialmicro-cathodearcthruster

XIONG Zichang1，GENG Jinyue2，CHEN Xuan3，WANG Haixing1,*

1.SchoolofAstronautics,BeihangUniversity，Beijing100191,China2.BeijingInstituteofControlEngineering，Beijing100094，China3.ScienceandTechnologyonSpacePhysicsLaboratory，Beijing100076，China

A modeling study was performed to investigate the working process of a co-axial micro-cathode arc thruster with particle-in-cell (PIC) method,and the self-similar method was applied to simplify the simulation model.The distribution of electron number density,ion number density,electric potential and axial average ion velocity in the thruster plume region was obtained.The influences of magnetic field on plasma motion characteristics and thruster performance were studied by changing magnetic flux density and magnetic field topology.The results show that electrons are constrained by the applied magnetic field,cycling around magnetic field lines,and that ions are accelerated by the ambipolar diffusion electric field produced by slower ions and faster electrons.With the same mass flow rate,the ion reflux is severe when the magnetic flux density is 0.02 T,while the ion velocity is less affected by the magnetic flux density when the magnetic flux density is 0.05～0.30 T.Magnetic field topology has great influence on the ion motion and thruster performance.The ion velocity decreases obviously when the angle between the axis and magnetic field lines is small,while the ion reflux is severe when the angle is large.

co-axial micro-cathode arc thruster; particle-in-cell; self-similarity; magnetic flux density; magnetic field topology

http://zgkj.cast.cn

10.16708/j.cnki.1000-758X.2017.0081

V439

A

2017-05-04；

2017-07-24；录用日期2017-09-12；< class="emphasis_bold">网络出版时间

时间：2017-09-25 14:07:15

http://kns.cnki.net/kcms/detail/11.1859.V.20170925.1407.002.html

国家自然科学基金(11275021，11575019,11702021)；民用航天项目(混合模式推进系统优化技术)

熊子昌(1995-)，男，硕士研究生，xiongzichang@buaa.edu.cn，研究方向为空间电推进

*通讯作者：王海兴(1969-)，男，教授，whx@buaa.edu.cn，研究方向为空间电推进

熊子昌，耿金越，陈轩，等.同轴微阴极电弧推力器的粒子网格法数值模拟[J].中国空间科学技术,2017,37(5)：81-88.XIONGZC,GENGJY,CHENX,etal.Particle-in-cellsimulationofaco-axialmicro-cathodearcthruster[J].ChineseSpaceScienceandTechnology,2017,37(5)：81-88(inChinese).