离子推力器羽流热效应仿真分析

2018-10-30张建华李晶华尤凤仪郑鸿儒

北京航空航天大学学报 2018年10期

张建华，李晶华，尤凤仪，郑鸿儒

(北京航空航天大学宇航学院, 北京 100083)

电推力器以其比冲高、寿命长和可控性强等优点在深空探测中发挥着重要作用，已经被应用于姿态控制、轨道转移、深空探测等空间任务中[1-2]。电推进技术的在轨应用为航天器空间环境研究带来了许多无法忽视的新难题，电推进羽流与航天器的相互作用是其中重要的内容之一。

离子推力器羽流主要由离子推力器高速喷出的高能工质离子与中和器发射的电子构成的束流等离子体及与电荷交换产生的交换电荷等离子体两部分组成。羽流会撞击航天器表面，对航天器产生力效应、热效应和污染效应等，影响航天器部件的正常工作和寿命，严重时会导致航天任务失败[3]。其中，热效应主要是由于离子推力器羽流会直接或者返流与航天器表面发生碰撞，这种碰撞会引起航天器敏感材料热变形，而热流的累积也会造成航天器重要元器件的损坏，进而导致航天器重要部件功能的丧失。因此，开展离子推力器的羽流热效应研究十分必要。

目前，国内外电推进的诊断测量主要利用静电探针、发射探针、阻滞势分析仪(RPA)、法拉第探针等接触式诊断设备直接放置在等离子羽流中，对电子密度、空间电位、离子能量分布、等离子体电势、电流密度等加以测量，然而关于等离子体的热效应研究正式发表的文献很少。2007年，上海交通大学周志雄等[4]对真空条件下气体与固体表面的相互作用过程的热适应系数进行了数值模拟仿真研究。2008年，美国奥斯汀大学Berisford等[5]以螺旋波等离子体为研究对象，对氩放电过程进行了热效应诊断实验研究，用热电偶和热辐射探针分别安装在气流的下游和真空舱尾部测量热流值。

对于推力器热效应的仿真研究主要集中于化学推力器。2014年，王黎珍等[6]对化学推力器真空羽流热效应的计算模型进行了修正和误差分析。目前，尚未见针对电推进羽流热效应的仿真研究，而仿真研究对于电推力器的实际在轨应用具有重要意义。如果分析过于保守，则导致热防护过度；估计不足，则会使热量积聚，可能导致设备失效。在电推进羽流热效应分析中，如何准确获得羽流流场及高能粒子与表面作用的精确模型是2个主要难点。本文针对LIPS-200型离子推力器地面真空舱内羽流热效应实验中的部分工况，采用三维粒子网格-直接模拟蒙特卡罗(PIC-DSMC)方法进行数值模拟研究，并与实验数据进行对比分析。

1 仿真方法

本文涉及到的仿真算例均采用EX-PWS(Extension of Plume Work Station)[7]软件实现。EX-PWS是羽流工作站(PWS)[8]软件的电推进扩展程序包，使用PIC方法[9]处理电场解算、粒子运动等过程，使用DSMC方法[10]处理粒子间碰撞。EX-PWS软件使用三维非结构网格适应复杂空间环境。此外，为提高计算效率，软件使用了MPI并行计算方法。

为了减少计算量，EX-PWS将等离子体中的电子视为流体，将离子及推力器出口喷出的Xe原子视为计算粒子。电子运动采用准中性假设[11]和Boltzmann关系式描述。电子温度使用等温模型，忽略磁场影响。Boltzmann关系式[12]可表述为

(1)

式中：ne为电子密度；nref为电势φ为零处的参考电子密度，这里取1×1016m-3；Te为电子温度，这里取为常数3.5 eV；k为Boltzmann常数。

粒子间碰撞包括弹性碰撞及电荷交换碰撞，其是影响流场的一个主要因素。在软件中，粒子间碰撞使用DSMC方法进行处理。忽略Xe原子、Xe离子与碳原子之间的弹性碰撞和电荷交换碰撞。各种粒子间的碰撞截面选取与文献[13]相同。

本文在模拟环境压强时采用虚拟粒子的方法，将背景气体视为虚拟粒子，仅参与碰撞计算。当计算粒子运动到某一网格中将要与背景Xe原子发生碰撞时，按照设置的环境压强给出数密度，按照环境温度根据Maxwellian 分布给出热运动速度，碰撞后背景粒子状态不变，计算粒子运动速度及方向，按照碰撞模型发生变化。计算粒子与背景粒子间的碰撞类型包括弹性碰撞及电荷交换碰撞。

2 实验系统

实验系统由真空羽流实验设备、离子推力器系统、热流传感器系统和Faraday探针系统组成，本文主要介绍测量热流的传感器系统。热流传感器系统主要由总热流计、辐射热流计、信号转换装置和电流信号采集系统组成。Schmidt-Boelter[14-16]热流计是基于空间温度梯度原理，利用热电堆测量温度差的方法，进而获得热流的大小。实验测量系统的示意图如图 1所示，本次实验过程中使用的热流传感器实物图如图 2所示。

本文进行离子推力器羽流热效应诊断采用的热流计为Schmidt-Boelter热流计。热流计的量程为10 kW/m2，精度满量程±3%，因此热流计的诊断误差在±0.3 kW/m2左右。同时需要注意的是，热流计最高容许传感器温度为204℃，在具体的实验过程中采用两路热电偶对热流传感器的温度进行监控，并采用降温处理使热流传感器的温度保持在合理的温度范围内。

图1 离子推力器羽流热效应测量系统示意图Fig.1 Schematic diagram of ion thruster plume thermal effect measurement system

图2 两种热流传感器实物图Fig.2 Photo of two kinds of heat flow sensor

3 网格及边界条件

仿真计算区域设定为包括推力器前方1.5 m和推力器后方0.5 m、直径2.0 m的卧式圆柱体，网格数为200万左右，沿x方向划分为24个并行分区。如图3所示，计算域内部有一个小圆柱体用来模拟推力器，模拟热流计为直径30 mm、长40 mm的小圆柱体，模拟热流计与推力器间的相对位置如图4所示。

在处理粒子与边界面之间的相互作用时，采用适应系数衡量来流粒子速度对壁面温度的适应程度。在粒子与表面发生碰撞时，程序会生成一个随机数，其大于适应系数，则粒子采用镜面反射，否则，采用漫反射。文献[17]中提出当离子能量高于400 eV时，适应系数将接近于1.0，为保守起见，本文中热流计表面适应系数取为0.9，即认为90%的粒子进行漫反射，10%的粒子进行镜面反射。计算域边界设置为自由边界，当粒子运动到边界后即从流场中移除该粒子。

图3 计算域示意图Fig.3 Schematic diagram of computational domain

图4 推力器与模拟热流计相对位置示意图Fig.4 Schematic diagram of relative position of thruster and simulated heat flow meter

实验中采用的推力器为兰州空间技术物理研究所研制的1.3 kW的LIPS-200型离子推力器[18]。推力器工质为氙气，工作时的流量为1.37 mg/s,栅极电压为1 000 V。仿真中设置推力器电离率为90%，羽流发散半角为16°，二价离子占离子总数的10%。推力器仿真参数设置如表 1所示。推力器出口射出的粒子包括Xe原子、一价离子和二价离子3种，出口流量分布设置为高斯分布[19]。计算中忽略了中和器带来的影响。

为了出口分布的准确性，针对离子推力器在实验中的工作情况进行了仿真，并与实验结果进行对比。仿真结果如图5所示，实验数据由Faraday探针测得，实验在真空羽流实验设备(PES)[7，20]中开展。可以看出，仿真结果与实验结果符合较好。因为近场羽流中的电流密度大部分来自于推力器直接喷出的高速离子，因此近场电流密度仿真结果与实验结果相符能够说明仿真所用的推力器出口参数是合理的。

表1 LIPS-200型离子推力器仿真基本参数

图5 羽流场电流密度仿真与实验对比Fig.5 Comparison of current density in plume flow field between simulation and experiment

4 计算结果分析

本节针对LIPS-200型离子推力器，开展电推进羽流热效应仿真分析，并与实验结果进行对比。本文对羽流热效应的计算综合考虑了热流计置于流场中对羽流的影响，尤其是高能粒子撞击热流计表面后被中和并损失能量成为慢速中性原子后，仍然参与Xe原子流场的计算，使得仿真条件与实验工况尽量吻合。计算中，时间步长设置为1×10-7s，粒子权重设置为5×108，环境及各壁面温度为300 K，环境压强为2 mPa，并认为推力器及热流计接地，所有表面电势设置为0 V。

图6给出了轴向位置上(距离推力器出口0.5～1.0 m)实验和仿真的对比。其中，平均滞止热流仿真值曲线显示粒子与热流计表面进行了完全的热交换(热适应系数为1.0)；平均热流仿真值曲线即带有热适应系数的仿真热流值(本文中的适应系数取值为0.9)。需要注意的是，仿真值中的平均指的是在30 mm的模拟热流计表面进行能量统计后在表面上进行平均，避免在单个网格中由于粒子数较少造成热流值畸高。

可以看出，在轴线上仿真值与实验测量值均沿着与推力器的距离升高而降低，实验测量值与平均滞止热流仿真值相近，最大误差出现在0.9 m处，与平均滞止热流仿真值的误差值为17.0%，与平均热流仿真值的误差为25.3%，在等离子体测量领域认为是可以接受的。但同时可以看出，本文算例中采用的适应系数偏低，当粒子能量全部传递(适应系数为1.0)，即粒子滞止时热流仿真值与实验测量值符合较好。仿真与实验误差产生的可能原因包括束流模型偏差、适应系数选取及实验误差等方面。

图7(a)给出了距离推力器出口0.5 m处径向热流随角度的分布情况。可以看出，在角度为0°～15°的羽流束流区域内，实验测量值基本小于平均滞止热流仿真值，而高于平均热流仿真值。平均滞止热流相对误差绝对值在11.1%以内，平均热流相对误差绝对值在21.1%以内，误差较小。平均滞止热流相对误差较大的位置在15°左右，此处为束流半角，即束流区和非束流区的交界处，此处的等离子体数密度较为稀薄，仿真和实验的难度都很大，因此相对误差较高。但同时也能看出，无论是平均滞止热流仿真值还是平均热流仿真值，其沿径向下降的速度与实验测量值相比较为平缓，说明在采用Maxwell模型处理粒子与表面能量交换时使用的适应系数不应是固定值，而是应该随着角度进行变化。从粒子撞击的角度来说，正入射时粒子垂直撞击热流计表面，绝大部分能量转化为热能，而斜入射粒子存在一定的出射速度，从而降低了能量交换。因此，当角度升高时，适应系数应随之降低。

图6 轴向位置热效应实验和仿真对比Fig.6 Comparison of thermal effect in axial direction between simulation and experiment

图7 距离推力器出口0.5 m和0.9 m处径向热流对比Fig.7 Comparison of radial heat flow at 0.5 m and 0.9 m from thruster exit

图7(b)给出了距离推力器出口0.9 m处径向热流随角度的分布。可以看出，在羽流束流区域内，实验测量值全面高于平均滞止热流仿真值和平均热流仿真值。平均滞止热流相对误差绝对值在18.4%以内，平均热流相对误差绝对值在26.5%以内，误差在可接受范围。但同时能够看出，在0.9 m处径向小角度内误差普遍较高，说明仿真中采用的束流模型在0.5 m后下降速度较快，与实际分布相差较大。

图8(a)给出了热流计处于距离推力器出口0.9 m处流场中的一价Xe离子分布情况。可以看出，对于一价Xe离子来说，热流计的影响范围仅在热流计后约0.1 m的范围内，热流计前方流场基本不受影响，且对流场整体影响较小。

图8 距离推力器出口0.9 m处放置热流计时的一价Xe离子和Xe原子分布Fig.8 Number density distribution of Xe atom and monovalent xenon ion when heat flow meter is located at 0.9 m away from thruster exit

图 8(b)给出了有热流计情况下流场中的Xe原子分布情况。可以看出，模拟热流计对Xe原子的影响同时存在于热流计的前方及后方部分区域。在热流计前方，滞止于热流计表面的Xe离子被中和成为Xe原子，使小范围内的粒子数密度升高，数密度级别与推力器出口相近。在热流计后方同时出现数密度较低的区域，但由于Xe原子运动速度较低，且粒子数密度相比Xe离子较高，易于实现空间积聚，因此并未出现粒子数密度为0的区域。

图9给出了距离推力器出口0.9 m处放置热流计后轴线上的一价Xe离子数密度、Xe原子数密度分布与无热流计时数密度的对比。可以看出，对于一价Xe离子来说，影响较大的区域在热流计测量平面后的0.1 m范围内，在热流计后方距离较近区域接近为0，然后逐渐恢复至正常量级，但在计算域内仍低于无热流计的情况。这是由于一价Xe离子的速度较高，约为3.8×104m/s，在无干扰的情况下基本沿直线运动，很难受到影响，当高速运动的Xe离子撞击热流计表面滞止后，被中和成为Xe原子，同时能量传递给热流计，因此造成了热流计后方小范围内粒子数为0。而对于Xe原子，由于一部分Xe离子被中和后返回Xe原子流场中，造成了热流计前方0.1 m范围内Xe原子数密度升高约2个量级，但由于在稀薄气体中粒子间的碰撞概率较低，小范围内Xe原子数密度的升高对流场影响不大。

图9 距离推力器出口0.9 m处放置热流计时的一价Xe离子、Xe原子及无热流计时数密度对比Fig.9 Comparison of xenon atom and monovalent xenon ion number density distribution on axis with or without heat folw meter located at 0.9 m away from thruster exit