自身场磁等离子体推力器数值仿真流场特性分析*

操乐晖，吴 鹏，罗梓浩，任军学，汤海滨,3,4

（1.北京航空航天大学 空间与环境学院·北京·102206；2.北京航空航天大学 宇航学院·北京·102206；3.航天器设计优化与动态模拟技术教育部重点实验室·北京·102206；4.空间环境监测与信息处理工业和信息化部重点实验室·北京·102206）

0 引 言

磁等离子体推力器（Magnetoplasmadynamic Thruster，MPDT）利用高电流对推进剂进行电离，并产生较高电离度的等离子体，再通过电磁场与等离子体相互作用产生洛伦兹力，将等离子体加速喷出。从目前已有的关于MPDT的研究来看，MPDT比冲的量级最高可达1万s以上，相较于传统的化学推进要高出10～20倍；此外，推力量级最高能达到100多N，远大于静电式推力器和电热式推力器[1]。MPDT的性能特点使得该推力器在未来严峻的深空探测任务中将扮演重要角色，由于MPDT需要的电功率偏大，结合未来空间核电技术，发展大功率的核电磁等离子体推力器将大大提高航天器的性能，减少发射成本，以及可进行更长时间更远距离的飞行任务[2]。

MPDT根据磁场的来源不同，可分为自身场磁等离子体推力器（Self-Field Magnetoplasmady-namic Thruster，SF-MPDT）和附加场磁等离子体推力器（Applied-Field Magnetoplasmadynamic Thruster，AF-MPDT）。SF-MPDT的磁场来源主要由施加在阴阳极之间的电流产生，该电流会在推力器内部及羽流处感应出角向磁场，该角向磁场与等离子体中的径向电流作用产生沿轴向的洛伦兹力jrBθ。通常SF-MPDT需要比较高的电流才能在推力器内部产生比较高的磁场，其电流需要在几千A甚至是十几千A才能获得比较好的性能参数，因此SF-MPDT的功率会比较高。AF-MPDT的磁场来源于外加线圈和永磁体，并在推力器内部及羽流处产生轴向和径向的磁场，该磁场与等离子体相互作用，产生更加复杂的推力产生机制以及等离子体流动行为。一般情况下，通电线圈或永磁体在推力器内部产生零点几T的磁场，此外阴阳极之间施加的电流较自身场要小很多，其量级在几百A左右。

对于SF-MPDT，其加速机制比较简单，主要有气动加速和电磁加速，阴阳极之间的高电流将推进剂加热并电离，形成高温高密度的等离子体，通过喷管碰撞喷出，该方式为气动加速；电磁加速主要来源于电流与磁场的相互作用，根据广义欧姆定律（式（1））可将电磁加速机制归结为电场效应（等式右边第一项）、流动效应（等式右边第二项）、电子压力效应（等式右边第三项）和霍尔效应（等式右边第四项）四种作用。

（1）

其中，J为等离子体电流；E为等离子体内部电场；ne为电子数密度；B为磁感应强度；σ为电导率；e为电子电荷量；pe为电子压力梯度；u为流体速度。

附加场MPDT由于存在径向磁场和轴向磁场，其推力产生机制有如下几种：

1）气动加速。其机制与SF-MPDT气动加速类似。

2）自身场加速。由于AF-MPDT也需要在阴阳极之间施加电流，该电流也会在推力器内部感应出角向磁场，该角向磁场与电流相互作用产生轴向的洛伦兹力。但通常AF-MPDT的电流偏小，感应出的角向磁场偏小，与附加场相比几乎可以忽略，该加速机制在总体加速机制中占比相对较小。

3）涡旋加速。在附加磁场和电流的相互作用下，产生周向洛伦兹力jrBz和jzBr，会带动等离子体产生周向的运动。通过具有磁镜效应的磁喷管，将高速旋转等离子体的旋转动能转化为轴向动能。实验中观察到这一能量转化机制通常被认为是附加场的主要加速机制。

4）霍尔加速。在强磁场和低质量流量的工况条件下，根据广义欧姆定律，角向会感应出比较强的电流，该电流与附加磁场相互作用，产生的轴向洛伦兹力jθBr作用在等离子体上。对于霍尔加速，目前存在比较多的争议，主要是该机制产生的力对推力的贡献方向并不是很明确[3]。

目前已有理论对SF-MPD内部等离子体参数进行预测[4]，此外在实验上建立了国内大功率的附加场MPDT设计以及推力测量手段[5-6]，并在此基础上开展了推力器内部物理机制的研究[7]。但实验中由于推力器内部等离子体的恶劣环境，导致测量方面存在困难，并不能很好地解释推力器内等离子体产生的物理现象，如onset现象等[8]。因此，需要在数值模拟的基础上，对推力器内部的等离子体进行描述。目前，描述等离子体行为的方法主要有：粒子（Particlein Cell，PIC）方法以及磁流体（Magnetohydrodynamics，MHD）方法。PIC方法通过追踪推力器内的微观粒子，研究其在电磁场中的运动情况[9-11]。由于MPDT中数密度相对较高，该方法将会非常耗时。MHD方法假设等离子体为导电流体，通过求解流体方程得到等离子体的流场特性，虽然不能描述等离子体的微观特性，但对于宏观特性，则能够很好地进行预测。

本文使用MHD方法，对自身场磁等离子体推力器内部及羽流进行建模，利用高精度数值方法求解流场和磁场，并进行耦合求解。

1 物理模型

1.1 流体模型

当粒子的速度分布达到麦克斯韦速度分布的时间远小于流场的特征时间时，可以将等离子体作为流体进行处理，描述该近似可以用Knudsen数表示[12]，即

其中，λ为分子平均自由程;L为特征长度。

碰撞截面与碰撞类型以及温度有关，且推力器的几何构型是一定的，因此其大小主要与推力器中的等离子体数密度有关。自身场推力器的数密度能达到1021～1022/m3,且由于感应磁场垂直于流场，使得碰撞更加激烈，缩短了达到麦克斯韦速度分布的时间，更加符合流体近似。此外，波与粒子相互作用产生的不稳定性,导致反常碰撞进一步满足流体近似。

根据自身场MPDT内部等离子体流动情况，使用以下假设条件建立物理模型：

1）内部流动为定常、轴对称、层流、无黏流动；

2）准中性假设，不存在电荷分离的现象；

3）热力学非平衡假设，即双温度模型；

4）忽略辐射能量损失；

5）认为电导率为标量，并考虑反常电导率；

6）考虑霍尔效应，但忽略离子滑移。

对于流场的描述，使用以下控制方程：

连续性方程

（2）

动量方程

（3）

总能量方程

（4）

电磁场计算方法由麦克斯韦方程组导出。麦克斯韦方程组微分形式如下

（5）

其中，ε0为真空介电常数；μ0为真空磁导率。

由于准中性假设，认为无电荷分离。此外，对于绝大部分问题，等离子体频率通常高达几十MHz，位移电流比传导电流小很多，可以忽略不计，因此可以得到简化的麦克斯韦方程组如下

（6）

并结合广义欧姆定律（式（1）），消去电场项，得到磁感应方程如下

（7）

其中，σ为等离子体电导率，由于自身场磁场比较单一，该电导率可处理为标量。

流体与磁场的耦合主要通过洛伦兹力和焦耳加热联系起来，因此可以通过求解流场参数分布求解磁场。即通过磁感应方程求解磁感应强度B，再由安培环路定律求出电流密度J，并由广义欧姆定律得到电场强度E，代入流场的源项中，从而使得磁场与流场联系起来。

但以上方程组并不封闭，需要引入气体状态方程

p=nekTe+nikTi

（8）

其中，ne、ni分别为电子和重粒子数密度；k为玻耳兹曼常数；Te、Ti分别为电子温度和重粒子温度。

1.2 热力学非平衡模型

要想维持等离子体内部热平衡态，需要碰撞的时间远小于等离子体流动特征时间，才能使得等离子体维持麦克斯韦速度分布。在MPDT中，数密度的数量级只有o（21）～o（22），而等离子体宏观流速一般是超声速的，经计算MPDT内部等离子体的停留时间与达到平衡的时间相当[13]，不频繁的碰撞使得电子和离子并不能达到局部热平衡，因此需要将离子和电子分开考虑。建立非平衡模型，即双温度模型。此外，实验结果也观察到了这一非平衡态[14]。

假设电磁场通过对电子的作用对电子进行供能，电子通过碰撞将能量传递给离子和原子。由于电子质量很小，电子热能比动能大得多，忽略电子动能的影响，并考虑能量的起源与耗散。在该模型中，不考虑辐射能量损失，得到如下的电子能量方程

（9）

J=ene（u-ue）

代入式（9），并使用广义欧姆定律消去电场强度，得到

（10）

1.3 电离模型

在高电流的工况条件下，内部等离子体通过焦耳加热迅速升温，该温度可能突破粒子的二次电离甚至更高电离需要的能量。因此，需要引入电离模型对MPDT内部电离情况进行描述，考虑氩气电离形成1～6价的离子，每种离子均满足连续性方程，不考虑离子的迁移和激发态离子，得到离子的连续性方程如下

（11）

其中，ns为s价态离子数密度；源项ωs主要为电离与复合，其计算方式如下

（12）

其中，kf,i、kb,i（i=1，…，6）分别为反应速率和复合速率，可由文献[15]求得。

1.4 MPDT中等离子体的性质

（1）霍尔参数

在MPDT中，由于强磁场的存在以及低压电弧放电，使得等离子体中存在很强的霍尔效应，离子和电子的霍尔效应可以用离子霍尔系数βi和电子霍尔系数βe来评估

（13）

（14）

式中，ωci、ωce分别为离子和电子的回旋频率；vii、vin、ven分别为离子-离子碰撞频率、离子-原子碰撞频率和电子-原子碰撞频率；vjz为组分j和组分z的碰撞频率；mj为组分j的粒子质量；nz为组分z的数密度；Qjz为组分j和组分z的碰撞截面；Tj为组分j的温度。当Te=Ti=3 eV，ne=1021/m3时，βi和βe分别等于1.5×10-2和4。这表明，离子的霍尔效应相较于电子的霍尔效应是可以忽略的。此外，若工质数密度在某些区域显著降低，需要考虑离子滑移的影响。由于在MPDT中，推进剂几乎是完全电离的，因此离子滑移并不重要。

（2）碰撞频率

带电粒子之间的碰撞可由库伦碰撞进行描述[16]。其中碰撞类型包括电子-电子碰撞、电子-离子碰撞、电子-原子碰撞、离子-离子碰撞以及离子-原子碰撞。其电子-离子碰撞截面Qei可由式（15）表示

（15）

其中，Zi为离子价态；Z为有效电离分数。对于氩气，其电子-中性原子以及离子-中性原子的碰撞截面可认为是常数[17]，即Qeo≈4.0×10-20m2，Qio≈1.4×10-18m2。电子-电子碰撞可认为与电子-离子碰撞等价。离子-离子碰撞截面为

（16）

电子-电子以及原子和离子-电子碰撞频率如下

（17）

（18）

其中，Qes、Qee分别为电子-离子碰撞截面和电子-电子碰撞截面；me为电子质量。

离子和中性原子与离子的碰撞频率可表示为

（19）

其中，Qis为离子与原子或离子的碰撞截面；mi为离子质量。

（3）电阻率

根据经典等离子体理论，等离子体中的电阻率可由碰撞频率表示为

（20）

由于在MPDT中，交叉场引起的等离子体不稳定性形成的湍流波动,会导致等离子体电阻率异常波动，从而影响有效电导率[18]，因此需要考虑反常碰撞引起的反常电阻率。根据Choueiri的研究表明[19]，当电子漂移速度与离子热运动速度满足以下关系式时,需要考虑反常碰撞的影响，即

（21）

其中，ude为电子漂移速度；vti为离子热运动速度。

反常碰撞率van与经典碰撞率vej的比值取决于经典电子霍尔参数和离子与电子的温度比值，其表示如下

（22）

因此，有效电导率可以由反常碰撞频率和经典碰撞频率共同表示如下

（23）

（4）热导率

根据文献[16]，等离子体的重粒子热导率λh和电子热导率λe可以由碰撞频率表示为如下形式。认为重粒子温度与等离子体的黏性相关，而电子质量较小，其导热系数主要来源于粒子碰撞。

λh=μCv,

（24）

式中，μ是黏性系数；Cv是等离子体的定容比热。

μ=

（25）

（26）

其中，U为气体内能，可由如下表达式求出，Zj为组分j的配分函数

（27）

2 数值方法

对于等离子体流动控制方程，使用有限体积法对二维柱坐标下的守恒型方程进行描述，其一般形式可表示为

（28）

其中，U为守恒型物理量；F、G为对流项的物理通量；S为源项。

在该研究中，使用TVD Lax-Friedrich格式对对流项进行离散求解，在空间上使用MUSCL差值方法获得二阶精度，在时间上使用预测校正方法获得二阶精度。对源项使用二阶中心差分即可。

在时间n+1/2,空间网格（i,j）处构建式（28）的有限差分格式，其中预测步为

（29）

其中，Δξ、Δη、Δt分别表示ξ方向空间步长、η方向空间步长以及时间步长。

在校正步中对通量进行重建，使用minmod限制器构建左右守恒量通量（UL,UR）：

ξ方向的守恒量通量为

（30）

（31）

故校正步可表示为

（32）

TVD Lax-Friedrich形式

（33）

修正TVD Lax-Friedrich形式

（34）

其中，|aj+1/2|max为通量雅可比矩阵的最大特征值，该特征值为当地速度加上快磁声速。由于修正TVD Lax-Friedrich的耗散限制器效果较为出色，拟采用该修正形式。

同理，η方向采用相同的处理方式，得到该方向的数值通量，并代入校正步得到n+1时间步的守恒量。

使用交替方向隐式差分格式（Alternating Direction Implicit method，ADI）对磁感应方程进行求解，该方法无条件稳定并能较快收敛。

3 几何构型及边界条件

针对阴阳极同轴的自身场推力器建立对称轴形式的计算域，只考虑其中横截面一半的部分，其结构模型以及计算域如图1所示，其中阳极内径55mm，阳极长度与阴极长度相同，均为20cm，阴极半径10mm。为了更好地捕捉羽流中等离子体流动特性，将计算域的长扩展到45cm，径向扩展到10cm。推进剂从阴阳极之间的背板流入推进器内部，并假设推进剂以均匀的方式注入。在阴阳极高电流的作用下，推进剂注入推进器内部后在很短的距离被迅速电离，在该区域内准中性假设并不能严格成立，并且有理论认为该区域等离子体的不稳定性产生了非麦克斯韦速度分布的电子[20]，该物理过程并不能很好地用流体理论描述，因此舍弃等离子体被迅速电离的过程以及距离，认为推进剂以完全电离的形式注入。

图1 结构模型及计算域Fig.1 Structural model and computational domain

4 仿真结果及讨论

使用氩气作为推力器的工质，并以6g/s的质量流量注入到推力器内部，并且在阴阳极之间施加15kA的放电电流以产生较高磁场。通过求解非定常流场控制方程并耦合磁场，得到等离子体流场轴向速度分布如图2所示。

其中，图2所示为等离子体流场轴向速度分布，等离子体羽流最大速度达到了12000m/s，在径向洛伦兹力的作用下，羽流整体速度分布更加偏向于对称轴;并且由于膨胀以及轴向洛伦兹力的作用，在阴极表面等离子体逐渐被加速，在阴极出口处已经被加速到11000m/s。但由于羽流区的磁场离放电电流较远，磁场很难扩散到羽流，因此羽流区磁感应强度较小，使得轴向洛伦兹力减小，加速效果将明显减弱。

图3 磁场等值线分布（与最大值的比值）Fig.3 Magnetic field contour distribution （ratio to maximum value）

图3所示为磁场的空间分布状况，并以最大值的比值展示出来，其中阴极根部的磁感应强度最大，达到0.33T，并且可以很明显地看到沿轴线方向磁场的线性分布。此外。在阴阳极出口平面，磁感应强度已经降到了入口磁场的5%，意味着在推力器内部，磁场具有很大的梯度。羽流区的磁感应强度较小，意味着洛伦兹力的减弱。

电子温度随空间的分布如图4所示，阴极尖端展现出较强的加热情况，这与等离子体流动特性相关，由于对称轴的作用，等离子体在阴极尖端产生停滞作用，将动能转化为内能。此外，由于阴极尖端处于比较靠近对称轴的地方，该处电流受尖端放电等影响产生比较高的电流密度，因此，阴极尖端处产生的焦耳加热同样导致温度的升高，并向下游扩散。阳极离对称轴较远，其尖端放电电流较小，该处焦耳加热并不是很明显，因此并没有出现较为明显的加热情况。停滞作用对于大质量的离子来说比碰撞升温更为明显，因此在阴极尖端离子也会产生极高的温度，高达12eV。

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

图5 霍尔参数分布Fig.5 Hall parameter distribution

电子霍尔效应分布如图5所示，其主要分布在阴极上表面以及阳极侧表面。根据文献[21]中对于霍尔参数的理论分析，其主要与磁感应强度和温度有关，磁场越强，温度越高，电子霍尔参数越高，这在阴极表面得到了很好的验证。但阳极侧面附近的高电子霍尔参数主要与数密度有关，该处的数密度较低，碰撞频率较低，同样也会导致霍尔参数过高。霍尔参数过高会使得电流线更加弯曲，压降将会增加，不利于性能的提升。

图6所示为二价氩离子分布云图，由于阴极尖端的极高温度，导致该处的电离急剧加强，二次电离明显增强，但温度没有达到可以大规模产生三次电离的程度。因此，可认为在高电流的工况下，工质在推力器内部几乎完全电离并且能够产生高次电离。在羽流下半区域，温度较低，数密度也偏低，使得该区域复合形成原子数密度并不高，因此，复合绝大部分产生于温度较低以及数密度较高的阳极附近。

图6 二价氩离子分布Fig.6 Ar++ distribution

5 结 论

本文通过使用流体耦合电磁场的方法,对特定构型以及特定工况的自身场磁等离子体推力器内部及羽流进行数值仿真，对流体控制方程使用计算量较小的TVD Lax-Friedrich格式，对电磁场使用ADI方法进行求解，得到了推力器等离子体流场空间分布特性。仿真结果与实验较为吻合，能够解释SF-MPDT中出现的一些物理现象。

通过流场分布情况可以看出，洛伦兹力对流场的约束作用，使得流线更加偏向于对称轴，这是有利于产生有效轴向推力的。此外，由于磁场以及数密度的影响，阴极根部以及阳极侧面会产生较高的霍尔参数，该参数导致电流线的延长，增加了压降，这不利于推力器产生更高效率。针对推力器内部等离子体完全电离的假设基本成立，在高放电电流的作用下，产生较高温度使得推力器内部等离子体几乎完全电离，并且在阴极尖端出现较高的二次电离。

本文初步分析了自身场磁等离子体推力器内部流动状况，但由于模型使用较少，需要添加更符合等离子体运动规律的模型，以减少使用的假设使之接近实际发生的物理过程。