林骁雄，陶家生，温 正

（中国空间技术研究院通信卫星事业部，北京100094）

离子推力器羽流沉积对卫星热控影响研究

林骁雄，陶家生，温 正

（中国空间技术研究院通信卫星事业部，北京100094）

离子推力器的羽流是等离子体，等离子体的组成是带电粒子，这与传统的化学推进系统的羽流成分有很大不同，带电粒子有在卫星表面吸附的倾向，会形成羽流沉积污染。这种羽流沉积会改变卫星表面的吸收率和发射率，从而影响卫星的热控性能。为了预测离子推力器的羽流对卫星的热控性能的影响，建立了离子推力器羽流模型。所建模型采用了工程化离子推力器的在卫星上的布局位置和离子推力器的工作参数，模拟了离子推力器的正离子与中和电子束在工程化中分置的实际情况，使模型更为符合实际。通过数值模拟得到了离子、电子、中性粒子的空间分布，电场分布，得到了钼粒子在卫星表面的分布及沉积厚度，比较了模型计算的离子分布与实验获得的离子分布情况，说明了模型分析的正确性，给出了卫星表面热性能的变化及局部区域温升的最大包络可达二十多度的结果。

离子推力器；羽流建模；卫星热控；羽流沉积

0 引言

电推力器在卫星上作为主推力器在国外已经在轨应用[1-2]，国内已有化学推力器与电推力器混合应用的卫星，全电推进卫星也已在研[3]。总之电推力器应用于卫星，对卫星在轨工作环境的影响情况仍需研究，这是因为传统的化学推力器燃烧产物在理论上是电中性的气体和水[4]，它们没有在卫星表面易于吸附的理论基础，但电推力器的羽流产物是等离子体，是电离态的带电粒子，理论上有在卫星表面吸附的倾向。尤其工程实际的离子电推力器正离子束与电子束是分立设置，工作时通过分立的电子束与离子束的汇聚形成宏观电中性的等离子体，汇聚的过程及结果有待研究[5-6]。离子推力器的离子束是在束流系统的引导下射出的，如果离子束流不能完全聚焦，离子束流将侵蚀离子推力器的钼栅极组件而溅射出钼粒子[7]，在羽流中携带有钼的中性粒子，这些中性粒子一部分通过离子交换电荷碰撞生成钼离子，一部分通过与电子直接碰撞离化形成钼离子[8]。在各种粒子碰撞和推力器羽流局部电势的作用下返流到卫星表面形成污染而影响卫星热控性能[9]。羽流的等离子体不仅可以在卫星表面吸附，它的离子在卫星表面也可能激发出二次溅射，更为有害。文献 [8]中通过数值模拟得到了钼原子的空间分布，文献 [13]中研究了均匀厚度钼膜的光学特性，文献 [7]结合试验并根据经验模型得到了栅极材料的原子溅射率，但是没有实现从钼原子的出射到钼原子的沉积直到对卫星OSR（Optical Solar Reflector） 片热控性能影响的完整分析。为了研究离子电推力器羽流中各种成分的相互作用和各种粒子的空间分布情况、电场分布情况，本文以某实践卫星的LIPS-200氙离子推力器应用为研究对象，结合离子推力器在卫星中的布局参数和工作参数建立了离子束与电子束分立设置的工程模型，研究包含电子束与离子束的汇聚并形成宏观电中性的等离子体的过程，得到离子、电子、中性粒子的空间分布和电场分布，以及钼粒子在卫星表面的分布及沉积厚度，预示卫星表面热性能的变化，从而使模型的建立更加贴近工程实际，得到客观的分析结果。

1 离子电推力器羽流场实验研究的难点

卫星的热性能在地面要通过热真空实验加以确定，电推力器羽流中钼原子沉积对卫星热控性能的影响也可以通过地面热真空实验加以研究确定，但是这种地面热真空实验很难测得准确的结果。原因是，要想实验模拟太空中的真空环境需要至少10-3Pa的真空度，热真空实验要在密闭的真空罐内进行，空罐内体积有限，一旦电推力器点火喷出羽流，真空罐内的真空度难以维持，而电推羽流打到罐壁上产生的溅射对类似OSR片这种光学表面性能的影响无法回避[10]，因此，通过地面实验准确确定电推羽流对OSR片光学表面性能的定量影响难于实现。所以，本文拟采用数值建模模拟的方法，辅以实验结果在一定程度上考察模型的正确性，来研究电推羽流对卫星热控性能的影响。

2 卫星表面的光学性质

为了建立卫星内部仪器的工作环境，需要对卫星进行热控。一般卫星热控是以被动热控为主，主动热控为辅，即通过调控卫星表面对阳光的吸收率和其红外发射率来建立起卫星内部温度与卫星在轨的空间低温环境的平衡[11]，这就是被动热控的基本思想。主动热控是在被动热控的基础上，提高或降低局部温度。可见，卫星表面的发射率和吸收率的改变，将是对卫星热控基础的改变，其影响是长期的和单向的，需要给予足够重视。

现有卫星的热控表面主要有两种，一是热控多层材料，二是OSR片。OSR片（光学太阳反射镜）是一种以玻璃为透明薄膜表面层的第二表面镜型热控涂层，其薄膜表面层可用石英玻璃或铈玻璃，金属底层由玻璃背面真空蒸发-沉积得到，通常为银或铝。热控多层主要是绝热作用，OSR片是一种发射率和吸收率可设计的光学表面，是通过发射率和吸收率的设计来改变卫星内部热的流出和外部热的流入，来调节卫星的内部温度。如果热控涂层性能退化，卫星入轨后温度升高，用于卫星温度场数值计算的参数设计初值与参数的实际值之间就会存在一定程度的偏差，以致卫星不能正常工作[9]。发射率和吸收率是在OSR片制造中即确定了的，卫星热控则通过改变OSR片的铺设面积调节卫星的内部温度。可见，OSR片是卫星内部温度调节的热通道，而电推力器的安装位置一般也恰好在OSR片附近，这也在一定程度上加重了电推力器羽流对卫星热控性能的影响。

随束流飞出的钼粒子，部分沉积在卫星表面及太阳翼表面，成为电推进羽流沉积污染的主要因素。从对电推力器羽流的初步分析可知[12]，电推力器羽流的钼粒子沉积厚度为几纳米到几百纳米，恰和光学膜系厚度的量值相当，可能会改变光学膜系的性质。因此，需要深入研究热控涂层性能在电推力器羽流钼原子沉积的影响下的退化规律。

3 电推力器羽流模型的建立

3.1 羽流建模的理论基础

推进剂在离子发动机电离室中电离，在栅极的加速作用下高速喷出，并由中和器发射电子中和，以保持整体电中性，因此离子发动机羽流为等离子体射流[14-16]。部分未电离的推进剂以自由分子流的形式从栅极流出，在束流区域与相对高速的束流离子会发生交换电荷碰撞，产生低速交换电荷离子。同时由于交换电荷离子会撞击栅极造成栅极溅射腐蚀，产生非推进剂粒子。

离子发动机的羽流中，推进剂电离形成的束流离子在栅极加速作用下喷出速度一般在20 000～40 000 m/s，是羽流的主要组成部分，其运动轨迹呈射流状，受电磁场的影响很小。下游区域的束流扩张主要由航天器的电势与空间电磁场环境决定，实验测试的数据结果表明，大部分束流离子都集中在扩张半角为15～20°的范围内[17]。为了保持总体准电中性，束流离子会产生一个径向的电场，防止电子离开束流区。尽管这个电场对于高速束流离子的运动没有影响，但是会对低速的交换电荷离子产生径向加速作用并形成交换电荷离子回流效应，从而对航天器和太阳能电池帆板产生影响。

带电粒子在电磁场中的受力根据洛伦兹方程得：

式中：F为电磁力；q为粒子电荷；v为速度；B为感应磁场。

磁场较弱时，忽略磁场的影响，则受力方程简化为：

钼的电离能 (Mo→Mo+)是7e V左右，羽流中的电子和离子的动能可以用上述理论方程通过建模运算得到，并建立起钼的电离条件。

3.2 羽流模型的程序化实现

分析计算的时间框架分为如图1所示的六个部分：

1） 导入控制文件：相关控制参数导入计算空间；

2） 导入网格文件：梳理“**.msh”文件的数据格式，把相关原始数据导入计算空间；

3）网格文件预处理，产生各级网格链接关系；

4）粒子初始化：基于粒子存储的数据结构，初始化分析计算所用的各类粒子链表；

5） 主循环：在前期各种控制参数和网格数据准备好后，进入主循环。在每个时间步中，更新粒子的位置，计算电势和电场，计算粒子的碰撞，并处理边界条件。

6） 结果输出：根据控制参数要求的输出种类，逐个输出结果文件。

图1 数值模拟的时间框架Fig.1 Time frame of numerical simulation

分析计算的逻辑框架如图2所示：

图2 数值模拟的逻辑框架Fig.2 Logical frame of numerical simulation

3.3 羽流建模中的工程化参数

建立离子推力器在某实践卫星中的配置模型（图3）。图中显示了推力器与OSR片的相对位置，推力器轴线指向卫星质心。推力器为双栅极系统，栅极表面向外凸出。

图3 电推力器在某卫星中的空间位置Fig.3 Position of thruster in a satellite

按照离子推力器在某实践卫星的空间位置，建立了离子推力器卫星的局部模型（图4），在该模型的建立中主要位置关系参数是离子推力器与卫星表面的位置、角度关系。

图4 电推力器与OSR片的模型Fig.4 Model of electric thruster and OSR

在卫星离子推力器的建模参数中除了空间位置关系外还有推力器的主要工作参数（如表1所示）。模拟区域的边界采用等电势的固定电势边界，为了弱化等电势边界对模型运算结果的影响，模型中加大了计算的空间区域。

表1 离子推力器工作参数Tab.1 Working parameters of ion thruster

3.4 羽流模型的运算结果

在上述模型和初始参数和边界条件的情况下，钼原子生成钼离子有两种方式，一种为电子碰撞电离，一种为和高速运动的氙离子发生电荷交换碰撞。由于电子的能量相对较高，而钼原子流的粒子数密度相对较小，后者生成的钼离子远小于前者。本文采用质点网格法（PIC，Paritcle In Cell） 和直接蒙特卡洛方法（DSMC， Direct Simulation Monte Carlo） 自编软件进行计算，得到如下结果。

图5是推力器近场局部电势分布图，图中最高电势位于推力器喷口前方为21.5 V，最低为-2 V。

图5 推力器局部电势分布图(左图为侧视，右图为俯视，下同)Fig.5 Partial potential distribution around thruster

图6 电推力器羽流离子数密度分布图Fig.6 Plume ion number density distribution around thruster

图6为离子数密度的空间分布，最高为6.33×1015个/m3，集中在推力器喷口前方附近，最低为0。图7为钼离子数密度的空间分布，最高为1.09×1011个/m3生成的钼离子在径向电场的作用下向推力器背部扩散，具有很大的返流倾向。由于采用四面体非结构网格进行计算，粒子流边缘处数值变化差异较大，粒子数在统计过程中也会有一定涨落和随机性，所以出现了锯齿状，下同。

图7 钼离子数密度空间分布Fig.7 Spatial distribution of Mo ions

图8是钼原子数密度的空间分布，最高为2.32×1012个/m3，集中在推力器喷口前方区域，最低为0。由于是中性粒子，仅受粒子碰撞扩散的影响，趋于与推力器羽流的方向一致。结合目前的分析，到达卫星的表面极少为钼原子，更多的是钼离子。

图8 钼原子数密度空间分布Fig.8 Spatial distribution of Mo atoms

图9是OSR片附近钼离子数密度分布，最高为4.19×108个/m3，集中在推力器附近的卫星表面，最低为0。

图9 OSR片附近钼离子数密度分布Fig.9 Spatial distribution of Mo ions around OSR

4 数值分析结果与实验结果的比较与分析

由于钼粒子空间分布较准确的实验结果不易得到，这里将电推力器羽流的离子分布模拟结果与实验结果进行比较来检验模型的正确与否。

图10为离子数密度模拟结果与试验结果对比，可以看出两者相当一致。由于实验在先，没有建模分析的预示，实验中布置测点较少，所以这里也引用文献 [17]的结果给出一个较为细致的比较。

图10 离子数密度模拟结果与试验结果对比Fig.10 Comparison between simulation results and experiment results of ion density

在图11中可见离子浓度在不同径向距离处的参数值都随轴向距离的增加而减小；在靠近推力器中心线位置处，参数值的变化幅度较大，在远离推力器中心线位置处，参数值的变化幅度较小。在径向140 mm曲线和径向180 mm曲线上均出现了明显偏离主束流特征的参数值，说明这些位置的点已经不在主束流区。模拟结果和试验结果的趋势和参数值范围基本一致，在此环节证明了建模分析的有效性。

图11 离子浓度的轴向分布的模拟结果和实验结果对比Fig.11 Comparison between simulation and experiment results of axial distributions of ion density

5 卫星表面热性能改变预示

通过计算得到钼离子在电推力器附近卫星表面10 000 h的沉积厚度如图12所示。可见，由于推力器两侧与推力器出口轴线的夹角相对较小，在推力器两侧会有较大量的沉积，沉积率最大值量级为2.2×10-7埃/s，10 000 h沉积量最大值为7.912埃，而在推力器安装位置的背面（即图3中下部中心位置处）会有较少的沉积。

图12 OSR表面钼离子10 000 h沉积厚度分布（已取对数，埃）Fig.12 Deposition thickness of Mo ions on OSR surface in 10 000 h

参阅文献 [11]得到钼膜玻璃数据，见表2。

表2 钼膜玻璃的光学性质Tab.2 Optical properties of molybdenum-film glass

结合上述分析，得到了基板的光学性质和钼膜的光学性质，复合膜的光学性质为复合模型每层膜的光学性质的组合。经过薄膜的光线会被薄膜反射和吸收，同样也会发生在玻璃-薄膜界面再次发射、透射、吸收。结合典型的OSR参数（吸收率为0.12，发射率为0.78）计算钼的复合膜特性，得到表3。

通过线性插值，得到OSR片复合膜吸收率分布如图13所示。可以看到，散热板OSR片上部分区域的吸收率有所上升，在推力器附近上升最为明显，这将引起卫星散热板平衡温度的上升。

表3 钼的复合膜特性Tab.3 Properties of molybdenum composite film

图13 复合膜吸收率Fig.13 Absorptance of composite film

结合已经得到的复合膜的光学性质对卫星的热平衡进行分析。OSR表面接收到的外热流有太阳辐射、地球红外辐射、地球反射辐射以及舱体辐射4种热流。忽略地球红外辐射、地球反射辐射以及舱体辐射，认为OSR表面接收到的外热流只有太阳辐射，OSR发出的热流主要是对外的热辐射。卫星南北板散热能力按照q=300 W/m2计算。根据能量守恒定律，可以列出卫星舱板温度计算表达式为

式中：传入的外热流量Q1=αS，其中α为膜的吸收率，太阳常数S为1 353 W/m2；内热源Q2=q；向深冷空间辐射散热量Q3=εσT4，其中ε为膜的发射率，斯忒藩-波尔兹曼常数σ为5.67×10-8W/（m2·K4），T为平衡温度。

将各热量的计算式代入上式，整理后得电池板能量守恒方程为

代入ε＝0.78，α＝0.12，得到卫星OSR片初始平衡温度为320 K。

由于真实卫星还有辐射换热等热导途径，所以实际结果会与这里的计算有差别，这里相当于绝热情况下的温度，目的是揭示羽流对热控的影响情况。

图14为离子电推力器工作10 000 h后卫星舱板上钼粒子沉积之后的温度升高情况，可以看到最高温度变化包络为25℃，绝大部分区域温度升高不超过10℃。由于一般仪器设备的环境温度设计余量为±11℃[11]，可见电推力器羽流将对卫星寿命末期的整星热控设计产生不利影响。因此，对温度变化较为敏感的设备应避免安装在推力器周围，另外通过采用局部加强导热等热控措施可以缓解局部小区域的温度升高。

图14 OSR片温度变化分布Fig.14 Distribution of OSR temperature variation

6 总结与讨论

为了研究离子推力器在卫星上的应用可能，对卫星在轨热性能的影响进行了建模分析。建模中尽可能采用了工程化参数，并且针对离子电推力器的离子流和中和电子流分置的实际情况进行了二者的分立式建模。建模中模拟了正负粒子流汇流成羽流的过程，成功获得了离子推力器羽流的工质离子、电子、中性粒子、钼离子的空间分布，电场分布，得到了钼粒子在卫星表面的分布及沉积厚度。比较了模型计算的离子分布与实验获得的离子分布情况，说明了模型分析的正确性，给出了卫星表面热性能的变化及局部区域温升最大可达二十多度的结果。预示了离子推力器卫星在寿命末期应关注羽流产生的对热控的不利影响。由于推力器两侧与推力器出口轴线的夹角相对较小，在推力器两侧会有较大量的沉积。此区域的热平衡温度会有较大升高，最大温升包络可能在二十摄氏度量级，一些对温度变化较为敏感的设备应避免安装在推力器周围。通过采用局部加强导热等热控措施可以缓解局部小区域的温度升高。

[1]DAVID J A,ERIC P,DANIEL V,et al.Products from NASA's in-space propulsion technology program applicable to low-cost planetary m issions[J].Acta astronautica,2014(93):516-523.

[2]胡照,王敏,袁俊刚.国外全电推进卫星平台的发展及启示[J].航天器环境工程,2015,32(5):566-570.

[3]周志成,高军.全电推进GEO卫星平台发展研究[J].航天器工程,2015,24(2):1-6.

[4]方杰,田辉,蔡国飙.N20单组元微推进系统及其喷管流场的初步研究[J].推进技术,2005,26(6):495-498.

[5]操慧珺,楚豫川,曹勇,等.离子推进器中离子束中和及推进过程的全粒子数值模拟[J].高电压技术,2014,40 (7):2119-2124.

[6]YUAN H,JOSEPH J W.Electron kinetic characteristics in plasma plumes:fully kinetic simulations[C]//Proceedings of 53rd AIAA Aerospace Sciences Meeting. Kissimmee,Florida:AIAA,2015:1-10.

[7]贾艳辉,李忠明,张天平,等.栅极系统电子返流对离子推力器寿命影响[J].真空与低温,2012,18(1):21-25.

[8]ROY R I S,HASTINGS D E,GATSONIS N A.Nume-rical study of spacecraft contamination and interactions by ion-thruster effluents[J].Journal of spacecraft and rockets,1996,33(4):535-542.

[9]马伟,宣益民,韩玉阁,等.长寿命卫星热控涂层性能退化及其对卫星热特征的影响[J].宇航学报,2010,31(2): 568-572.

[10]张天平,陈继巍,李小平,等.地面寿命试验中离子推力器表面的溅射沉积量计算 [J].航天器环境工程,2011,28 (5):4 36-439.

[11]侯增祺,胡金刚.航天器热控制技术-原理及其应用[M].北京:中国科学技术出版社,2007(1):368-370.

[12]KEMP R F,LUEDKE E E,HALL D F,et al.Effects of electrostatic rocket material deposited on solar cells[C] //Proceedings of the 9th NASA Electric Propulsion Conference.Washington D C:NASA,1972:1-10.

[13]范平,邵建达,易葵,等.纳米Mo膜的光学特性及最小连续膜厚研究[J].中国激光,2005,32(7):977-981.

[14]孙安邦,毛根旺,陈茂林,等.离子推力器羽流特性的粒子模拟[J].强激光与粒子束,2010,22(2):401-405.

[15]李娟,楚豫川,曹勇.离子推力器羽流场模拟以及Mo～＋CEX沉积分析[J].推进技术,2012,33(1):131-137.

[16]任军学,顾左,郭宁,等.离子发动机羽流特性的数值模拟[J].航空动力学报,2013,28(6):1372-1379.

[17]张尊,汤海滨.氙气离子推力器束流等离子体特征参数的Langmuir单探针诊断[J].高电压技术,2013,39(7): 1602-1608.

(编辑：王建喜)

Influence of plume deposit of ion thruster on thermal control properties of satellite

LIN Xiaoxiong,TAO Jiasheng,WEN Zheng
（Institute of Telecommunication Satellite,China Academy of Space Technology,Beijing 100094,China）

The ion electric thruster plume is plasma that consists of charged particles,which has a tendency to be adsorbed onto satellite surface.The deposition of the plume can change the absorptivity and emissivityof the satellite surface,which has a negative influence on thermal control property of the satellite.In order to predict this influence,a ion thruster plume model was built.The layout position and working parameters of an engineering ion thruster are adopted in the model to simulate the actual conditions of positive ion and electron beam,which makes the model more in accord with reality.The spatial distribution and electric field distribution of ions,electrons and neutral particles,and deposition distribution of molybdenum particles on satellite surface were obtained with numerical simulation.The ion distributions got by simulation calculation and experiment was compared,which validated the correctness of the model analysis.The thermal property of the satellite surface and the maximumenvelope oftemperature increase in local region ofthe satellite surface are given in this paper.

ion thruster;plume modeling;satellite thermal control;plume deposition

V434-34

A

1672-9374(2017)02-0009-09

2016-04-15；

2016-12-22

林骁雄（1990—），男，硕士研究生，研究领域为电推进总体技术