基于非结构网格的离子推力器羽流数值模拟方法研究

林骁雄,温 正,陶家生

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

基于非结构网格的离子推力器羽流数值模拟方法研究

林骁雄,温 正,陶家生

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

为适应复杂的航天器外形，文章采用基于非结构网格的数值模拟方法研究离子推力器真空羽流的基本特征。首先，介绍了非结构网格划分的方法，并分析了建模区域生成各种拓扑关系的特点；然后，提出了粒子穿越网格的具体算法，包括有利于减少计算量的结合矩阵重标号的带宽压缩算法；最后，开展了模拟结果和试验结果的对比研究，验证了以上数值模拟方法的可行性。

离子电推力器；羽流；数值模拟；非结构网格

0 引言

采用高比冲的电推进系统（如霍尔推力器、离子推力器）来实现星箭分离后的卫星变轨、入轨后的位置保持等任务，具有良好应用前景[1-2]。但是，电推进系统的应用中存在羽流污染问题。羽流中由于离子和中性原子碰撞所引起的低能量电荷交换（CEX），使离子速度降低，较易受到羽流自洽电场的影响，进而改变运动方向，撞击到太阳电池阵以及航天器表面等，产生污染。国内外学者均对如何减小羽流对航天器的污染影响进行过研究[3-8]。目前，数值模拟研究多采用二维轴对称网格或三维结构网格，较少采用非结构网格[9]。

结构网格中每一内部节点周围的网格单元数是完全相同的，并能通过任一节点立即找到与其邻近的节点。因此，结构网格用在推力器局部羽流特征分析时有较大优势。然而，多个电推力器的应用，加上复杂的航天器外形，使得航天器周围环境越来越复杂。对于多推力器阵列的羽流分析[10]、推力器羽流形成的航天器周围等离子体环境分析[11]等就不太适宜用结构网格的模拟。而非结构网格中内部节点周围的单元数是不固定的，更容易离散外形复杂的区域，对任意复杂外形的高度适体性是其他计算网格所无法比拟的。因此，本文采用基于非结构网格的方法模拟离子推力器羽流，并与试验结果对比，验证这一方法的可行性。

1 物理及数值模型

1.1 羽流物理模型

羽流等离子体中存在着大量的带电粒子，每个粒子同时受到其他粒子的作用，如洛伦兹力和碰撞等。本文主要对推力器喷射的羽流粒子在电场、磁场和碰撞作用下的运动状况进行数值模拟。对于电场的作用，需要计算等离子体运动产生的电势，并与外加电场产生的静态电势进行叠加以得到电势分布。

电场作用方程为

式中：ε0为真空介电常数；ρ为空间电荷密度；E为电场强度。

引入电位Φ，由于E=-∇Φ，得到泊松（Poisson）方程为

通过求解此方程可以得到空间电势，进而得到电场强度。

电磁力作用方程为

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

牛顿运动方程为

式中：P为动量，有P=mv，m为粒子质量。

离子电推力器产生的磁场在其外部很微弱，因此可以忽略磁场的影响[14]。当外部磁场较弱时，带电粒子的运动方程简化为

本文没有对电推力器羽流的碰撞模型、入口分布等问题进行详细讨论，具体分析可见文献[3-7]。

1.2 计算参数

本文以某卫星的LIPS-200氙离子推力器为研究对象。该推力器为双栅极系统，栅极表面向外凸出。推力器的主要工作参数如表1所示。模拟区域的边界采用固定电势边界。

表1 离子推力器工作参数Table 1 Ion thruster parameters

表1（续）

2 模拟计算

2.1 计算区域和计算网格

图1为卫星上离子推力器的安装位置示意图，有2个推力器固定在矢量推力机构上。

图1 卫星上推力器安装位置Fig.1 Installation location of thrusters on satellite

计算区域为包含推力器出口和中和器出口的一个立方体区域，其网格采用Gmsh软件通过脚本语言生成，如图2所示。推力器加速栅极曲面的弧度为15°。建立网格的目的：一方面是便于粒子碰撞时的采样，另一方面便于应用 PIC（Particle in Cell）方法求解空间电势。结构网格可以方便建立泊松方程的差分格式，但是对于复杂边界，比如极小的中和器出口和较细的天线结构，会遇到网格生成程序无法收敛的情况；另外在结构网格中确定粒子所在网格单元需要花费较大计算量。非结构网格能自适应生成符合要求的贴体网格，而且根据网格面片和体积单元的链接关系可以方便确定粒子所在的网格单元。因此，本文采用了非结构网格，并在电势计算中将体积网格统一选作四面体网格。选取四面体网格的原因是一方面它可以模拟各种几何形状，另一方面在四面体网格中可以非常方便地实施粒子穿越网格算法。这个算法是保证在非结构网格中粒子运动定位的关键步骤，如果选用其他类型的非结构网格可能导致这个关键步骤无法实施，具体分析见2.3节。

图2 模拟区域和推力器及中和器局部网格Fig.2 Simulation domain and local grids around thruster and neutralizer

2.2 网格预处理

网格预处理主要用来生成非结构网格间各种层次关系，而层次关系是指求解区域（Group）、四面体单元（Cell）、三角形表面（Face）、棱边（Edge）、节点（Node）等拓扑元素的链接关系。在羽流的数值模拟中主要包含以下几种：

1）在缩减矩阵带宽中需要用到计算区域内边到节点的映射Edge2Node。

2）在求解电势的有限元计算中需要用到四面体单元到节点的映射Cell2Node、推力器表面三角形单元到节点的映射 ThrFace2Node、航天器表面三角形单元到节点的映射ScFace2Node，以及边界表面三角形单元到节点的映射 BdFace2Node等。其中，ThrFace指推力器表面三角形单元、ScFace指航天器表面三角形单元、BdFace指模拟区域边界表面三角形单元。这些映射已经由Gmsh生成网格时提供。

3）在粒子穿越网格算法中需要用到四面体单元到三角形单元的映射Cell2Face、三角形单元到四面体单元的映射Face2Cell、表面三角形单元局部标号到全局三角形单元标号的映射 ThrFace2Face、ScFace2Face、BdFace2Face。这些拓扑关系数组在结果后处理过程中也会用到，提前生成这些拓扑关系数组就可以在粒子运动过程中直接查表，极大提高了计算速度。拓扑关系如图3所示。

图3 网格间的拓扑关系Fig.3 Topological relationship among grids

2.3 粒子穿越网格算法

粒子穿越网格是非结构网格中模拟粒子运动的关键步骤之一。在结构网格中，可以借助网格下标直接索引粒子将要穿越的相邻网格，而在非结构网格中只能凭借网格层级关系进行搜索。

粒子穿越网格表面的过程和表面的特征关系密切。如果粒子穿越的是计算区域内用于计算的辅助网格表面（如图4中的粒子3），将不对粒子的运动产生任何影响；如果粒子穿越的是删除表面（如图4中的粒子4），则粒子运动到此表面后，在维护粒子的列表中将被删除；如果粒子穿越的是反弹表面（如图4中的粒子5），则粒子运动到此表面之后，将和壁面发生碰撞，反弹之后仍位于计算区域内。本文还提出了一种粒子更新算法，如图5所示。该更新算法主要根据粒子当前位置、粒子在模拟步长内移动的距离和当前单元边界的关系进行判断，粒子在模拟步长内移动的距离通过公式(5)得到的速度对时间积分而得到。

图4 粒子穿越网格模型Fig.4 Model of particle traversing through the meshs

图5 粒子运动更新算法Fig.5 Update algorithm of particle motion

2.4 Poisson方程的求解

采用有限元方法求解 Poisson方程（式(2)），形成的代数方程为Ax=b，其中A为系数矩阵，x为待求解各节点的电势，b为结合边界条件形成的矢量。最终形成的系数矩阵A为对称和稀疏的矩阵，其中稀疏矩阵是指矩阵中绝大部分元素皆为0。本文采用带状压缩存储的方法存储矩阵A中的上三角部分。为了使矩阵重排序后的带宽及轮廓较小，借助于RCM算法[12]寻找这样的一个排序数组。RCM示例网格如图6中左图，图中4、5为内部节点，1、2、3、6、7、8为边界节点，将各个节点和边形成的无向图定义为图G，算法过程简单概括为：

1）找出图G的边界顶点r，并赋予x1（在此示例中r取为节点1）；

2）以顶点r为根，生成图G的层次结构{x1,x2, …,xn}（如图中虚线区分的3个层，括号中的数值为节点在该层中的度）；

3）对i=1, 2, …,n，找出图G的层次结构{x1,x2, …,xn}中顶点xi的所有未编号相邻顶点，并按其度的增加顺序编号；

4）将上述顶点编号按逆序排列，得到RCM排序数组（如图6中右图所示）。

图6 RCM算法示例网格Fig.6 Grids of RCM algorithm

使用RCM前、后的非零矩阵元素及半带宽如图7所示。

图7 使用RCM算法前、后的非零矩阵元素及半带宽Fig.7 Non zero elements and semi-bandwidth before and after using RCM algorithm

可以看到，原来矩阵的半带宽为 5，变换后的半带宽为 4。由于本计算为示例，半带宽并没有明显减少。但是在实际的应用中，因为内部节点数远远大于边界节点数，半带宽可以减少到原来的 1/5甚至更多。采取这种方法形成矩阵可以极大减少方程的运算次数，而且在用稀疏带状矩阵进行Cholesky分解时引入的非零元素都在带宽及轮廓范围内，不会过多占用额外内存，可以提高运算效率。

3 计算结果及对比分析

图8给出了计算所得的束流离子数密度分布，基本分布在轴线±15°附近，而边缘附近的锯齿状图案是由于粒子数偏少导致的局部统计涨落。

图8 计算所得离子数密度分布Fig.8 Simulation results of ion density distribution

图9给出的是本文模拟计算结果与文献[13]试验测量结果的比较。可见离子数密度在不同径向距离处其值都随轴向距离的增加而减小；在靠近推力器中心线位置处，数密度值的变化幅度较大，在远离推力器中心线位置处，变化幅度较小。在径向距离140 mm的曲线和径向距离180 mm的曲线上均出现了明显偏离主束流特征的数值，说明这些位置点已经不在主束流区。模拟结果和试验结果的趋势与数值范围基本一致。

图9 离子数密度的轴向分布的模拟结果和试验结果对比Fig.9 Comparison of axial profiles of ion density between simulation and experimental results

图10给出的是在整星状态下计算得到的羽流电势分布，左侧为卫星，右侧为太阳电池阵。对于这种复杂结构组合形成的计算区域中，非结构网格方法更显示其优势。

图10 羽流电势Fig.10 Plume potential

4 结束语

本文介绍了采用三维非结构网格对离子推力器羽流进行数值模拟的方法。模拟结果与试验测量数据进行了比较，数据基本一致，证明了基于非结构网格方法进行羽流模拟的可行性。本方法在求解边界情况更加复杂的航天器构型的周围等离子体环境时会更有优势。

（References）

[1]周志成, 高军.全电推进 GEO卫星平台发展研究[J].航天器工程, 2015, 24(2): 1-6 ZHOU Z C, GAO J.Development approach to all-electric propulsion GEO satellite platform[J].Spacecraft Engineering, 2015, 24(2): 1-6

[2]胡照, 王敏, 袁俊刚.国外全电推进卫星平台的发展及启示[J].航天器环境工程, 2015, 32(5): 566-570 HU Z, WANG M, YUAN J G.A review of the development of all-electric propulsion platform in the world[J].Spacecraft Environment Engineering, 2015, 32(5): 566-570

[3]孙安邦, 毛根旺, 陈茂林, 等.离子推力器羽流特性的粒子模拟[J].强激光与粒子束, 2010, 22(2): 401-405 SUN A B, MAO G W, CHEN M L, et al.Particle simulation of ion thruster’s plume characteristics[J].High Power Laser and Particle Beams, 2010, 22(2): 401-405

[4]李娟, 楚豫川, 曹勇.离子推力器羽流场模拟以及Mo~＋CEX 沉积分析[J].推进技术, 2012, 33(1): 131-137 LI J, CHU Y C, CAO Y.Numerical simulations of ion thruster plume contamination interactions on spacecraft[J].Journal of Propulsion Technology, 2012, 33(1): 131-137

[5]任军学, 顾左, 郭宁, 等.离子发动机羽流特性的数值模拟[J].航空动力学报, 2013, 28(6): 1372-1379 REN J X, GU Z, GUO N, et al.Numerical simulation of characteristics of ion thruster plume[J].Journal of Aerospace Power, 2013, 28(6): 1372-1379

[6]任军学, 王艳, 仇钎, 等.离子发动机羽流二维轴对称数值模型与验证[J].北京航空航天大学学报, 2011, 37(12): 1498-1503 REN J X, WANG Y, QIU Q, et al.Validation of two-dimensional axi-symmetric plume model for ion thruster[J].Journal of Beijing University of Aeronautics and Astronautics, 2011, 37(12): 1498-1503

[7]HU Y, WANG J J.Electron kinetic characteristics in plasma plumes: fully kinetic simulations[C]//53rdAIAA Aerospace Sciences Meeting.Kissimmee, Florida, 2015: 1-10

[8]王虹玥, 贺碧蛟, 蔡国飙.稳态等离子体推力器羽流仿真及其对微波的衰减和相移作用[J].推进技术, 2011, 32(6): 839-844 WANG H Y, HE B J, CAI G B.Numerical simulation of the stationary plasma thruster plume and its effects on the microwave attenuation and phase shifts[J].Journal of Propulsion Technology, 2011, 32(6): 839-844

[9]王学德, 伍贻兆, 夏健.三维非结构网格 DSMC方法的实现及其应用[J].计算力学学报, 2007, 24(2): 231-235 WANG X D, WU Y Z, XIA J.Implementation of 3D unstructured DSMC method and its application[J].Chinese Journal of Computational Mechanics, 2007, 24(2): 231-235

[10]FOSTER J E, PATTERSON M, PENCIL E, et al.Neutralizer characterization of a next multi-thruster array with electrostatic probes[C]//42ndAIAA/ASME/ SAE/ASEE Joint Propulsion Conference.Sacramento, California, 2006: 1-21

[11]杨集, 陈贤祥, 周杰, 等.伸杆对星载电场探测仪的影响研究[J].电子与信息学报, 2009, 31(4): 1010-1012 YANG J, CHEN X X, ZHOU J, et al.Study of the influence of booms on spaceborne electric field sensors[J].Journal of Electronics and Information Technology, 2009, 31(4): 1010-1012

[12]于二青, 王春江, 赵金城.矩阵重排序算法在结构分析快速求解中的应用[J].空间结构, 2010, 16(1): 45-50 YU E Q, WANG C J, ZHAO J C.Application of improved RCM algorithm in fast solution of structural analysis[J].Spatial Structures, 2010, 16(1): 45-50

[13]张尊, 汤海滨.氙气离子推力器束流等离子体特征参数的 Langmuir单探针诊断[J].高电压技术, 2013, 39(7): 1602-1608 ZHANG Z, TANG H B.Single Langmuir probe diagnostics for characteristic parameters of the beam plasma in Xe ion thruster[J].High Voltage Engineering, 2013, 39(7): 1602-1608

[14]任军学, 李娟, 仇钎, 等.离子发动机交换电荷离子返流的粒子模拟[J].强激光与粒子束, 2011, 23(7): 1929-1934 REN J X, LI J, QIU Q, et al.Particle simulation of charge-exchange ions back flow in ion thruster plume[J].High Power Laser and Particle Beams, 2011, 23(7): 1929-1934

（编辑：肖福根）

A method for numerical simulation of ion thruster's plume based on unstructured mesh

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

A method of particle simulation based on unstructured mesh is proposed in this paper to reveal the basic characteristics of ion thruster’s plume with complex spacecraft configurations.The method of grid discretization and the topological relationship among the generation of grids are discussed.An algorithm of particle motion is developed to push the particles in grids, then the method of matrix permutation based on the RCM algorithm is used to reduce the bandwidth of the matrix.Finally, a simulation example is compared with experiment results to validate this method.

ion electric thruster; plume; numerical simulation; unstructured mesh

V439.2

：A

：1673-1379(2016)05-0484-06

10.3969/j.issn.1673-1379.2016.05.005

林骁雄（1990—），男，硕士研究生，研究方向为电推进系统总体设计。E-mail: linxiaoxiong_1990@163.com。

2016-03-16；

：2016-09-03