轨道平行栅推力器工作机理的数值研究

汪 忠，李艳武，车清论

（兰州空间技术物理研究所 真空技术与物理重点实验室，兰州 730000）

0 引言

电推进技术是目前空间推进的先进技术，已成为航天器性能先进性的标志[1-2]。电推进的原理为空间推进技术的发展提供了进一步的空间。对近地轨道有一种可能性，就是利用地球高层大气（电离层中F层）中自然存在的稀薄等离子体作为空间发动机的推进剂。太阳辐射的能量致使高层大气部分电离，形成稀薄的等离子体。原理上如果航天器上的适当装置可以直接加速这些轨道环境中的等离子体，则航天器将受到一个反向推力。有关这个设想的基本方案是在航天器上添加一个加有电压的平行双栅结构，当航天器在电离层飞行时，电离层中的等离子体以6 000～8 000 m/s的相对速度进入航天器上的平行栅，其中正离子被双栅之间的电场加速，电子则被栅间电场排斥，由于电子质量远小于正离子，所以产生的阻力可以忽略。2006年美国Northrop Grumman公司提出了利用平行栅加速轨道环境离子的电推进技术方案[3]（The Ambient At⁃mosphere Ion Thruster，AAIT）。此后，美国乔治亚州理工学院研究小组报告了研究结果，建立了利用轨道环境离子的平行栅推力器的数学模型，并计算了推力、阻力等[4-5]。平行栅电场加速轨道离子能否产生推力是研究平行栅推力器要解决的一个核心问题，用二维PIC方法对轨道离子被双栅引出的过程进行数值模拟[6]。

1 计算模型

加有电压的平行双栅的结构类似于离子推力器的屏栅和加速栅，双栅上打有规则分布的同轴小孔，因栅极上的每个孔都具有轴对称性，采用二维轴对称模型[7-8]对单孔进行模拟，如图1所示。

图1 平行双栅电场加速电离层等离子体图

计算域r、z分别为径向和轴向位置，R0、Z0分别为计算域径向和轴向长度，rf为前栅孔半径，rb为后栅孔半径，Φ为电势，下边界为栅孔轴线，坐标原点选在计算域左下角，计算区域的大小为7.5 mm×1.1 mm，划分为300×44的正交等距网格，网格大小为0.025 mm，满足△z≤λD，△r≤λD的要求，其中λD为等离子体的德拜长度。计算域的边界条件分为电场边界条件和粒子运动边界条件两个部分，如图2所示。计算域左右边界的电场边界条件取为Φ=0，前栅电势取为0，后栅电势取为-1 000 V，上下边界的电场边界条件取为Neumann型边界条件。粒子运动边界条件分为两种。上下边界为反射边界，左右边界和栅极为吸收边界。模拟过程中，如果粒子超出左右边界或者撞到栅极，则将粒子删除，如果达到上下边界，则将粒子镜面反射回计算域，即粒子z向速度不变，r向速度相反。

双栅产生的静电势和运动等离子体产生的自洽电势可通过求解Possion方程来得到：

式中：e为电子电量；ε0为真空介电常数；ni为离子数密度；ne为电子数密度。ni通过以网格点为中心的一个网格大小的方形区域的离子个数除以网格面积求得。将电子视为热平衡流体，ne通过Boltzmann方程求得。采用有限差分法求解出电势后，可求出电场强度：

图2 计算域及边界条件图

PIC方法的基本思想为：通过求解电场Poisson方程得到电势分布，用插值法求得粒子所在位置处的电场力，然后加速粒子，再用新的粒子位置求解新的电场分布，上述过程循环进行，直至收敛。

粒子以初始位置和初始速度从计算域左边界进入，经电场作用后再从边界离开。在粒子初始化中，须确定粒子的初始位置与初始速度，粒子初始位置都位于计算域左边界发射面上，其在发射面上位置随机给定，并保证粒子在整个圆形发射面上均匀分布，粒子初始发射速度z向分量由航天器在轨速度和Maxwell速度确定：

发射速度径向分量满足Maxwell速度分布：

式中：vM为Maxwell速度；θ是在0到π间均匀分布的随机数。

粒子加速方法采用蛙跳格式，粒子运动差分方程为：

式中：m为粒子质量；v为粒子速度；F为电场力；q为电量；x为粒子位置，电场强度由与粒子邻近的4个网格点上的电场强度插值得到。

时间步长的选取必须满足：

式中：ωp为等离子体频率；选取tp=1.0×10-9s。

为简化问题，在粒子模拟过程中，只考虑轨道等离子体环境中的一价单原子氧离子、电子，不考虑中性原子。等离子体中存在多种粒子之间的碰撞，如带电粒子之间的库仑碰撞，离子与中性粒子间的电荷交换碰撞等，但在低轨道环境中，粒子密度较低，分子自由程较大，因此，不考虑粒子间的碰撞。

2 结果分析

模拟程序需要输入的参数包括：栅极的几何参数、栅极电压、等离子体参数。所取参数如表1、2所列，其中等离子体参数以300 km高度轨道的电离层等离子体为例[9]。

表1 栅极参数

表2 等离子体参数

在300 km高度，主要的离子成分为单原子氧离子，在模拟过程中，只考虑O+而忽略其他离子成分（如He+）。平行栅上游离子密度niu是正常电离层300 km高度的数值，下游离子密度根据上游粒子密度计算得到，具体如式（7）：

式中：vu为上游离子速度，即离子与航天器的相对速度为7 700 m/s；vd为下游离子速度；ΔV为双栅电势差；m为氧离子质量。

从图3、4中可以看出，栅极间的电势线几乎垂直于z轴，因此在z轴方向上形成极强的电场将进入的离子加速引出，而当离子被引出后栅，电势则开始上升，离子受到与z轴相反方向的电场开始减速，这是因为加有负高压的后栅在等离子体中将形成较大尺寸的鞘层（几分米至几米），鞘层处的电位与空间电位相当，可认为是0电位，后栅和鞘层之间形成与z相反方向的电场。值得说明的是在模拟中，计算域z轴方向的尺度应与鞘层尺度相当，但这将耗费相当长的机时。因此，模拟过程未将计算域z轴取到鞘层尺度，但z轴方向尺度已经取得较大，可认为右边界处即为鞘层，这样的做法对模拟结果的影响可以忽略。

图3 电势分布图

图4 粒子轨迹图

如图5所示，左边界电位为0，在平行栅中，电位下降△U，到平行栅下游，电位开始上升，到右边界，电位恢复到0，离子在这一过程中电场势能没有改变，因此并没有从平行栅电场中获得能量，即离子在这一过程中没有被加速。

图5 轴线处电势随z轴变化曲线图

King等[5]提出在平行栅前增加一个离子势能提升区，用来提升离子势能，这样，当离子到达平行栅下游的鞘层处，多余的势能转化为动能，即离子最终被加速喷出，航天器获得反向推力。

图6和图7是在离子势能提升区将电势提升到50 V后的模拟结果。比较图4和图7可以看出，加了离子势能提升区后，离子在加速栅中的聚焦变的明显。区的办法，平行栅则可以加速轨道中的离子以获得反向推力。

图6 加离子势能提升区的电势分布图

图7 加离子势能提升区的粒子轨迹图

如图8所示，由于离子势能提升，到下游鞘层处，电势恢复到0时离子还剩余部分动能，即离子被加速喷出，航天器将获得反向推力。

图8 加离子势能提升区后轴线处电势随z轴变化曲线图

3 结论

模拟了不加离子势能提升区和加上离子势能提升区两种情况下离子被平行栅电场引出的情况。模拟结果表明：加电压的平行栅无法直接加速轨道离子，利用在平行栅前增加一个离子势能提升

[1]张天平，张雪儿.空间电推进技术及应用新进展[J].真空与低温，2013，19（4）：187-193.

[2]张天平，周昊澄，孙小菁，等.小卫星领域应用电推进技术的评述[J].真空与低温，2014，20（4）：187-192.

[3]Dressler G A.Spacecraft propulsive device using ambient up⁃per atmospheric constituents for reaction mass[C]//42nd Joint Propulsion Conference&Exhibit，2006.

[4]King S T，Walker M L R，Chianese S G.Ambient Atmosphere Ion Thruster（AAIT）Proof-of-Concept Modeling[C]//47thJoint Propulsion Conference&Exhibit，2011.

[5]King S T，Walker M L R，Silvio G.Chianese.Atmospheric Elec⁃tric Propulsion Mission Performance Tool[J].Journal of Space⁃craft and Rockets，2014，51（3）：931-937.

[6]邵福球.等离子体离子模拟[M].北京：科学出版社，2002：9-11.

[7]刘畅，汤海滨，顾佐，等.基于PIC方法的离子发动机光学系统粒子模拟[J].北京航空航天大学学报，2006，32（4）：382-386.

[8]刘洋，李娟，楚豫川，等.考虑电荷交换的栅极区离子流数值模拟[J].真空与低温，2011，17（4）：198-208.

[9]庄洪春.宇航空间环境手册[M].北京：中国科学技术出版社，2000：116-117.