离子推力器羽流主要由离子推力器高速喷出的高能工质离子与中和器发射的电子构成的束流等离子体及与电荷交换产生的交换电荷等离子体两部分组成。羽流会撞击航天器表面，对航天器产生力效应、热效应和污染效应等，影响航天器部件的正常工作和寿命，严重时会导致航天任务失败[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×1016 m－3；Te为电子温度，这里取为常数3.5 eV；k为Boltzmann常数。
粒子间碰撞包括弹性碰撞及电荷交换碰撞，其是影响流场的一个主要因素。在软件中，粒子间碰撞使用DSMC方法进行处理。忽略Xe原子、Xe离子与碳原子之间的弹性碰撞和电荷交换碰撞。各种粒子间的碰撞截面选取与文献[13]相同。
本文在模拟环境压强时采用虚拟粒子的方法，将背景气体视为虚拟粒子，仅参与碰撞计算。当计算粒子运动到某一网格中将要与背景Xe原子发生碰撞时，按照设置的环境压强给出数密度，按照环境温度根据Maxwellian分布给出热运动速度，碰撞后背景粒子状态不变，计算粒子运动速度及方向，按照碰撞模型发生变化。计算粒子与背景粒子间的碰撞类型包括弹性碰撞及电荷交换碰撞。
2 实验系统 实验系统由真空羽流实验设备、离子推力器系统、热流传感器系统和Faraday探针系统组成，本文主要介绍测量热流的传感器系统。热流传感器系统主要由总热流计、辐射热流计、信号转换装置和电流信号采集系统组成。Schmidt-Boelter[14-16]热流计是基于空间温度梯度原理，利用热电堆测量温度差的方法，进而获得热流的大小。实验测量系统的示意图如图 1所示，本次实验过程中使用的热流传感器实物图如图 2所示。
 |
| 图 1 离子推力器羽流热效应测量系统示意图 Fig. 1 Schematic diagram of ion thruster plume thermal effect measurement system |
| 图选项 |
 |
| 图 2 两种热流传感器实物图 Fig. 2 Photo of two kinds of heat flow sensor |
| 图选项 |
本文进行离子推力器羽流热效应诊断采用的热流计为Schmidt-Boelter热流计。热流计的量程为10 kW/m2，精度满量程±3%，因此热流计的诊断误差在±0.3 kW/m2左右。同时需要注意的是，热流计最高容许传感器温度为204℃，在具体的实验过程中采用两路热电偶对热流传感器的温度进行监控，并采用降温处理使热流传感器的温度保持在合理的温度范围内。
3 网格及边界条件 仿真计算区域设定为包括推力器前方1.5 m和推力器后方0.5 m、直径2.0 m的卧式圆柱体，网格数为200万左右，沿x方向划分为24个并行分区。如图 3所示，计算域内部有一个小圆柱体用来模拟推力器，模拟热流计为直径30 mm、长40 mm的小圆柱体，模拟热流计与推力器间的相对位置如图 4所示。
 |
| 图 3 计算域示意图 Fig. 3 Schematic diagram of computational domain |
| 图选项 |
 |
| 图 4 推力器与模拟热流计相对位置示意图 Fig. 4 Schematic diagram of relative position of thruster and simulated heat flow meter |
| 图选项 |
在处理粒子与边界面之间的相互作用时，采用适应系数衡量来流粒子速度对壁面温度的适应程度。在粒子与表面发生碰撞时，程序会生成一个随机数，其大于适应系数，则粒子采用镜面反射，否则，采用漫反射。文献[17]中提出当离子能量高于400 eV时，适应系数将接近于1.0，为保守起见，本文中热流计表面适应系数取为0.9，即认为90%的粒子进行漫反射，10%的粒子进行镜面反射。计算域边界设置为自由边界，当粒子运动到边界后即从流场中移除该粒子。
实验中采用的推力器为兰州空间技术物理研究所研制的1.3 kW的LIPS-200型离子推力器[18]。推力器工质为氙气，工作时的流量为1.37 mg/s, 栅极电压为1 000 V。仿真中设置推力器电离率为90%，羽流发散半角为16°，二价离子占离子总数的10%。推力器仿真参数设置如表 1所示。推力器出口射出的粒子包括Xe原子、一价离子和二价离子3种，出口流量分布设置为高斯分布[19]。计算中忽略了中和器带来的影响。
表 1 LIPS-200型离子推力器仿真基本参数 Table 1 Basic parameters of LIPS-200 ion thruster simulation
| 粒子种类 | 流率/s－1 | 温度/K | 速度/(m·s-1) |
| Xe | 5.69×1017 | 300 | 325 |
| Xe+ | 4.609×1018 | 46400 | 39000 |
| Xe++ | 5.12×1017 | 46400 | 55154 |
表选项
为了出口分布的准确性，针对离子推力器在实验中的工作情况进行了仿真，并与实验结果进行对比。仿真结果如图 5所示，实验数据由Faraday探针测得，实验在真空羽流实验设备(PES)[7, 20]中开展。可以看出，仿真结果与实验结果符合较好。因为近场羽流中的电流密度大部分来自于推力器直接喷出的高速离子，因此近场电流密度仿真结果与实验结果相符能够说明仿真所用的推力器出口参数是合理的。
 |
| 图 5 羽流场电流密度仿真与实验对比 Fig. 5 Comparison of current density in plume flow field between simulation and experiment |
| 图选项 |
4 计算结果分析 本节针对LIPS-200型离子推力器，开展电推进羽流热效应仿真分析，并与实验结果进行对比。本文对羽流热效应的计算综合考虑了热流计置于流场中对羽流的影响，尤其是高能粒子撞击热流计表面后被中和并损失能量成为慢速中性原子后，仍然参与Xe原子流场的计算，使得仿真条件与实验工况尽量吻合。计算中，时间步长设置为1×10-7 s，粒子权重设置为5×108，环境及各壁面温度为300 K，环境压强为2 mPa，并认为推力器及热流计接地，所有表面电势设置为0 V。
图 6给出了轴向位置上(距离推力器出口0.5~1.0 m)实验和仿真的对比。其中，平均滞止热流仿真值曲线显示粒子与热流计表面进行了完全的热交换(热适应系数为1.0)；平均热流仿真值曲线即带有热适应系数的仿真热流值(本文中的适应系数取值为0.9)。需要注意的是，仿真值中的平均指的是在30 mm的模拟热流计表面进行能量统计后在表面上进行平均，避免在单个网格中由于粒子数较少造成热流值畸高。
 |
| 图 6 轴向位置热效应实验和仿真对比 Fig. 6 Comparison of thermal effect in axial direction between simulation and experiment |
| 图选项 |
可以看出，在轴线上仿真值与实验测量值均沿着与推力器的距离升高而降低，实验测量值与平均滞止热流仿真值相近，最大误差出现在0.9 m处，与平均滞止热流仿真值的误差值为17.0%，与平均热流仿真值的误差为25.3%，在等离子体测量领域认为是可以接受的。但同时可以看出，本文算例中采用的适应系数偏低，当粒子能量全部传递(适应系数为1.0)，即粒子滞止时热流仿真值与实验测量值符合较好。仿真与实验误差产生的可能原因包括束流模型偏差、适应系数选取及实验误差等方面。
图 7(a)给出了距离推力器出口0.5 m处径向热流随角度的分布情况。可以看出，在角度为0°~15°的羽流束流区域内，实验测量值基本小于平均滞止热流仿真值，而高于平均热流仿真值。平均滞止热流相对误差绝对值在11.1%以内，平均热流相对误差绝对值在21.1%以内，误差较小。平均滞止热流相对误差较大的位置在15°左右，此处为束流半角，即束流区和非束流区的交界处，此处的等离子体数密度较为稀薄，仿真和实验的难度都很大，因此相对误差较高。但同时也能看出，无论是平均滞止热流仿真值还是平均热流仿真值，其沿径向下降的速度与实验测量值相比较为平缓，说明在采用Maxwell模型处理粒子与表面能量交换时使用的适应系数不应是固定值，而是应该随着角度进行变化。从粒子撞击的角度来说，正入射时粒子垂直撞击热流计表面，绝大部分能量转化为热能，而斜入射粒子存在一定的出射速度，从而降低了能量交换。因此，当角度升高时，适应系数应随之降低。
 |
| 图 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×104 m/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 |
| 图选项 |
5 结论 本文针对LIPS-200型离子推力器羽流热效应进行了三维数值模拟，采用准中性假设和Boltzmann关系式描述电子分布，采用Maxwell模型对粒子与表面能量交换过程，得到如下结论：
1) 采用混合PIC-DSMC方法可以对电推进羽流热效应进行有效的数值模拟分析，推力器出口轴线上滞止热流仿真与实验对比，误差小于17.0%。
2) 在进行径向热效应仿真时，由于粒子具有一定的出射速度，热适应系数不是恒定不变的，而是应随着角度的升高而减小。
3) 在进行三维粒子模拟时，模拟热流计的介入对流场整体影响较小。对Xe离子的影响主要集中在热流计后方0.1 m范围内，并出现了小范围内数密度为0的区域。热流计后方Xe离子数密度低于无热流计时的流场。
4) 模拟热流计对于Xe原子的影响同时存在于热流计前方和后方，在热流计前方中性粒子数密度升高约2个量级，与推力器出口附近的数密度相近。热流计后方Xe原子数密度低于无热流计时的流场。
综上所述，本文在验证羽流流场的基础上，对比了推力器出口1.0 m范围内的羽流热效应实验测量结果，验证了仿真模型精确性，可为进一步分析卫星及敏感仪器受热流的影响提供可信的输入条件。
参考文献
| [1] | KORKUT B, LEVIN D A, TUMUKLU O. Simulations of ion thruster plumes in ground facilities using adaptive mesh refinement[J]. Journal of Propulsion and Power, 2017, 33(3): 681-696. DOI:10.2514/1.B35958 |
| [2] | BROPHY J. Advanced ion propulsion systems for affordable deep-space missions[J]. Acta Astronautica, 2003, 52(2): 309-316. |
| [3] | HU Y, WANG J. Electron properties in collisionless mesothermal plasma expansion:Fully kinetic simulations[J]. IEEE Transactions on Plasma Science, 2015, 43(9): 2832-2838. DOI:10.1109/TPS.2015.2433928 |
| [4] | 周志雄, 魏蔚, 汪荣顺. 真空下气-固界面热适应系数的数值计算[J]. 低温与超导, 2007, 35(1): 36-40. ZHOU Z X, WEI W, WANG R S. Themathematic calculation of thermal accommodation coefficients at gas-solid interface in vacuum[J]. Cryogenics, 2007, 35(1): 36-40. DOI:10.3969/j.issn.1001-7100.2007.01.010 (in Chinese) |
| [5] | BERISFORD D F, BENGTSON R D, RAJA L L, et al. Heat flow diagnostics for helicon plasmasa[J]. Review of Scientific Instruments, 2008, 79(10): 10F515. DOI:10.1063/1.2955710 |
| [6] | 王黎珍, 史纪鑫, 郑世贵. 推力器真空羽流热效应计算模型修正及误差分析[J]. 航天器环境工程, 2014, 31(5): 483-488. WANG L Z, SHI J X, ZHENG S G. The error analysis and the improvement of heating effect model of the thruster vacuum plume[J]. Spacecraft Environment Engineering, 2014, 31(5): 483-488. DOI:10.3969/j.issn.1673-1379.2014.05.005 (in Chinese) |
| [7] | SHANG S F, CAI G B, ZHU D Q, et al. Design of double-layer anti-sputtering targets for plume effects experimental system[J]. Science China:Technological Sciences, 2016, 59(8): 1265-1275. DOI:10.1007/s11431-016-6037-y |
| [8] | HE B J, ZHANG J H, CAI G B. Research on vacuum plume and its effects[J]. Chinese Journal of Aeronautics, 2013, 26(1): 27-36. DOI:10.1016/j.cja.2012.12.016 |
| [9] | BIRDSALL C K, LANGDON A B. Plasma physics via computer simulation[M]. Bristol: ADAM Hilger, 1991. |
| [10] | BIRD G A. Molecular gas dynamics and the direct simulation of gas flows[M]. Oxford: Oxford University Press, 1994. |
| [11] | GARNER C E, POLK J R, BROPHY J R, et al.Methods for cryopumping xenon: AIAA-96-3206[R].Reston: AIAA, 1996. |
| [12] | BYOD I D. A review of hall thruster plume modeling[J]. Journal of Spacecraft & Rockets, 2000, 3(38): 381-387. |
| [13] | ZHENG H R, CAI G B, LIU L H, et al. Three-dimensional particle simulation of back-sputtered carbon in electric propulsion test facility[J]. Acta Astronautica, 2017, 132: 161-169. DOI:10.1016/j.actaastro.2016.12.016 |
| [14] | KIDD C T, NELSON C G.How the Schmidt-Boelter gage really works: CONF-9505201[R].Research Triangle Park: Instrument Society of America, 1995. |
| [15] | KELTNER N R.Heat flux measurements: Theory and applications[M].AZAR K.Thermal measurements in electronic cooling.Boca Raton: CRC Press, 1997: 273-320. |
| [16] | GIFFORD A, HOFFIE A, DILLER T, et al. Convection calibration of Schmidt-Boelter heat flux gages in shear and stagnation air flow[J]. Journal of Heat Transfer, 2010, 132(3): 031601. DOI:10.1115/1.3211866 |
| [17] | VLADIMIR K, ALEXANDER S, IGOR S.Investigation of the accelerated ions energy accommodation under their impingement with solid surfaces: AIAA-2002-4110[R].Reston: AIAA, 2002. |
| [18] | ZHANG T P.Initial flight test results of the LIPS-200 electric propulsion system on SJ-9A satellite[C]//33rd International Electric Propulsion Conference, 2013: 47-54. |
| [19] | VANGILDER D B, FONT G I, BOYD I D. Hybrid Monte Carlo-particle-in-cell simulation of an ion thruster plume[J]. Journal of Propulsion and Power, 1999, 15(4): 530-538. DOI:10.2514/2.5475 |
| [20] | CAI G B, LING G L, HE B J. An introduction to the novel vacuum plume effects experimental system[J]. Science China-Technological Sciences, 2016, 59(6): 953-960. DOI:10.1007/s11431-016-6024-3 |
