多级会切场推力器振荡特性模拟研究

，，，

哈尔滨工业大学 先进动力研究所，哈尔滨 150001

多级会切场推力器振荡特性模拟研究

牛翔，伍环，刘辉*，于达仁

哈尔滨工业大学 先进动力研究所，哈尔滨 150001

目前多级会切场推力器(High Efficiency Multi-Plasma Thruster，HEMPT)数值模拟研究采用的主要手段是全粒子(Paticle In Cell，PIC)模拟方法，但是这种全粒子模拟对计算机的要求较高，因此文章通过建立HEMPT的一维流体模型，再现HEMPT的振荡现象，对放电电流和等离子体参数之间的关联性及等离子体参数分布规律进行研究。沿着两条典型的电子传导路径,利用MATLAB中Simulink模块建立一维流体模型。结果分析表明，推力器内部等离子体参数的振荡周期与放电电流振荡周期一致。推力器的主电离区在两个磁尖端之间，电势降主要集中在出口附近。

多级会切场推力器；一维流体模型；振荡；放电电流；电子传导；电势分布

与化学推进相比，电推进具有高比冲、小推力等特点。霍尔推力器是一种常见的电推力器，离子与壁面剧烈的碰撞是导致其寿命下降的重要原因，而会切场的磁场位型可以避免这一点，因而其具有长寿命的特点[1]。此外试验表明，会切场推力器可实现推力大范围连续可调[2]。它的工作原理如图1所示。

图1 多级会切场推力器结构原理Fig.1 Illustrative diagram of HEMPT

多级会切场推力器(High Efficiency Multi-Plasma Thruster，HEMPT)内部存在一系列的复杂物理过程，在研究这些物理过程中，数值模拟因其成本低，能够直接获得难以测量的参数等优点而被广泛应用。基于微观机理的全粒子(Paticle In Cell，PIC)模拟是研究推力器内部复杂物理过程的一种重要手段，文献[3]用PIC模拟了霍尔推力器与HEMPT推力器，并得到了HEMPT推力器内部电子与壁面碰撞少的结论。

PIC作为一种以粒子为对象的模拟方法，能够反映等离子体参数在空间的二维分布。文献[4]就利用PIC模型关于发射电子的径向分布对等离子体参数空间分布的影响进行了研究。但是目前PIC模拟存在对计算机性能要求较高、计算量较大等难点。与PIC相比，一维流体模型虽然不能反映等离子体参数的二维分布，但是算例的运行时间明显减少。另外，文献[5]在稳态等离子体推进器(Stationary Plasma Thruster，SPT)的研究中也验证了一维流体模型研究放电电流振荡特性的准确性，并且利用该模型得到的结论与胡俊锋等的试验结论相符[6]。因此在放电电流振荡模拟研究方面，使用一维流体模型可以保证在一定准确程度的情况下，显著降低算例运行时间。流体模型中电子采用稳态假设，可以很快得到推力器内部参数值，极大地减小了运算量。文献[7-10]利用流体模型对霍尔推力器内部呼吸模式振荡进行了模拟，文献[11-12]不仅利用流体模型对霍尔推力器的模式转换进行研究，还利用混合直接动力学模拟以获得更为准确的离子速度分布。

刘辉等利用PIC研究表明在有两个磁尖端的多级会切场推力器内部，电子有两条典型的电子路径[13]。当电子经过出口磁尖端进入通道内部后，在中部尖端区域电子路径开始分叉，有一部分电子会沿着磁感线运动至磁尖端，然后跨越磁感线继续运动，这类电子称之为捕获电子；还存在一部分电子并没有沿磁感线运动，而是直接通过通道中央从而到达阳极，这部分电子称之为漏失电子。漏失电子的产生是由于一部分电子在跨越第一个磁尖端后变成高能电子，在做拉莫尔回旋运动时，所受的离心力较大，因而能够摆脱后面磁尖端的束缚，从中央中轴线直接逃逸。本文沿着两条典型电子路径建立一维流体模型，在抓住电子主要的传导路径基础上简化了模型，减小了运算量。本模型不仅可以研究会切场推力器的低频振荡，还可以通过计算通道内部等离子体参数的轴向分布，为更精确的模型提供初始分布。但是，模型计算得到的等离子体参数分布是通道截面的平均值，因而无法获得在某一截面处的更为精确的等离子体参数分布，这是目前一维流体模型的不足之处。另外，由于一维流体模型主要用来研究会切场推力器的低频振荡，为了缩短收敛时间，这里没有用高维度的模型进行研究。

1 流体模型的建立过程

本文通过建立流体模型来对等离子体内部各个参数进行模拟。电子路径如图2所示，通过沿着两条电子路径建模可以将问题简化为一维。根据电子、原子、离子的运动过程作出如下假设：

1)为了避免双流体模拟带来的时间尺度问题，流体模拟采用对电子进行稳态假设，即忽略电子的时间偏导项。

2)电子在运动过程中，始终沿着磁感线运动。

3)模型仅考虑氙原子的一次电离碰撞，不考虑高价氙离子的产生、原子激发、库伦碰撞。

图2 电子路径Fig.2 Electron path diagram

对玻尔兹曼方程取速度的零阶矩，可以得到原子、离子的连续性方程；对玻尔兹曼方程取速度的一阶矩可以得到离子动量方程：

(3)

式中：n为等离子体体密度；na为原子密度；Va为原子速度，取常数Va=200 m/s；Vi为离子速度；mi为氙离子的质量，取值2.2×10-25kg；e为元电荷，取值为1.6×10-19C；E为电场强度；βi为电离率。

回路方程及电流方程如下：

式中：I为放电电流；l为通道长度，取值为0.1 m；A为通道面积，取值为5.3×10-4m2；Lc为通道电感，取值为10-6H；R为内阻，取值为10Ω；Vex为电子速度。

由于假设电子一直沿着磁感线运动，因此电子运动方程如下：

式中：me为电子质量，取值为9.1×10-31kg；Te为电子温度；B为磁感应强度；μ为比例系数，取值为1/3；νeff为有效电子碰撞频率：

式中：βm为电子与氙原子的动量碰撞交换系数；αB为波母系数，取值为1/32；ωce为电子自旋频率，取值(eB)/me。

由图2可知，通道内存在两条电子路径，为了将二维问题简化为一维问题，提出下面合并电子路径的公式：

式中：α为给定的参数，用来表示在中部电子通过中央路径引起的电流密度与总电流密度之比。jex1表示逃逸电子引起的电子电流密度，jex2表示捕获电子引起的电子电流密度。合并式(6)、(8)可得：

式中：νeff1为逃逸电子有效碰撞频率；νeff2为捕获电子有效碰撞频率。

由之前的电子处于稳态可忽略偏导项的能量方程：

(10)

式中：Δε为单次电离碰撞所引起的能量损失，由于只考虑一次碰撞，这里取值为12.1 eV。右边第一项为欧姆加热项，第二项为电离碰撞项，由于忽略电子碰壁及热传导造成的损失，因此在第二项前面加一个系数γi，使模拟结果更接近真实情况，这里取系数γi为3，文献[8]也做了相同处理。

边界条件处理如下：阳极位置的原子密度，离子密度和离子速度分别为2×1019m-3，6×1017m-3，2×103m/s。阴极处的电子温度为1 eV。模型的初始条件即初始原子和离子密度分布，离子速度分布分别如图3、图4所示。

图5为计算所用磁场。利用捕获电子与漏失电子数量之比α反映两种不同类型电子的相对大小。非尖端区不考虑波母传导即计算电子有效碰撞频率时不包括电子自旋频率的平方项。而在尖端区则考虑该项，尖端区与非尖端区的分界是以磁场强度导数为0定义的。

图3 初始原子、离子密度分布Fig.3 Initial atom and ion distribution

图4 初始离子速度分布Fig.4 Initial ion velocity distribution

图5 磁场分布Fig.5 Magnetic field distribution

2 计算结果及分析

在电压Us=300 V，α=0.1的条件下模拟推力器的振荡问题。如图6所示，放电电流产生剧烈振荡，时均放电电流约为7 A，电流峰值约为10.8 A，振幅约为5.8 A，振荡频率约为40 kHz。

图6 电流的振荡曲线Fig.6 Current oscillation

图7、图8分别为电子温度及归一化的方程源项的分布情况，其中点画线为磁分界面，虚线为电流振荡为峰值的时刻。方程源项是表征推力器内部电离程度的参数。源项越大，电离程度就越大。从图8中可以看出方程源项在时间上具有周期性，且周期与放电电流的振荡周期一致；在空间上，方程源项在归一化的轴向位置为0.4处达到最大，所以这个位置应为推力器的主电离区。

从时间上可以看出，当源项恰好处于极小值时，如在0.17 ms左右，振荡电流达到峰值，这是之前主电离区大量电离的结果。之后从图8中可以看出，电离程度变得很弱，因此电流也随之下降。在一个振荡周期里，主电离区内存在两个源项极大值，且其数值一小一大，也就是说电离会有两次不同程度增强。在第一次电离增强区域，电流随之增加，但是由于电离程度较小，所以电流所达到的极大值小于峰值;而第二次由于电离程度更大，导致第二次电流达到峰值，随之进入下一个周期。

从空间上看，主电离区的位置位于两个磁尖端之间，与高电子温度区域一致。在主电离区上游，由于电子温度较低，电离程度下降，所以源项较小。而在主电离区下游，一方面电子温度比较小，另一方面大部分原子已经被电离，所以在这一区域电离程度也很弱。

值得注意的是，电子温度分布也存在类似于源项分布的周期性，并且在一个周期内存在一大一小两个极值。

图9为瞬时电场分布。可以看出，在出口磁尖端附近电场强度总是处于极大值，这主要是由于出口磁尖端存在磁镜效应，电子跨越磁感线的运动很难实现，因而在出口处存在较大的电势降。除此之外，在通道中部至出口处存在随时间周期性变化的场强，而通过与电子温度的分布图7对比可以发现，在场强取得极大值处电子温度均较小，这可能是由于局部电子温度较低导致电离不充分而形成的电势降。

图7 电子温度分布Fig.7 Electron temperature distribution

图8 归一化方程源项Fig.8 Normalized source term

值得注意的是，在中部磁尖端之前存在着较小的电势降，而非磁尖端。而通过图2可以看出，中部磁尖端电离比较强，因而电势降提前可能是由于中部磁尖端附近电离增强造成的。

原子、离子密度分布如图10、11所示。可以看出，原子密度在归一化的轴向位置为0.4左右时会有振荡现象发生，而粒子密度分布从归一化轴向位置为0.4后也有明显振荡现象。从放电电流振荡图像中容易得到，多级会切场推力器的振荡频率在5 kHz左右，这与霍尔推力器内部低频振荡频率类似，研究表明霍尔推力器内部的振荡可能是由于原子供给与电离过程的不匹配造成的[14-17]。当局部电离增强时离子密度上升，原子密度下降，从而放电电流上升。此后离子经电场加速喷出，使离子密度下降，待原子供给充足后，电离继续增强，进入到下一轮循环。所以这也可能是多级会切场推力器内部低频振荡发生原因。

从图10和图11通道中部就可以看出原子密度下降的同时,离子密度上升，同时对应放电电流上升并且达到最大值。

图10 原子密度分布Fig.10 Atom density distribution

图11 离子密度分布Fig.11 Ion density distribution

3 结束语

本文沿着两条典型电子路径建立了HEMPT的一维流体模型，复现了放电电流振荡现象。模拟结果分析表明，HEMPT放电电流与等离子体参数振荡频率一致，推力器的主电离区在两个磁尖端之间，主要电势将集中在通道出口附近。利用建立的一维流体模型可以在保证一定精度的前提下，较快地得到推力器内部振荡参数及等离子体参数分布，为推力器振荡特性的研究提供便利。但是一维流体模型各个位置处均采用了截面平均的假设，这在一定程度上限制了模拟精确程度的提高，因此后续可以通过建立HEMPT的二维流体模型，在更高精度上对HEMPT内部振荡现象进行研究。

References)

[1] GILDEA S R，MATLOCK T S，MARTNEZSNCHEZ M，et al.Erosion measurements in a low-power cusped-field plasma thruster[J].Journal of Propulsion & Power，2013，29(4):906-918.

[2] GENOVESE A，LAZURENKO A，KOCH N，et al.Endurance testing of HEMPT-based ion propulsion modules for small GEO[C].32nd International Electric Propulsion Conference,Wiesbaden，Germany：IEPC，2011: 15.

[3] MATYASH K,SCHNEIDER R,MUTZKE A,et al.Kinetic simulations of SPT and HEMP thrusters including the near-field plume region[J].IEEE Transactions on Plasma Science,2009,38(9):2274-2280.

[4] BRANDT T，TROTTENBERG T，GROLL R，et al.Particle-in-cell simulation of the plasma properties and ion acceleration of a down-scaled HEMP-Thruster[C].30th International Spacecraft Propulsion Conference,Cologne,Germany：ISPC，2014.

[5] BOEUF J P，GARRIGUES L.Low frequency oscilla-tions in a stationary plasma thruster[J].Journal of Applied Physics，1998，84(7):3541-3554.

[6] 胡俊锋，刘辉，李建置，等.会切磁场推力器低频振荡特性[J].中国空间科学技术,2016，36(1)：26-34.

HU J F，LIU H，LI J Z，et al.Low frequency oscillation characteristics in a cusped field thruster[J].Chinese Space Science and Technology，2016，36(1):26-34(in Chinese).

[7] MOROZOV A I，SAVEL′EV V V.One-dimensional hybrid model of a stationary plasma thruster[J].Plasma Physics Reports，2000，26(10):875-880.

[8] MOROZOV A I，SAVEL′EV V V.One-dimensional hydrody-namic model of the atom and ion dynamics in a stationary plasma thruster[J].Plasma Physics Reports，2000，26(3):219-224.

[9] BARRAL S，AHEDO E.Theoretical study of the breathing mode in Hall thrusters[C]//42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference，Scacramento,California,July 9-12,2006.

[10] BARRAL S，AHEDO E.Low-frequency model of breathing oscillations in Hall discharges[J].Physical Review E,2009，79(4):046401.

[11] HARA K，BOYD I D，KOLOBOV V I.One-dimens-ional hybrid-direct kinetic simulation of the discharge plasma in a Hall thruster[J].Physics of Plasmas，2012，19(11):530-541.

[12] HARA K，SEKERAK M J，BOYD I D，et al.Mode transition of a Hall thruster discharge plasma[J].Journal of Applied Physics，2014，115(20):203304-2033010.

[13] LIU H，WU H，ZHAO Y，et al.Study of the electric field formation in a multi-cusped magnetic field[J].Physics of Plasmas，2014，21(9):214.

[14] MOROZOV A I，SAVEL′EV V V.One-dimensional hydrodynamic model of the atom and ion dynamics in a stationary plasma thruster[J].Plasma Physics Reports，2000，26(3): 219-224.

[15] MOROZOV A I，SAVEL′EV V V.One-dimensional hybrid model of a stationary plasma thruster[J].Plasma Physics Reports，2000，26(10): 875-880.

[16] FIFE J M.Hybrid-PIC modeling and electrostatic probe survey of Hall thrusters[D].Massachusetts: Massachusetts Institute of Technology，1998:189-195.

[17] BARRAL S，AHEDO E.Theoretical study of the breathing mode in Hall thrusters[C]//42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit,Scacramento,California,July 9-12,2006: 5172.

(编辑：车晓玲)

Simulationresearchonoscillationcharacteristicinhighefficiencymulti-plasmathruster

NIU Xiang，WU Huan，LIU Hui*，YU Daren

InstituteofAdvancedPower，HarbinInstituteofTechnology，Harbin150001，China

Now the main numerical simulation method of high efficiency multi plasma thruster (HEMPT) is particle in cell.But this method puts higher demands on computer.Therefore，this method aims at reappearing the oscillation phenomenon in high efficiency multi-plasma thruster by building up one dimension fluid model and studying the plasma parameters distribution and the relationship between the current and plasma parameters.A one dimension fluid model was built up along two distinguishing typical electron paths in thrusters by Simulink.The simulation results show that the plasma parameters oscillation period is consistent with that of the current.The main ionization zone is between the two cusped fields.The main area where potential drops is near the channel exit.

high efficiency multi-plasma thruster；one dimension fluid model；oscillation；current；electron conductivity；potential distribution

http://zgkj.cast.cn

10.16708/j.cnki.1000-758X.2017.0072

V430

A

2017-05-04；

2017-08-03；录用日期2017-09-12；< class="emphasis_bold">网络出版时间

时间：2017-09-24 16:01:06

http://kns.cnki.net/kcms/detail/11.1859.V.20170924.1601.008.html

国家自然科学基金(51421063，11505041，11275055)

牛翔(1995-)，男，硕士研究生，15776867067@163.com，研究方向为振荡与噪声分析

*通讯作者：刘辉(1981-)，男，副教授，huiliu@hit.edu.cn，研究方向为空间电推进

牛翔，伍环，刘辉，等.多级会切场推力器振荡特性模拟研究[J].中国空间科学技术，2017，37(5)：75-80.NIUX，WUH，LIUH，etal.Simulationresearchonoscillationcharacteristicinhighefficiencymulti-plasmathruster[J].ChineseSpaceScienceandTechnology，2017，37(5)：75-80(inChinese).