考虑电荷交换的栅极区离子流数值模拟

刘 洋，李 娟，楚豫川，曹 勇

(1.哈尔滨工业大学深圳研究生院，广东深圳518055;2.兰州空间技术物理研究所，真空低温技术与物理重点实验室，甘肃兰州730000)

1 引言

离子推力器是利用电能电离、加速推进剂离子并高速喷出后，获得反推力的一种电推进装置。其比冲是常规化学火箭的数倍到数十倍［1］。与常规化学火箭相比，电推进系统具有比冲高、寿命长等显著特点，可有效减小推进剂的总质量，增加有效载荷，适应了空间飞行器的发展对推进系统提出的要求。因此，离子推力器系统越来越多地受到人们的关注，目前已成为我国航天领域的迫切需求［2］。

受电场等因素的限制，离子推力器的推力通常较小，因而，要完成一定的任务，离子推力器需要更长的工作时间。这就需要离子推力器有更长的工作寿命，一般要达到上万小时。离子推力器光学系统是推力器的关键部件之一，它决定离子推力器的性能和寿命。离子光学系统由栅极及其支撑结构组成，本文中离子光学系统栅极包括屏栅极和加速栅极。其失效形式主要有3种:屏栅极和加速栅极短路，栅极结构失效和栅极系统下游等离子体中的电子回流［3］。栅极的腐蚀是影响栅极寿命的主要因素。其中加速栅极下游区域电荷交换碰撞形成的电荷交换离子对加速栅极的溅射腐蚀最严重［4］。栅极长时间工作后，加速栅极可能产生穿孔、加速栅极孔变大等现象，严重影响离子推力器的工作性能和使用寿命。因此研究光学系统的腐蚀情况对于离子推力器的寿命评估有重大意义。

目前研究光学系统腐蚀主要采用实验和数值模拟2种方法。实验方法可以直接准确的获得真实的实验数据，但需投入大量的资金和时间，消耗巨大。与实验方法相比，数值方法结果直观，实验资金投入较小且大大降低了时间消耗，因此数值方法已被广泛采用。作者基于嵌入式有限元－质点网格［5～7］与蒙特卡洛［8，9］(IFE－PIC－MCC)的三维等离子体离子模拟程序，对离子推力器光学系统中离子密度分布及加速栅极腐蚀情况进行了研究。

2 离子运动计算模型设置

2.1 计算域和边界

离子推力器的栅极是很薄的电极板，栅极上开有数千个栅极孔，栅极孔呈六边形分布，具有对称结构。在计算模拟中，根据栅极的对称性，从栅极区域中选取包含栅极所有信息的最小区域进行计算。图1按照灰色区域将光学系统划分为长方形区域。计算域中包含了2个1/4的栅极孔，这样选取的计算域包含六边形分布的栅极孔的所有信息，也就是将相同的计算域拼接在一起后会形成完整的栅极排列。这样便可以认为所选取的计算域包含了整个栅极孔的信息。计算模型的参数都经过无量纲化。栅极几何参数等如表1所列。

图1 计算域二维简图

对于中性原子，只需考虑粒子运动边界，无需考虑电场条件。当中性原子越过上游边界或下游边界时，将该原子删除，如果中性原子到达4个侧面中的任意一个或栅极，则将该原子镜面反射回计算区域。

表1 栅极几何参数

2.2 CEX 的产生与加速栅极腐蚀

由于加速栅极下游表面的腐蚀主要是由交换电荷离子造成的，在离子束引出的过程中，如果不考虑离子与原子之间的碰撞，就忽视了交换电荷离子，只能计算得到部分束流对加速栅极孔壁的腐蚀情况，而无法获得加速栅极下游表面的腐蚀情况［10］。因此，要计算加速栅极下游表面的腐蚀，必须考虑离子之间的碰撞。离子之间的碰撞，主要包括离子与中性原子之间的电荷交换碰撞(CEX)，离子之间的库仑碰撞，中性原子之间的动量交换碰撞等。其中，中性原子与离子之间的电荷交换碰撞是加速栅极电流产生及栅极腐蚀的主要原因［11］。因此，在计算加速栅极腐蚀过程中，本文只考虑中性原子与离子之间的CEX碰撞，CEX碰撞是经过电场加速的高速离子和热运动的中性原子碰撞，产生慢速的电荷交换离子和快速的中性原子。即

作者用IFE－PIC方法模拟离子在离子光学系统中的运动，加入了蒙特卡洛［12］碰撞模块，来模拟考虑CEX离子在离子光学系统中的产生和运动。空间任意一个带电离子与空间分布的中性原子在每个步长Δt内的碰撞频率［12］为

碰撞概率为

式中 nn(xi)为离子位置的中性原子的数密度;vi为带电离子的速度;σT(vi)为离子电荷交换碰撞截面，采用经验公式［13］

式中 k1=－1.303×10－10;k2=30×10－10;v为带电离子与中性原子的相对速度，由于带电离子的速度远大于只有以热运动为主的中性原子的速度，所以本文用带电离子的速度来代替v。

在模拟计算中，离子在一个时间步长内运动的距离不超过最小网格长度的1/2，即

式中 vmax为计算域内模拟离子的最大速度;Δz为z方向的网格长度。

带电离子与中性原子的碰撞概率很小，但需要计算所有离子的碰撞概率与碰撞截面，非常浪费时间，所以本文采用零碰撞技术［13］来简化蒙特卡洛法的计算过程。首先，找出所有离子的最大碰撞频率νm

则可得到离子的最大碰撞概率为

假设在计算域内总的带电离子数为N，随机选择NPnull个离子，分别计算每个离子的碰撞频率ν。选一个均匀分布的随机数RAN，如果就认为该离子发生了电荷交换碰撞。离子发生电荷交换碰撞后，将该位置的中性原子的速度赋给带电离子，产生的交换电荷离子将参与下一时间步的电场计算和离子运动过程。

离子撞击到材料的表面，使材料脱离本体的效应称为溅射效应，它是离子发动机加速栅极下游表面电荷交换碰撞腐蚀形成的物理原因。每个入射离子所轰击出的材料原子的个数称为溅射产额，用Y表示。溅射产额计算公式［14］如下

单个模拟离子的入射能量EXe计算公式为

入射离子的入射角θi对溅射产额的影响是通过上述公式得到的溅射产额进行修正而获得:

3 计算流程

IFE－PIC－MCC三维等离子体离子模拟程序计算离子(包括CEX离子)和中性原子在离子光学系统中的运动。程序可分为3个模块。首先在不考虑电荷交换碰撞的情况下，分别计算离子和中性原子在计算空间的分布。再利用离子和中性原子的空间分布，根据MCC碰撞模块计算CEX离子的产生，最后得到包含CEX离子的所有离子在计算空间的分布。

(1)模拟不考虑电荷交换碰撞的离子在离子光学系统中的运动

在每个时间步长内都有一定数目的模拟粒子从计算域的左边界(z=0平面)进入计算域，在此过程中，程序需要给定的参数包括，每个时间步长进入计算域的模拟粒子的数量，初始位置，初始速度等。进入计算域后，通过给定的边界条件和粒子所带的电量，用体积权重法将离子所带电量向所在网格的节点进行插值，得到网格节点上的电荷密度;利用网格节点上的电荷密度计算出网格节点上的电势和电场;由网格节点上的电场向单个粒子的位置进行插值，计算出该离子位置的电场强度进而计算出离子所受的力;用牛顿第二定律计算出粒子在下一时刻的位置。在下一个步长中，重复以上步骤，直到计算域中离子的数目达到稳定状态。

(2)模拟中性原子在离子光学系统中的运动

中性原子不带电，因此中性原子不受电场的影响，在离子光学系统中匀速运动。在一个时间步长内，通过给定的边界条件，用体积权重法将粒子向所在网格的节点进行插值，计算出网格节点上的粒子密度。计算下一个步长中原子的位置，重复以上步骤，直到计算域中原子的数目达到稳定状态。最后统计出中性原子的粒子密度。

(3)CEX在离子光学系统中的产生

在IFE－PIC中加入蒙特卡洛碰撞模块模拟CEX离子的产生和运动。束流中部分经过电场加速的高速离子和以热速度运动的中性原子之间会发生电荷交换碰撞，产生CEX离子。生成的CEX离子将参与下一个时间步长的离子运动。

4 模拟结果与分析

栅极孔中心线沿z轴方向的电势如图2，从图中可以看出，屏栅极和加速栅极之间存在很高的电势差，可以使带电离子以很高的速度喷出。加速栅极的电势要低于下游边界的电势。图3为栅极孔中心线在z轴方向速度变化。离子在计算域内的初速度是以波尔兹曼分布给出的。经过屏栅极与加速栅极之间的强电场加速后，离子的速度急剧增加，最终以很高的速度喷出，最高速度能到几万米/秒。结合图2的电势分布可以看出，在屏栅极和加速栅极之间离子电势急速下降，离子急剧加速，加速栅极下游区，电势是上升的，离子在该区域稍有减速。

中性原子的粒子密度见图4。中性原子不带电，粒子在计算域中的运动不受电场的影响，做匀速运动。

图2 栅极孔中心线在z轴方向的电势变化(电势以5 V无量纲)

图3 栅极孔中心线上离子速度Z方向分量的变化(速度以1 917.4 m/s无量纲)

图4 中性原子的粒子数密度(以1×1017m－3无量纲)

图5 电荷交换离子的离子数密度(以1×1017m－3无量纲)

图6 2种情况的离子密度比较

图5 是初始离子数密度(即放电室等离子体密度)n0=1.5×1017m－3时的电荷交换离子的数密度，电荷交换离子的数密度比离子数密度小得多。图6是2种情况对应工况下的离子数密度在计算域内的空间分布。上半部分是不考虑电荷交换离子时，带电离子数密度的空间分布;下半部分为考虑电荷交换的离子数密度的空间分布。从图中可以看出，由于栅极系统的聚焦效果，从屏栅极开始到加速栅极，主束流的半径逐渐减小。在正常情况下，主束流半径在加速栅极孔附近小于加速栅极孔半径。由于栅极系统这种聚焦作用，使得主束流边缘的离子数密度比轴线附近的离子数密度小得多。

在加速栅极上游区，2种情况下的离子束密度分布基本没有差别。在加速栅极下游区，不考虑电荷交换情况下，离子在电场作用下聚焦通过加速栅极孔，没有撞击加速栅极。考虑电荷交换情况下，有离子撞击加速栅极。因此与加速栅极碰撞的不是主束流离子而是电荷交换离子;离子与加速栅极碰撞产生的加速栅极电流，主要是由电荷交换离子引起的。2种情况下的栅极电流如表2所列。考虑电荷交换后加速栅极上的栅极电流由0增大到10.68。

表2 栅极电流对比(以8.489×10－8A无量纲)

图7为加速栅极下游表面电流密度分布数值计算结果。电流密度参数体现的是与加速栅极碰撞的离子在栅极表面的分布情况。在孔中心连线形成的等边三角形中心附近的电流密度最大。图8为离子入射角度分布，单位为弧度。图9为加速栅极下游表面溅射产额分布的数值计算结果。图10为加速栅极下游表面的腐蚀深度计算结果。腐蚀深度的分布与电流密度的分布相似，与溅射产额分布趋势相同。图10中的腐蚀深度为发动机工作时间为1 000 h时的腐蚀深度。

图7 加速栅极表面电流分布(A/m2)

图8 离子的入射角分布(rad)

图9 加速栅极表面溅射产额分布

图10 加速栅极下游表面腐蚀深度(m)

图11 为加速栅极孔壁电流密度分布数值计算结果。电流密度参数体现的是离子对栅极孔壁的冲击情况。图中显示靠近加速栅极下游处的离子数密度较大，撞到孔壁的离子数也会较多，在靠近下游的孔壁上产生的电流密度就会增大。图12为离子入射角度分布。离子的入射角度不同，会影响离子的溅射产额。图13为加速栅极孔壁溅射产额分布的数值计算结果，溅射产额与入射能量的分布相近，图中显示溅射产额在孔壁下游部分会比较大。图14为加速栅极孔壁的腐蚀深度计算结果。腐蚀深度的分布与电流密度的分布相似，与溅射产额分布趋势相同。图14中的腐蚀深度为离子发动机工作时间为1 000 h时的腐蚀深度，最大腐蚀深度靠近加速栅极下游，在孔壁可能会形成凹槽。

图11 加速栅极孔壁电流分布(A/m2)

图12 加速栅极孔壁离子的入射角分布(rad)

图13 加速栅极孔壁溅射产额分布

图14 加速栅极孔壁腐蚀深度(m)

5 结论

采用IFE－PIC－MCC三维等离子体离子模拟程序模拟了考虑电荷交换的带电离子在离子光学系统中的运动。电荷交换离子可以使加速栅极冲击电流增加，电荷交换离子是产生栅极电流和栅极腐蚀的主要因素。由于电荷交换离子对加速栅极的冲击作用，加速栅极下游表面会形成“坑”形和“凹槽”形腐蚀，孔壁的腐蚀会造成孔径增大。

［1］SUTTON R G，BIBLARZ O.Rocket Propulsion Elements［M］.John Wiley ＆ Sons，2001.

［2］吴汉基，蒋远大，张志远.电推进技术的应用与发展趋势［J］.推进技术，2003，24(5):385～392.

［3］KAFAFY R.Immersed Finite Element Particle－In－Cell Simulations of Ion Propulsion［D］.Virginia Polytechnic Institute and State University，2005.

［4］PENG X，RUYTEN W，KEEFER D.Three dimensional particle simulation of grid erosion in ion thrusters［R］.IEPC 91－119，1991.

［5］ARAKAWA Y，NAKANO M.An efficient three dimensional optics code for ion thruster research［R］.AIAA 96－3198，1996.

［6］BOND R，LATHAM P.Ion thruster action extraction grid design and erosion modeling using computer simulation［R］.AIAA 95－2923，1995.

［7］HAYAKAWA Y.Three－dimensional numerical model of ion optics system［J］.Journal of Propulsion and Power，1992，8(1):110～117.

［8］GROTE D.Three dimensional simulations of space charge dominated heavy ion beams with applications to inertial fusion energy［D］.University of California，Davis，1994.

［9］SHIRAISH T，KUNINAKA H，SATORI S，et al.Numerical simulation of grid erosion for ion thruster［R］.IEPC 95－90，1995.

［10］WANG J，CAO Y，Kafafy R.Numerical and experimental investigations of crossover ion impingement for subscale ion optics［J］.Journal of Propulsion and Power，1994，24(3):562 ～570.

［11］钟凌伟，刘宇，王海兴，等.电荷交换离子对栅极系统束流影响的数值研究［J］.航空动力学报，2009，24(8):53～56.

［12］熊家贵，王德武.离子引出的二维PIC－MCC模拟［J］.物理学报，2000，49(12):2420～2426.

［13］BIRDSALL C.K，LANGDON A B.Plasma Physics via Computer Simulation［M］.Mcgraw－Hill，1985.

［14］ROSENBERG D，WEHNER G K.Sputtering Yield for Low Energy He+，Kr+，and Xe+Ion Bombardment［J］.Journal of Applied Physics，1962，33(5):1842 ～1845.

［15］李娟，顾左，江豪成，等.氙离子火箭发动机补偿栅极设计［J］.真空与低温，2005，11(1):29～33.