自身場磁等離子體推力器數值仿真流場特性分析*

操樂暉，吳 鵬，羅梓浩，任軍學，湯海濱,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中出現的一些物理現象。

通過流場分布情況可以看出，洛倫茲力對流場的約束作用，使得流線更加偏向于對稱軸，這是有利于產生有效軸向推力的。此外，由于磁場以及數密度的影響，陰極根部以及陽極側面會產生較高的霍爾參數，該參數導致電流線的延長，增加了壓降，這不利于推力器產生更高效率。針對推力器內部等離子體完全電離的假設基本成立，在高放電電流的作用下，產生較高溫度使得推力器內部等離子體幾乎完全電離，并且在陰極尖端出現較高的二次電離。

本文初步分析了自身場磁等離子體推力器內部流動狀況，但由于模型使用較少，需要添加更符合等離子體運動規律的模型，以減少使用的假設使之接近實際發生的物理過程。