碘工质霍尔推力器内部过程数值模拟

2022-04-28梁子轩王平阳徐宗琦杭观荣

上海航天 2022年2期

王尚，梁子轩，王平阳，徐宗琦，杭观荣

（1.上海交通大学机械与动力工程学院，上海 200240；2.上海空间推进研究所上海空间发动机工程技术研究中心，上海 201112）

0 引言

霍尔推力器具有比冲高、推力精确可调、推力功率比大、寿命长等优势，能够给航天器带来质量、寿命、经济等增益，作为航天器的推进装置得到了广泛应用。随着霍尔推进技术的不断发展，碘作为氙的替代工质，因其自身升华温度低、电离能低、相对原子质量高、储量丰富、价格低廉等特点，得到了广泛研究。由于碘电离后，等离子体组分不同于氙，导致电离区位置和加速效应均不同，为提高推力器的性能，需重点研究其内部电离和加速过程。通过数值模拟方法可以得到推力器内部等离子特性的详细参数分布，得到放电通道内部电势、等离子体数密度、电子温度等参数，以及推力器的性能参数，进而深入分析放电通道内部等离子体的生成机理。上述工作为推力器性能和寿命的优化提供了详细的理论支撑，具有重要参考价值。

对于碘工质推力器放电过程的数值模拟，国内外学者均开展了相关研究。法国的GRONDEIN等采用流体动力学模型对碘工质离子推力器进行了数值模拟，获取了离子、电子数密度和电子温度等特征参数，将碘和氙工质情况下的推力、比冲和效率进行了比较；NASA 的MARIA采用混合单元粒子法模型对200 W 霍尔推力器进行了模拟，对比了碘和氙工质情况下羽流区的中性气体和离子的数密度分布；法国的LUCKEN 等针对离子推力器进行数值模拟，研究了沿推力器轴向方向的磁场强度变化时，碘工质和氙工质情况下推力器的推力和效率变化；哈尔滨工业大学的牛翔等采用一维流体动力学模型对碘工质会切场推力器进行数值模拟，研究了碘工质会切场推力器内部的等离子体参数分布。相较于上述研究成果，目前对百瓦级碘工质霍尔推力器放电通道内部过程的数值模拟研究工作还鲜有报道。

为了配合200 W 级碘工质霍尔推力器的设计优化和实验研究工作，本文基于单元粒子法/直接模拟蒙特卡罗法/蒙特卡罗碰撞模型（Particel-in-Cell/Direct-Simulation-Monte-Carlo/Monte-Carlo-Collision，PIC/DSMC/MCC）混合方法，对一定结构的200 W 级碘工质霍尔推力器进行了二维数值模拟仿真研究，得到了设计工况下推力器放电通道内部等离子体的各项参数分布；根据数值模拟的结果分析推力器的内部过程，和氙霍尔的等离子体参数分布进行对比。此项研究工作可为百瓦级碘工质霍尔推力器的性能、寿命优化、地面实验的开展提供理论依据和技术参考。

1 数值仿真

1.1 计算基本模型

采用PIC/DSMC/MCC 方法对推力器的内部过程进行数值仿真，PIC 用于求解自洽电场和计算带电粒子的加速度；DSMC 用于求解相对速度较小的中性分子-离子之间、中性分子-中性分子和离子-离子之间的碰撞，主要处理碘原子、碘离子之间的动量交换和电荷交换碰撞；MCC 用于求解电子参与的碰撞过程，主要处理电子和碘分子、碘原子之间的碰撞，三者结合使用完成粒子间的碰撞输运和运动过程。文献［10］已对此混合方法的可行性和正确性进行了验证，在此基础上继续开展工作。

基于PIC 计算时采用了静电求解模型，电场则通过差分求解泊松方程得到，对于二维轴对称模型，在柱坐标系中，泊松方程可以简化为

式中：为电势；为电荷密度；为真空介电常量；为径向距离；为纵向深度。

泊松方程的求解采用五点差分格式，使用超松弛迭代法，通过改变松弛因子以加速收敛。为了减少仿真计算量，采用虚拟阴极，位于推力器的中部。选择准中性模型作为虚拟阴极模型，在计算过程中的每个时间步长内，通过区域边界向等离子体发射一定量的电子，以维持阴极边界内的电子和离子满足准中性条件。计算阴极边界所有网格的总的净离子数Δ，具体的处理方式如下：

在推力器运行过程中，原子的速度约为百米每秒量级，小于离子和电子的运动速度，可以认为仿真中的原子为近似背景场。原子的分布情况满足通道等离子体扩散准则。由于二价以上离子含量极少，可以忽略，在模型中电子和原子碰撞只考虑生成一价离子，电子和原子碰撞的概率为

式中：为碰撞概率；为中性原子的数密度；为电子和中性原子的相对速度（可近似为电子速度）；为电子和中性原子的碰撞截面的总和，包括弹性碰撞截面、激发碰撞截面和电离碰撞截面，与电子能量相关；Δ为时间步长。

内部碰撞过程主要为电子和原子、离子之间的碰撞，忽略中性粒子相互之间的弹性碰撞。离子相比原子有着较高的速度和能量，加速喷出的离子即为推力器推力来源。粒子的运动和加速通过解耦计算获得，求解电场后将电场分配至粒子，得到粒子所受的电场力，通过静态磁场分布得到粒子所受的磁场力，进一步计算粒子的加速度，更新粒子的速度和位置。

和氙气相比，碘蒸汽在放电通道内部和阴极发射的电子相互作用，存在电离、解离、复合等多种不同的碰撞类型，产生的粒子种类更多，放电过程更为复杂。I和I的电离能分别为10.45 eV 和9.40 eV，Xe的电离能为12.10 eV。单个粒子的I和I的电离碰撞截面约为6.0×10cm和12.3×10cm，大于Xe 的4.8×10cm。I的解离能约为1.57 eV。在仿真中，主要实现了解离-电离过程，主要考虑的反应如下：

式中：e 为电子。

相关碘原子的弹性、激发、电离碰撞截面参数从文献［14-16］中获取。在模型中主要考虑的电荷交换碰撞（CEX 碰撞）反应如下：

上述CEX 碰撞截面表达式在文献［17］中已经给出：

式中：为CEX 碰撞截面；为相对碰撞速度；为碘的电离能，约10.45 eV；为氢的电离能，约13.60 eV；系数为1.81×10；系数为2.12×10。文献［7］中也采用了相同的公式进行模拟。

整体仿真流程可简述为：初始化计算模型，对计算区域进行网格划分，将磁场数据耦合加载到计算程序中，初始化物理参数，载入初始粒子。根据工质流量，在计算的每个时间步长内射入对应数目的中性碘分子，根据准中性条件射入对应数目的电子，依据粒子的现有速度计算其运动状态。通过PIC 方法建立粒子间的电磁场作用模型；通过DSMC 方法求解原子、离子间的碰撞和动量交换、电荷交换；通过MCC 方法求解电子和原子间的碰撞；通过求解泊松方程得到电势、电场分布，再根据加载的磁场数据计算粒子的加速度和速度。当流场中的粒子净流量小于设置值时，视为计算进行至推力器完全稳定工作状态，输出计算结果；如不稳定则继续进行计算，待计算至稳定状态后，再统计相关参数。

1.2 计算区域

推力器在周向不同方位角各类参数基本分布均匀，为了简化计算，仿真模型采用二维轴对称型结构。数值模拟时所选取的计算区域主要包括推力器的放电室、出口和部分羽流区域，左侧边界为推力器阳极边界，以推力器中轴线为对称轴，相关边界条件设置如图1 所示。碘原子和内外壁面之间的作用采用漫反射模型。根据推力器的尺寸，放电通道的长度=20.00 mm，内壁面处于=13.38 mm 处，外壁面处于=25.00 mm 处，在加速通道外部选取半径为38.00 mm、轴向长度为16.00 mm 的圆柱形羽流场区域。

图1 推力器结构及流场仿真区域Fig.1 Construction and simulated zone of the thruster

1.3 加速方法

碘离子和电子的质量相差约5 个数量级，在相同的能量条件下，电子速度将是碘离子速度的百倍以上，如果按照电子的时间步长同时推进碘离子，需花费大量时间和计算资源以达到稳定，采取降低碘离子质量和提高介电常数的办法来提高计算速度。具体处理方式如下：将碘离子质量降低到1/2 500，则其运动速度相应提高50 倍；将介电常数提高100 倍，则德拜长度相应提高10 倍，可以降低达到稳定时所需要的时间步，减少网格的总数量。

2 数值模拟以及分析对比

2.1 磁场分布

推力器磁场考虑为静态磁场，由ANSYS 商用软件计算得到分布结果，以输入文件形式将磁场数据加载至模拟计算程序。计算所得放电通道中轴线上磁感应强度分布如图2 所示，从阳极到羽流近场，整体呈先上升再下降趋势，最大值253.51×10T出现在通道出口附近，符合霍尔推力器内部磁场设计准则和分布规律。

图2 通道轴线磁感应强度分布Fig.2 Distribution of the magnetic flux density along the axis

2.2 设计工况下内部参数分布及性能参数

对设计工况下200 W 级推力器工作情况进行仿真，得到放电通道内的粒子数密度分布、电子温度和通道内离子轴向速度等结果。设计工况的输入参数包括：功率为200 W、放电电压为200 V、工质流量为13.83 sccm。初始电场通过设置左侧阳极边界为200 V、右侧自由边界为0 V 求解得到；稳态电场通过使用文献［19］中的方法，迭代求解方程组得到。离子数密度分布如图3 所示，最大值出现在通道中间区域，最大值约为4.12×10m，随着通道内电磁场对碘离子加速，数密度迅速降低，在通道出口位置的碘离子数密度约为7.66×10m。沿通道轴向离子数密度先增加后降低，分布规律和HOFER 等的结果相似，通道出口处碘离子数密度和文献［7］中的结果相似。推力器放电通道内分为近阳极区、电离区和加速区，满足放电通道粒子的分布特征。通道内的速度分布如图4 所示，离子速度在通道后半段迅速增加，并在电场作用下保持加速效果至通道出口，最大速度约为16 860 m/s。

图3 设计工况下离子数密度分布Fig.3 Distribution of ion number density under design

图4 设计工况下离子轴向速度分布Fig.4 Distribution of the ion axial velocity under the design condition

通道内碘原子数的密度先增大后减少分布情况如图5 所示，在内外壁面处积附了大量的碘原子，碘原子的数密度达到最高值；在电离区则由于大量碘原子在此处发生电离，碘原子的数密度较低。碘原子轴向的速度分布情况如图6 所示，离子被电场加速，在加速区速度较高，和中性原子发生碰撞后，区域内中性原子的速度也有所增加，达到千米每秒量级；其他区域内碘原子的速度量级则和文献［12］中一致。

图5 碘原子数密度分布Fig.5 Distribution of the number density of iodine atoms

图6 碘原子轴向速度分布Fig.6 Distribution of the axial velocity of iodine atoms

2.3 碘工质与氙工质对比

使用氙工质代替碘工质，在设计工况下进行数值模拟，比较2 种工质的性能。通道内氙离子的数密度分布如图7 所示，其分布趋势和碘离子大致相同，最大值出现在通道中央，约为1.38×10m，出口处最大轴向速度达到了约15 150 m/s。相较于碘工质，在通道内沿通道中央的氙原子数密度，从阳极开始就一直呈降低趋势，由于近阳极区电离率较低，下降趋势较平缓，在电离区则快速降低。对于碘工质，其沿通道中央的碘原子数密度呈先增大后降低的分布趋势，在阳极附近碘原子数密度较低，在轴向1～3 mm 处碘原子数密度明显增加，再在电离区快速降低，如图7 所示。说明碘作为双原子分子，其机理和氙原子存在差异。由于存在解离过程，并且解离阈值较低，在中性碘分子通过气体分配器进入到放电室后，需先和电子碰撞解离为碘原子，之后再进行下一步反应。碘分子数密度分布也证明了这一点，如图8 所示。在阳极附近处碘分子数密度较高，碘原子数密度较低，因为碘分子刚进入通道，还未充分解离；然后碘分子数密度快速降低，同处的碘原子数密度则明显上升，说明解离过程在此处频繁发生。对于碘工质霍尔推力器而言，除了传统的近阳极区、电离区和加速区外，还存在解离区，在此区域内，大量的碘分子解离成碘原子，其位置在近阳极区之后、电离区之前。在此计算条件下，其宽度略大于2 mm。为了提高推力器的性能，应设法提前解离区的位置，使碘原子充分电离，提高推力器效率。