基于环境工质的超低轨吸气式螺旋波电推进器仿真分析

丁 亮，彭毓川，郑慧奇，唐振宇，任琼英，赵 华，秦 玮

（北京卫星环境工程研究所，北京 100094）

0 引言

大气阻力是制约航天器在超低轨道工作寿命的主要原因之一。而如果航天器能够充分利用轨道环境中的残余大气作为工质，经过收集、电离、加速后为航天器提供空间推进动力，维持航天器的轨道速度，将大大减少航天器所需携带的推进剂量，延长航天器的在轨运行寿命。Nishiyama[1]提出利用电子回旋共振（ECR）技术的吸气式离子发动机概念，以期维持航天器的超低轨长时间运行。Pekker[2]研究霍尔电推进技术，以超低轨道上残余气体为工质，理论上可以将航天器维持在80～90 km高度的轨道上。这2种方案都是先将轨道残余大气经过压缩收集到高压容器中，然后向推进器供气。其主要难点一是如何高效收集轨道残余大气并实现上万倍的密度压缩；二是轨道残余大气的主要成分为氧原子，而在离子电推进和霍尔电推进这些有电极加速机制的推力器上，氧原子对电极的剥蚀作用将导致电推进器的寿命大幅缩短，难以实现长时间的轨道维持[3-5]。

螺旋波电推进采用射频功率耦合，无电极腐蚀过程，从原理上能够避免电推进器的上述寿命缺陷。国内外针对螺旋波电推进技术已开展了大量研究[6-12]：Chen最先开展了螺旋波等离子体的研究，分析了螺旋波电离、放电模式、功率吸收与模式激发过程；Charles等人将螺旋波等离子体应用到了电推进领域，并分析了推进模式和加速机制。但目前还未见关于螺旋波等离子体在吸气式电推进技术中应用的相关报道。

为了验证超低轨道高度环境下吸气式螺旋波电推进的可行性，本文进行以原子氧为工质，综合考虑碰撞、激发、电离等过程，通过13.56 MHz射频加热，由磁喷口完成推力输出全流程的仿真模拟；建立特定形态吸气式螺旋波电推进结构模型，施加不同功率的输入，分析功率沉积、等离子体参数分布和推力输出。

1 模拟分析的原理

在超低轨道环境中，以180 km高度为例，气体主要成分为原子氧（O）。收集的轨道残余气体通过射频波加热后呈等离子体态，故可采用等离子体物理模型模拟该气体流。在等离子体模型中考虑电子（e）、氧离子（O+）、激发态氧原子（Os）以及本底氧原子（O）4种粒子的电化学反应：

对于弹性碰撞，有

对于激发过程，有

对于退激过程，有

对于电离过程，有

同时，取吸附系数为1时壁面电化学反应为：

在所采用的等离子体模型中，电子密度ne、平均电子能量nε的漂移‒扩散方程为：

式(7)～(9)中：μe为电子迁移率；E为等离子体电场强度；De为电子扩散系数；Re为电子源项；με为电子能量迁移率；Dε为电子能量扩散系数；Γe为电子通量；Rε为非弹性碰撞损失项；Qrh为波加热项贡献；J为等离子电流；E*为射频波场。其中，De、Dε、με的计算式为

假设有M种反应净生成电子，则电子源项为

假设有P种非弹性碰撞反应导致能量变化，则非弹性碰撞损失项为

式(10)～(12)中：Te为电子温度；xj为碰撞目标粒子的份额；kj为类反应的速度系数；Nn为总的中性粒子密度；Δεj为j类型反应转移的能量。

对于氧离子，可通过等离子体电化学反应计算获得其密度分布。对于非带电的激发态原子氧，采用重粒子的对流扩散方程进行模拟计算，

其中：ρ为电荷密度；wk为总粒子数的份额；u为等离子体宏观流体速度矢量；jk为激发态原子氧扩散通量；Rk为激发态原子氧的生成速率源项，由电化学反应决定。假设有n种带电重粒子，则等离子体电场的计算方程为

其中：ε0为真空介电常数；εr为相对介电常数；V为电势；电荷密度ρ由等离子体化学反应过程耦合计算获得；q为单位电荷；Zk为带电粒子电荷数；nk为带电重粒子密度。

采用Maxwell方程组模拟加热等离子体的电磁场：

式(15)～(18)中：B为磁感应强度；J为等离子体电流；A为磁矢势；D为电位移矢量；假设等离子体各向同性，电导率σ为

其中：me为电子静止质量；γe为电子‒中性粒子碰撞频率。

采用可膨胀层流模型方程组模拟原子氧流体变化：

式(20)～(22)中：U为原子氧流体的速度矢量；P为流体压力；I为单位张量；μ为等离子体黏度系数，通过各粒子的份额综合计算获得；Cp为原子氧流体热容；Q为热源项，由等离子体电化学反应方程耦合获得；ρ是由各粒子的份额综合计算获得的密度值。

整体上将洛伦兹、Maxwell方程组耦合的热流磁流体方程对装置吸入气体的热场、流速场进行模拟，并通过流体的物性与前述的等离子体模型中的粒子密度、温度参数进行耦合，可形成一个评估波加热等离子体效应对宏观流体力学性能影响的物理模型，进而得出装置的推力效果。该模型通过有限元COMSOL-MULTIPHYSIC软件平台进行求解。模型网格划分为：等离子体壁面边界采用四边形网格，线圈区域采用四边形网格，其他区域采用三角形网格；总网格数 750 000，时间步长 1×10-5s。

当达到热平衡时，推进器产生的推力Fthrust由以下方程计算：

质量流量为

沉积功率为

式 (23)～(25)中：Vo、ρo、So分别为出口的流体速度、密度、截面积；ΔW为沉积在原子氧流体中的功率。

2 几何模型及计算条件

吸气式螺旋波电推进器的几何模型如图1所示，其装置尺寸见表1。

图1 吸气式螺旋波电推进器仿真结构Fig. 1 Simulated structure of the air-breathing helicon plasma propulsion

表1 吸气式螺旋波电推进器仿真结构装置尺寸Table 1 Sizes of simulated structure of the air-breathing helicon plasma propulsion

推进器运行轨道高度设计为180 km。在该轨道环境下的残余气体粒子密度为15.4×1015m-3，飞行器的飞行速度为7800 m/s，则气体收集口的粒子通量为 1.2×1020m-2·s-1。设计装置的入口截面积为0.020 1 m2，相应的，单位时间吸入的粒子数为 2.4×1018s-1，约为 4.0×10-6mol/s，对应为 178 sccm。根据粒子密度推算大气压力约为1.9×10-5Pa，环境温度约为650 K。用以约束粒子的偏转磁场，最强处设计为1500 Gs，磁场位型如图2所示。轴向0点设定为收集口喉部位置，为在直观反映计算结果的同时减小计算量，计算过程示以半剖面。

图2 约束磁场位型Fig. 2 The constrained magnetic field

3 不同驱动功率下的仿真计算结果

200 W的输入驱动功率水平下，沉积在原子氧流体中的功率为2.478 W，出口处质量流量为7.511 9×10-6kg/s，计算得到推力值约为 6 mN。相应的，达到稳态时功率沉积、电子密度分布、电子温度分布、流体速度分布如图3所示，电子密度峰值1.24×1018m-3，电子温度峰值3.25 eV。

800 W的输入驱动功率水平下，沉积在原子氧流体中的功率为6.452 7 W，出口处质量流量为7.511 9×10-6kg/s，计算得到推力值为 9.85 mN。相应的，达到稳态时功率沉积、电子密度分布、电子温度分布、流体速度分布如图4所示，电子密度峰值3.52×1018m-3，电子温度峰值 3.45 eV。

图3 200 W 输入功率仿真结果Fig. 3 Simulation result with input power of 200 W

图4 800 W 输入功率仿真结果Fig. 4 Simulation result with input power of 800 W

1200 W的输入驱动功率水平下，沉积在原子氧流体中的功率为8.285 4 W，出口处质量流量为7.511 9×10-6kg/s，计算得到推力值为 11.16 mN。相应的，达到稳态时功率沉积、电子密度分布、电子温度分布、流体速度分布如图5所示，电子密度峰值5.15×1018m-3，电子温度峰值 3.61 eV。

图5 1200 W 输入功率仿真结果Fig. 5 Simulation result with input power of 1200 W

2000 W的输入驱动功率水平下，沉积在原子氧流体中的功率为11.71 W，出口处质量流量为7.511 9×10-6kg/s，计算得到推力值为 13.23 mN。相应的，达到稳态时功率沉积、电子密度分布、电子温度分布、流体速度分布如图6所示，电子密度峰值7.53×1018m-3，电子温度峰值 3.72 eV。

流体沉积功率、推力随驱动功率的变化如图7所示。

图6 2000 W 输入功率仿真结果Fig. 6 Simulation result with input power of 2000 W

图7 流体沉积功率、推力随驱动功率的变化Fig. 7 Deposited power and thrust versus input power

4 结束语

通过以上仿真分析发现，随着驱动功率的增加，推进器的推力和流体电子密度均显著增加——输入功率从200 W增加至2000 W的过程中，推力从6.00 mN增加至13.23 mN。而吸气式螺旋波电推进器的设计应用轨道环境截面阻力约0.226 mN，因此该推进器能够满足阻力补偿的要求。电子密度峰值从 1.24×1018m-3增加至 7.53×1018m-3，增幅较大；电子温度峰值缓慢增加，从3.25 eV增加至3.72 eV。电磁辐照的功率沉积量较低，但沉积功率的密度较高，在 2000 W 时到达了最大值 105W/m3。以上数据表明：随着输入功率的增加，较多功率用于电离，而温度升高得较少，说明能量获得能力较低，从而导致推力涨幅不大，仅增加了1倍多。其中可能的原因是通道直径较小，粒子复合率较高，输入的功率大部分消耗在将复合后的粒子重新电离上。因此，需要考虑加大加热段的容积来增加功率的沉积值，以增大装置的推力。

总之，根据上述推力计算结果，在下一步工作中，将考虑加大入口截面积、增加气体的质量流量、增大沉积功率，继而增大推进器的推力。