Investigation on the Numerical Design of the Beam Current Production Process of Ion Thruster
-
摘要:
离子推力器具有比冲高、寿命长、推力连续可调等优点，已成为未来空间探测任务推进系统的首选。离子推力器产品的优化设计通常采用实验测量的手段来实现，现有的仿真软件均是针对离子推力器的各个组成组件来开发，忽略了组件之间的相互耦合，为此，亟待开发一款综合考虑各组成部件相互作用的离子推力器工作全过程数值仿真软件。离子推力器结构复杂，且涉及等离子体场、流场、电场、粒子间碰撞以及粒子与壁面相互作用等过程，现有的仿真手段很难实现在同一时空下的数值模拟。提出了一种近似耦合的数值仿真方法，即仍将离子推力器的各组成部件作为一个独立的个体来仿真，但在数值建模时，考虑相关联部组件的影响。分别构建了两两耦合的数值仿真模型，开展了数值仿真计算，仿真结果显示部组件单独工作和耦合工作时组件内部的等离子体微观参数将发生较大变化。
Abstract:Ion thruster has become the first choice of the propulsion system for the future space exploration missions for its high specific impulse, long lifetime and continuous thrust adjustment. Optimization of the ion thruster is usually realized by large numbers of experimental tests. Current simulation models developed only focus on the components themselves and ignore the interactions between components. The professionally numerical model that takes these interactions relationship into account should be developed. However, the complicated structures and coupled interactions between plasma, electric potential, magnetic field, particle collisions, particles and walls, it isimpossible to simulate the working processes simultaneously. A similar coupling simulation method is proposed, where the component itself is treated as a whole but it is simulated by giving the coupling interaction as the boundary conditions. Coupling simulation models of two components are established independently and their interactions are simulated. The results show that the plasma micro-parameter changes tremendously when it works alone or not.
-
Keywords:
- ion thruster /
- working process /
- simulation
-
0. 引言
与化学推进相比，离子推力器具有推力连续可调的特点，意味着完成同样的轨道控制任务，离子推力器的控制精度要远远高于化学推进，同时，由于离子推力器高比冲的特点，使其在执行航天任务时所消耗的推进剂量将显著低于化学推进，航天器的有效载荷大幅提升[1-3]。
离子推力器[4-6]是一种静电型电推进类型，电子和中性气体碰撞在放电腔体内产生高浓度等离子体，等离子体中的离子再经过强电场加速喷出，产生卫星所需的推力。相比其他类型，离子推力器具有比冲高、寿命长、推力高、分辨率连续可调等显著技术优势[7-10]，可以大幅提升卫星的有效载荷，因而自1960年美国科学家戈达德首次提出离子电推进概念至今，已有上百颗离子推力器产品用于执行卫星的位置保持、轨道转移和轨道提升等任务。随着未来空间探测任务的逐年增多，离子推力器产品的应用需求越来越多，推力器产品的性能和工作寿命要求也越来越高。
针对离子推力器的性能和工作寿命，国内外已开展了相关的理论分析[11]、数值仿真[12]和实验测量[13]研究。总体来看，不管是哪一方面的研究，其主要关注点均为部组件，即通过优化部组件的设计，来提升离子推力器单机产品的性能和寿命。对于数值仿真，现有的仿真软件缺乏组成部件之间的相互耦合，因此需在软件开发过程中对部组件之间的相互耦合关系进行解耦处理。事实上，这些组成部组件在离子推力器的工作过程中是相互作用、相互影响且缺一不可的。而现有的仿真软件却无法很好地研究部组件耦合作用下的物理过程及其影响，因而急需开发一套集成化的数值仿真软件，以期通过仿真手段，再现离子推力器工作过程，旨在从等离子体的微观角度，描述离子推力器工作时内部等离子体在等离子体流场、热场、电场、磁场等多物理场耦合下的产生及演化过程，快速分析影响产品性能和工作寿命的关键影响因素，从理论角度提出产品的优化设计方向和解决措施。
本文将针对这一问题，开展离子推力器工作全过程的数值仿真方案研究。
1. 离子推力器结构组成和基本原理
按照电推进的工作原理，将电推进分为电热式、静电式和电磁式三种，三种电推进类型的主要区别在于等离子体的产生方式。离子推力器是一种静电式电推进类型，通过电子和中性气体电离碰撞产生等离子体，等离子体中的离子在栅极组件静电场作用下加速引出，形成束流，产生推力。
离子电推进系统包括离子推力器、电源处理单元、推进剂贮供单元和控制单元等，其中离子推力器是离子电推进系统的核心单机，其他单机则负责离子推力器的供电、供气、工作时序的控制。传统离子推力器由空心阴极、放电室、栅极组件和中和器组成，其结构组成如图1所示。从组成结构来看，放电室的前端是空心阴极，而放电室又是栅极组件的输入，栅极组件的输出则是羽流区束流形成的输入，羽流区形成束流的大小又决定着空心阴极的发射电流。因此，四个组成部件之间相互作用、相互影响。
空心阴极是离子推力器的电子源，向放电室内部提供电离中性气体所用的电子，电子在放电室内部电磁场作用下沿着磁力线做螺旋加速运动，并和其中的中性气体碰撞，使其电离产生等离子体，等离子体中的离子在扩散运动和电场作用下向栅极组件方向运动，大部分离子穿过屏栅孔进入栅极组件，经过屏栅和加速栅强电场的作用加速聚焦引出，和中和器发射的电子中和，形成准中性的等离子体。
2. 离子推力器工作全过程集成化数值仿真
首先从离子推力器的性能参数和各组成部组件的几何结构和工作参数之间的变化关系来进行分析。由离子推力器的性能参数计算公式可知，只要确定离子推力器工作参数，即可获得其性能参数。这些工作参数为空心阴极、放电室、栅极组件和中和器的电压、电流，具体如图2所示。
从目前的实验情况来看，对于离子推力器产品，开展不同阶段，即性能调试、力学测试、热真空实验、PPU联试和长寿命实验，同时利用实验测量手段对实验测试过程中四个部组件的工作参数进行测量，比如性能调试实验过程中空心阴极的阴极电压和电流，放电室的阳极电压和阳极电流，栅极组件的屏栅电压、电流，加速栅电压、电流，中和器的电压和电流，说明了推力器的整个过程是四个部组件共同作用完成，每个部组件都有其各自的作用。
由部组件的性能计算公式可知，离子推力器四个组成部件的性能参数同样也与其各自的几何结构参数和工作参数有关。当离子推力器的几何结构参数恒定时，即可获得离子推力器以及各个部组件的性能参数。
图3给出了四个组成部件的输入输出参数之间的关系。由图3可知：连接空心阴极和放电室之间的参数为电子的微观特征参数；连接放电室和栅极组件的参数为等离子体的微观特征参数，如离子速度、密度、位置和能量；连接栅极组件和羽流的参数为离子的微观特征参数，如离子速度和密度。
![]()
图 3 离子推力器四个组成部件工作参数之间的耦合关系Figure 3. Coupling relationships between operational parameters of an ion thruster对于空心阴极，其输入工作参数为阴极触持极电压，输出参数为电子微观特征参数和放电室的输入参数。放电室其他输入参数为阳极电压、屏栅电压，输出参数为等离子体微观特征参数和栅极组件的输入参数。栅极组件的其他输入参数为屏栅电压、加速栅电压，输出参数为离子微观参数和羽流输入参数。羽流的输入参数为中和器电位，输出参数为束流。由此，对离子推力器产品而言，要完成单机工作全过程的数值仿真，在确定产品几何结构参数和供气参数的情况下，仅需要输入各部件工作参数即可。
但实际上，这一过程难以实现。第一，离子推力器的每个组成部件的几何结构组成和工作过程非常复杂，而且每个组成部件所关注的点或者带电粒子有所不同，比如空心阴极内部电子为主体，放电室内部离子和电子均为主体，栅极组件内部的主体为离子，羽流区的主体又为离子和电子，因此根据关注点的不同，选择的数值仿真手段也不同，现有的仿真手段还无法做到多种仿真方法在不同工作区域的随意切换，比如电子在放电室中通常采用粒子云方法，通过跟踪每个电子的运动轨迹得到它在气体放电瞬间的动力学行为，而在栅极组件中电子仅存在于屏栅上游和加速栅的下游，在没有电子反流发生的情况下，它的存在几乎不会影响栅极组件中离子的聚焦和引出过程，因而在栅极组件的数值计算时通常采用流体方法对其进行模拟。在放电室和栅极组件之间的离子直接关系着放电室的性能，以及栅极组件的引出性能和工作寿命，为了更好地得到离子的瞬时特征参数，通常采用粒子云的方法对其进行模拟。不同的带电粒子在不同的计算区域都有其分别对应的适用数值仿真方法。第二，各个组成部件的几何尺寸相差非常大，比如栅极组件的双栅间距仅为放电室轴向尺寸的千分之一，因而穿过屏栅孔的离子将高速穿过栅极组件，若要通过仿真手段模拟这一过程，就需要对栅极组件的计算区域进行网格精细化处理，上千甚至上万个栅极孔的束流引出模拟对计算机的计算量来说也是一种挑战。因此，从数值仿真的手段和计算机的计算能力角度来看，离子推力器工作全过程的集成化仿真很难实现。
本文提出了一种近似耦合的数值仿真方案，即各个组成部件仍然作为独立的个体进行数值仿真，但在仿真时考虑关联部组件的影响，关注部组件之间耦合对各个部组件性能或工作寿命的影响。图4给出了近似耦合数值仿真方案示意图。
下面将详细介绍各个部组件之间的耦合仿真。
2.1 空心阴极和放电室耦合数值仿真
对于空心阴极和放电室的耦合问题，在建立空心阴极数值仿真模型时，扩大了阴极下游放电室内部的计算区域，以考虑放电室阳极电压对阴极出口处电子运动行为的影响，具体的仿真模型如图5所示。模型的计算边界包括上、下、左和右边界，其中上边界、下边界和左边界设置为开放边界，右边界为定值边界，代表阳极，下边界也称对称边界。在阴极点火计算过程中，发射体温度设为恒定值，为
1860 K，点火过程计算至稳态，温度设为浮动，之后再进行热稳态计算。在定值边界上，由于$ \varphi (z,r + \Delta r) = \varphi (z,r) $，如果使用中心差分计算得到的电场值，就会比真实值小大约一半，作为一种简易的近似，这里直接用边界下方相邻节点的电场值作为边界上的电场值。![]()
图 5 空心阴极和放电室耦合仿真模型Figure 5. Coupling simulation model of the hollow cathode and discharge chamber采用粒子云方法通过物理场实现粒子间作用，将粒子携带的物理量权重代入计算网格中得到各类物理场，再差值得到粒子所在位置的场值，由此对粒子进行相应处理，当右边界吸收的离子和电子个数不再随时间发生变化时，模型达到了稳态。此时可获得阴极内部以及阴极出口处的电子微观特征参数。
图6所示分别为空心阴极内部的电子能量和速度分布。从计算结果来看，当给定空心阴极几何结构参数和放电室阳极电压时，可以模拟得到空心阴极内部的电子微观特征参数。
![]()
图 6 空心阴极内部电子能量和速度分布计算结果Figure 6. Eimulation results of electron energy and velocity distribution in the hollow cathode之后，采用该仿真计算模型，进一步扩大放电室的仿真计算区域，在阴极的下游距离25 cm的位置处选择一片阳极，给定阳极电压，模拟空心阴极单独放电情况，而在阴极出口给定与实际产品结构相同的阳极结构，同样给定阳极电压，模拟空心阴极和放电室耦合放电的情况。
图7和图8分别给出了空心阴极单独放电和耦合放电时的空间电势分布和电子密度分布。结果显示，独立放电时空间电势分布较为平坦，变化均匀，梯度较小，特别是在远离边界的中间部分，其等势线近似平行于平行阳极板边界，而耦合放电电势向放电室径向方向扩散。而电子密度的等势线也说明了耦合放电时电子离开空心阴极后运动轨迹呈发散状，沿着轴向和径向均有位移。
图9所示为空心阴极和放电室耦合作用下触持极表面的电场流场分布和离子运动轨迹。结果显示，空心阴极下游对应的放电室轴线区域的电势最低，使得整个放电室内的电场流场呈现全部指向阴极壁面的趋势，阴极下游的电场基本由垂直电场组成。在这些电场的作用下，位于阴极下游的离子直接撞击到阴极端面上。由此可见，当空心阴极和放电室耦合工作时，放电室内的气体放电过程将极大地影响阴极下游的电子运动行为及其特征参数，进而影响到放电室内部等离子体的整体运动行为。利用该模型可以实现空心阴极单独工作、空心阴极和放电室耦合工作过程的数值模拟。
图10所示为仿真得到的触持极表面的溅射腐蚀形貌。由结果可以看出，触持极表面出现了凹凸不平的刻蚀痕迹，这进一步证实在图9（b）所示的离子作用下，触持极表面会受到大量离子的轰击，这会影响空心阴极的工作寿命。
2.2 放电室和栅极组件耦合数值仿真
在放电室和栅极组件的耦合数值仿真方面，本文构建了两个数值仿真模型，即不考虑屏栅和考虑屏栅，具体如图11（a）和图11（b）中红色实线和红色虚线所示。不考虑屏栅，指的是在放电室数值仿真建模时，假设屏栅表面仅为一个边界条件，碰撞到它表面的离子均被其吸收，电子被反射，此时不会考虑从屏栅孔泄漏出去的离子个数。而考虑屏栅的方案是指放电室仿真计算时，虽然将屏栅表面作为了仿真计算的右边界，但同时考虑从该边界泄漏出去的离子和原子个数，这里离子和原子的个数根据屏栅表面的离子透明度和原子透明度来确定。模型的输入条件是图8（b）所示阴极出口处的电子速度分布，计算边界为空心阴极、阳极和屏栅，其中阴极电位为触持极电压，阳极电位为实验测到的阳极电压，屏栅电位为屏栅电压。放电室内部的磁场分布由磁体系统产生，采用商用软件对其进行计算并读入。计算采用粒子云方法和蒙特卡罗方法相结合的方式，通过对带电粒子运动轨迹的跟踪和概率算法对粒子间碰撞过程进行处理，得到带电粒子的微观特征参数。整个计算过程中性原子假设为均匀分布，它和原初电子、二次电子发生弹性碰撞、激发碰撞和电离碰撞，被电离和穿过屏栅孔的原子均做删除处理。根据删除的原子个数对放电室内部的中性原子密度进行实时更新，当右边界统计的离子电流达到设计值时，计算达到稳态。
图12为仿真计算得到的两种情况下的放电室内部电势分布。结果显示，当放电室内部考虑屏栅和不考虑屏栅时，其放电室内部的电势分布发生变化，相比于不考虑屏栅的情况，考虑屏栅时，其电势最小值从原来的206.1减小至203.1，而且靠近屏栅表面的电势从原来的206.5增大至207.6。从推力器出口即z方向来看，考虑屏栅后，该方向的电势变化梯度增大，这点从放电室的径向方向电势变化也能看出，径向方向的电势梯度增大，这意味着将有大量的离子将沿着电势梯度变化大的方向运动。
![]()
图 12 放电室内部电势分布Figure 12. Electric potential distribution inside of the discharge chamber图13为两种情况下仿真计算得到的放电室内部等离子体密度分布。结果显示，考虑屏栅后，放电室内部的等离子体沿着放电室的径向和轴向方向扩散，使得腔体内的等离子体向周围运动，而充满整个空间，放电区域变大。屏栅上游的等离子体沿着放电室的径向方向范围扩大，但栅极中心区域的等离子体密度没有明显增大。
图14为仿真计算得到的两种磁拓扑结构对应的放电室内部等离子体密度在径向和轴向方向的分布。该结果进一步说明了在放电室仿真计算时，是否考虑屏栅的离子透明度对放电室内部的等离子体特征参数及其分布有很大影响，这就导致两种情况下进入栅极组件内部的离子密度有所不同。但从仿真情况来看，考虑屏栅的离子透明度，不会改变栅极中心区域的等离子体密度，但会扩大屏栅上游的等离子体密度分布区间，使得栅极最边缘将有离子穿过栅孔。
![]()
图 14 放电室等离子体在径向和轴向方向的分布Figure 14. Plasma distribution in the radial and axial directions将图14（a）和图14（b）所示放电室内部且距离屏栅上表面一个德拜长度的等离子体密度分布结果作为栅极组件的仿真计算输入条件。考虑到上千甚至上万个栅极孔的计算量，本文选择一个栅极孔作为仿真计算的研究对象，因栅极组件中心区域离子极易出现欠聚焦而边缘区域则发生过聚焦现象，因而在中心区域和边缘区域分别选择一个栅极孔，其上游放电室内的等离子体密度则为图14计算结果对应的中心区域和边缘区域的计算结果。图15（a）和图15（b）分别给出了在栅极组件中心区域，放电室内部等离子体密度分别为2.2×1017 m−3和8.0×1016 m−3时的三栅极组件内等离子体密度分布。计算过程中采用了粒子云方法跟踪了离子在电场作用下的加速过程，运用流体方法来处理屏栅上游等离子体中的电子和栅极组件下游中和器发射电子的运动行为，粒子间的碰撞类型考虑为离子和中性原子之间的弹性碰撞和电荷交换碰撞。模型假设中性原子呈均匀分布。当计算区域所有节点电场变化小于0.05%时认为程序达到收敛。图15的计算结果显示，当屏栅上游的等离子体密度从8.0×1016 m−3增大至2.2×1017 m−3时，束流离子的聚焦状态发生了改变，但是两种情况都属于正常聚焦，没有任何离子被加速栅或减速栅所截获。
![]()
图 15 栅极组件中心区域的束流离子密度分布Figure 15. Beam ions density distribution at the center region of ion optics对于栅极组件的最边缘栅孔，最担心的是由于其上游等离子体密度较低，离子很容易发生过聚焦，因此本文也模拟了屏栅上游等离子体密度较低时的束流离子密度分布。图16为屏栅上游等离子体密度为1.0×1015 m−3时的束流离子密度分布。结果显示，束流离子聚焦正常，且无任何离子被加速栅或减速栅所截获。
![]()
图 16 栅极组件边缘区域的束流离子密度分布Figure 16. Beam ions density distribution at the edge region of ion optics2.3 栅极组件和羽流耦合数值仿真
对于栅极组件和羽流耦合的数值仿真，这部分通常的处理是将栅极组件作为一个整体，但是事实上，栅极组件由上千甚至上万个栅极孔组成，如果同时考虑上千甚至上万个孔的数值仿真，这不管是对数值仿真手段还是对计算机的能力都是一种巨大的挑战。而且栅极不同的区域因其屏栅上游等离子体密度的差异导致束流离子在栅极组件内的密度分布区别较大，因而通常在栅极组件数值仿真时仅考虑单个栅极孔或两个栅极孔。但因羽流区束流离子和电子的中和过程为离子束和电子束的整体行为，因而在构建栅极组件和羽流耦合仿真模型时，仍将栅极组件和中和器分别作为独立的整体来处理。图17为仿真计算模型。图中蓝色区域为中和器，r方向红色实线所示的区域为栅极出口，其他区域为羽流，zmax和$ {r_{\max }} $分别为计算区域轴向和径向最大值。图17给出了栅极组件和羽流耦合的二维数值仿真模型，根据推力器的结构组成，中和器沿推力器的轴向方向（zmax方向）离推力器出口面有一定的距离。因而在中和器建模时在横向和径向方向均给了一定的长度，见图17中蓝色实线所示区域，栅极组件则为图17中红色实线所示部分。离子从栅极组件区域进入计算区域，每个时间步长进入的离子个数由推力器的束流来决定。电子则从中和器出口（图17蓝色实线部分的红色区域）喷出，每个时间步长进入计算区域的电子个数由中和器的发射电流来获得。
![]()
图 17 栅极组件和羽流耦合仿真计算模型Figure 17. Coupling simulation model of the ion optics and the plume计算流程是每个时间步长离子和电子分别从栅极组件出口和中和器出口喷出进入计算区域内部，采用电荷权重法分别将离子和电子所带的电量权重代到网格节点上，得到该节点上的离子密度和电子密度，将其代入泊松方程，得到计算区域内的电势分布。根据牛顿第二定律分别对离子和电子进行加速，并根据速度结果更新粒子位置。采用蒙特卡罗方法处理电子和离子之间的弹性碰撞和中和反应过程。利用边界条件对超出计算边界的带电粒子进行边界处理。重复以上步骤，直到计算区域所有节点的电场变化小于0.05%时认为程序达到了稳定状态。
图18给出了仿真计算得到的中和器电位为0.1 V，栅极组件内等离子体密度分别为2.37×1017 m−3和1.19×1017 m−3时，计算区域内的电势分布。结果显示，当栅极组件内的等离子体密度减小时，羽流区对应的电势也随之降低，距离栅极出口下游10 mm的位置处，其电势由原来的35 V减小为20 V，而且此时区域内的电势线向栅极组件方向靠拢。
图19为仿真计算得到的两种密度下的束流离子密度分布。结果显示，羽流区电势的减小并不会影响该区域束流离子的密度分布，分析认为这是由于从栅极口喷出的离子运动速度过快导致的。
3. 总结
本文采用近似拟合的方法对空心阴极、放电室、栅极组件和中和器这四大部组件的两两组件耦合进行了仿真建模和数值仿真，从计算结果来看，模型计算结果能够一定程度上反映组件耦合对等离子体微观参数的影响，由此可以进一步获得部组件的性能参数变化。该方法可以有效地分析相邻两个部组件之间的相互作用及其影响，结果能够为未来离子推力器产品或部组件优化设计以及产品研制提供数据支持和技术支撑。
-
![]()
图 3 离子推力器四个组成部件工作参数之间的耦合关系
Figure 3. Coupling relationships between operational parameters of an ion thruster
![]()
图 5 空心阴极和放电室耦合仿真模型
Figure 5. Coupling simulation model of the hollow cathode and discharge chamber
![]()
图 6 空心阴极内部电子能量和速度分布计算结果
Figure 6. Eimulation results of electron energy and velocity distribution in the hollow cathode
![]()
图 12 放电室内部电势分布
Figure 12. Electric potential distribution inside of the discharge chamber
![]()
图 14 放电室等离子体在径向和轴向方向的分布
Figure 14. Plasma distribution in the radial and axial directions
![]()
图 15 栅极组件中心区域的束流离子密度分布
Figure 15. Beam ions density distribution at the center region of ion optics
![]()
图 16 栅极组件边缘区域的束流离子密度分布
Figure 16. Beam ions density distribution at the edge region of ion optics
![]()
图 17 栅极组件和羽流耦合仿真计算模型
Figure 17. Coupling simulation model of the ion optics and the plume
-
[1] CHOUEIRI E Y. A Critical history of electric propulsion:The first 50 years (1906-1956)[J]. Journal of Propulsion and Power,2004,20(2):193−203. doi: 10.2514/1.9245
[2] JORNS B,MIKELLIDES I,MAZOUFFRE S,et al. Physics of electric propulsion[J]. Journal of Applied Physics,2022,132(11):110401. doi: 10.1063/5.0118076
[3] JAHN R G. Physics of electric propulsion [M]. New York:McGraw-Hill,1968.
[4] STUHLINGER E. Ion propulsion for space flight [M]. New York:McGraw-Hill,1964.
[5] BREWER G R. Ion propulsion technology and applications [M]. New York:Gordon and Breach,1970.
[6] KAUFMAN H R. Technology of electron-bombardment ion thrusters[J]. Advances in Electronics and Electron Physics,1975,36:265−373.
[7] SHIMADA S,SATO K,TAKEGAHARA H. 20 mN class xenon ion thruster for ETS-VI [C]//19th International Electric Propulsion Conference,AIAA-1987-1029,Colorado:Colorado Springs,1987.
[8] BEATTIE J R. XIPS keeps satellites on track[J]. The Industrial Physicist,1998,4(2):24−26.
[9] BEATTIE J R,MATOSSIAN J N,ROBSON R R. Status of xenon ion propulsion technology[J]. Journal of Propulsion and Power,1990,6(2):145−150. doi: 10.2514/3.23236
[10] GOEBEL D M,KATZ I. Fundamentals of electric propulsion:ion and hall thruster [M]. California Institute of Technology,Jet Propulsion Laboratory,2008.
[11] BROPHY J R,WILBUR P J. Ion thruster performance model [R]. Colorado State of University,Department of Mechanical Engineering,NASA CR-174810,1984.
[12] MAHALINGAM S. Particle based plasma simulation for an ion engine discharge chamber [D]. Ohio:Wright State University,2002.
[13] BROPHY J R. NASA's deep space 1 ion engine[J]. Review Scientific Instruments,2002,73(2):1071−1078. doi: 10.1063/1.1432470
下载:
计量
- 文章访问数: 17
- HTML全文浏览量: 3
- PDF下载量: 5
