操樂暉，吳 鵬，羅梓浩，任軍學(xué)，湯海濱,3,4

（1.北京航空航天大學(xué) 空間與環(huán)境學(xué)院·北京·102206；2.北京航空航天大學(xué) 宇航學(xué)院·北京·102206；3.航天器設(shè)計優(yōu)化與動態(tài)模擬技術(shù)教育部重點實驗室·北京·102206；4.空間環(huán)境監(jiān)測與信息處理工業(yè)和信息化部重點實驗室·北京·102206）

0 引 言

磁等離子體推力器（Magnetoplasmadynamic Thruster，MPDT）利用高電流對推進劑進行電離，并產(chǎn)生較高電離度的等離子體，再通過電磁場與等離子體相互作用產(chǎn)生洛倫茲力，將等離子體加速噴出。從目前已有的關(guān)于MPDT的研究來看，MPDT比沖的量級最高可達1萬s以上，相較于傳統(tǒng)的化學(xué)推進要高出10～20倍；此外，推力量級最高能達到100多N，遠大于靜電式推力器和電熱式推力器[1]。MPDT的性能特點使得該推力器在未來嚴峻的深空探測任務(wù)中將扮演重要角色，由于MPDT需要的電功率偏大，結(jié)合未來空間核電技術(shù)，發(fā)展大功率的核電磁等離子體推力器將大大提高航天器的性能，減少發(fā)射成本，以及可進行更長時間更遠距離的飛行任務(wù)[2]。

MPDT根據(jù)磁場的來源不同，可分為自身場磁等離子體推力器（Self-Field Magnetoplasmady-namic Thruster，SF-MPDT）和附加場磁等離子體推力器（Applied-Field Magnetoplasmadynamic Thruster，AF-MPDT）。SF-MPDT的磁場來源主要由施加在陰陽極之間的電流產(chǎn)生，該電流會在推力器內(nèi)部及羽流處感應(yīng)出角向磁場，該角向磁場與等離子體中的徑向電流作用產(chǎn)生沿軸向的洛倫茲力jrBθ。通常SF-MPDT需要比較高的電流才能在推力器內(nèi)部產(chǎn)生比較高的磁場，其電流需要在幾千A甚至是十幾千A才能獲得比較好的性能參數(shù)，因此SF-MPDT的功率會比較高。AF-MPDT的磁場來源于外加線圈和永磁體，并在推力器內(nèi)部及羽流處產(chǎn)生軸向和徑向的磁場，該磁場與等離子體相互作用，產(chǎn)生更加復(fù)雜的推力產(chǎn)生機制以及等離子體流動行為。一般情況下，通電線圈或永磁體在推力器內(nèi)部產(chǎn)生零點幾T的磁場，此外陰陽極之間施加的電流較自身場要小很多，其量級在幾百A左右。

對于SF-MPDT，其加速機制比較簡單，主要有氣動加速和電磁加速，陰陽極之間的高電流將推進劑加熱并電離，形成高溫高密度的等離子體，通過噴管碰撞噴出，該方式為氣動加速；電磁加速主要來源于電流與磁場的相互作用，根據(jù)廣義歐姆定律（式（1））可將電磁加速機制歸結(jié)為電場效應(yīng)（等式右邊第一項）、流動效應(yīng)（等式右邊第二項）、電子壓力效應(yīng)（等式右邊第三項）和霍爾效應(yīng)（等式右邊第四項）四種作用。

（1）

其中，J為等離子體電流；E為等離子體內(nèi)部電場；ne為電子數(shù)密度；B為磁感應(yīng)強度；σ為電導(dǎo)率；e為電子電荷量；pe為電子壓力梯度；u為流體速度。

附加場MPDT由于存在徑向磁場和軸向磁場，其推力產(chǎn)生機制有如下幾種：

1）氣動加速。其機制與SF-MPDT氣動加速類似。

2）自身場加速。由于AF-MPDT也需要在陰陽極之間施加電流，該電流也會在推力器內(nèi)部感應(yīng)出角向磁場，該角向磁場與電流相互作用產(chǎn)生軸向的洛倫茲力。但通常AF-MPDT的電流偏小，感應(yīng)出的角向磁場偏小，與附加場相比幾乎可以忽略，該加速機制在總體加速機制中占比相對較小。

3）渦旋加速。在附加磁場和電流的相互作用下，產(chǎn)生周向洛倫茲力jrBz和jzBr，會帶動等離子體產(chǎn)生周向的運動。通過具有磁鏡效應(yīng)的磁噴管，將高速旋轉(zhuǎn)等離子體的旋轉(zhuǎn)動能轉(zhuǎn)化為軸向動能。實驗中觀察到這一能量轉(zhuǎn)化機制通常被認為是附加場的主要加速機制。

4）霍爾加速。在強磁場和低質(zhì)量流量的工況條件下，根據(jù)廣義歐姆定律，角向會感應(yīng)出比較強的電流，該電流與附加磁場相互作用，產(chǎn)生的軸向洛倫茲力jθBr作用在等離子體上。對于霍爾加速，目前存在比較多的爭議，主要是該機制產(chǎn)生的力對推力的貢獻方向并不是很明確[3]。

目前已有理論對SF-MPD內(nèi)部等離子體參數(shù)進行預(yù)測[4]，此外在實驗上建立了國內(nèi)大功率的附加場MPDT設(shè)計以及推力測量手段[5-6]，并在此基礎(chǔ)上開展了推力器內(nèi)部物理機制的研究[7]。但實驗中由于推力器內(nèi)部等離子體的惡劣環(huán)境，導(dǎo)致測量方面存在困難，并不能很好地解釋推力器內(nèi)等離子體產(chǎn)生的物理現(xiàn)象，如onset現(xiàn)象等[8]。因此，需要在數(shù)值模擬的基礎(chǔ)上，對推力器內(nèi)部的等離子體進行描述。目前，描述等離子體行為的方法主要有：粒子（Particlein Cell，PIC）方法以及磁流體（Magnetohydrodynamics，MHD）方法。PIC方法通過追蹤推力器內(nèi)的微觀粒子，研究其在電磁場中的運動情況[9-11]。由于MPDT中數(shù)密度相對較高，該方法將會非常耗時。MHD方法假設(shè)等離子體為導(dǎo)電流體，通過求解流體方程得到等離子體的流場特性，雖然不能描述等離子體的微觀特性，但對于宏觀特性，則能夠很好地進行預(yù)測。

本文使用MHD方法，對自身場磁等離子體推力器內(nèi)部及羽流進行建模，利用高精度數(shù)值方法求解流場和磁場，并進行耦合求解。

1 物理模型

1.1 流體模型

當粒子的速度分布達到麥克斯韋速度分布的時間遠小于流場的特征時間時，可以將等離子體作為流體進行處理，描述該近似可以用Knudsen數(shù)表示[12]，即

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

碰撞截面與碰撞類型以及溫度有關(guān)，且推力器的幾何構(gòu)型是一定的，因此其大小主要與推力器中的等離子體數(shù)密度有關(guān)。自身場推力器的數(shù)密度能達到1021～1022/m3,且由于感應(yīng)磁場垂直于流場，使得碰撞更加激烈，縮短了達到麥克斯韋速度分布的時間，更加符合流體近似。此外，波與粒子相互作用產(chǎn)生的不穩(wěn)定性,導(dǎo)致反常碰撞進一步滿足流體近似。

根據(jù)自身場MPDT內(nèi)部等離子體流動情況，使用以下假設(shè)條件建立物理模型：

1）內(nèi)部流動為定常、軸對稱、層流、無黏流動；

2）準中性假設(shè)，不存在電荷分離的現(xiàn)象；

3）熱力學(xué)非平衡假設(shè)，即雙溫度模型；

4）忽略輻射能量損失；

5）認為電導(dǎo)率為標量，并考慮反常電導(dǎo)率；

6）考慮霍爾效應(yīng)，但忽略離子滑移。

對于流場的描述，使用以下控制方程：

連續(xù)性方程

（2）

動量方程

（3）

總能量方程

（4）

電磁場計算方法由麥克斯韋方程組導(dǎo)出。麥克斯韋方程組微分形式如下

（5）

其中，ε0為真空介電常數(shù)；μ0為真空磁導(dǎo)率。

由于準中性假設(shè)，認為無電荷分離。此外，對于絕大部分問題，等離子體頻率通常高達幾十MHz，位移電流比傳導(dǎo)電流小很多，可以忽略不計，因此可以得到簡化的麥克斯韋方程組如下

（6）

并結(jié)合廣義歐姆定律（式（1）），消去電場項，得到磁感應(yīng)方程如下

（7）

其中，σ為等離子體電導(dǎo)率，由于自身場磁場比較單一，該電導(dǎo)率可處理為標量。

流體與磁場的耦合主要通過洛倫茲力和焦耳加熱聯(lián)系起來，因此可以通過求解流場參數(shù)分布求解磁場。即通過磁感應(yīng)方程求解磁感應(yīng)強度B，再由安培環(huán)路定律求出電流密度J，并由廣義歐姆定律得到電場強度E，代入流場的源項中，從而使得磁場與流場聯(lián)系起來。

但以上方程組并不封閉，需要引入氣體狀態(tài)方程

p=nekTe+nikTi

（8）

其中，ne、ni分別為電子和重粒子數(shù)密度；k為玻耳茲曼常數(shù)；Te、Ti分別為電子溫度和重粒子溫度。

1.2 熱力學(xué)非平衡模型

要想維持等離子體內(nèi)部熱平衡態(tài)，需要碰撞的時間遠小于等離子體流動特征時間，才能使得等離子體維持麥克斯韋速度分布。在MPDT中，數(shù)密度的數(shù)量級只有o（21）～o（22），而等離子體宏觀流速一般是超聲速的，經(jīng)計算MPDT內(nèi)部等離子體的停留時間與達到平衡的時間相當[13]，不頻繁的碰撞使得電子和離子并不能達到局部熱平衡，因此需要將離子和電子分開考慮。建立非平衡模型，即雙溫度模型。此外，實驗結(jié)果也觀察到了這一非平衡態(tài)[14]。

假設(shè)電磁場通過對電子的作用對電子進行供能，電子通過碰撞將能量傳遞給離子和原子。由于電子質(zhì)量很小，電子熱能比動能大得多，忽略電子動能的影響，并考慮能量的起源與耗散。在該模型中，不考慮輻射能量損失，得到如下的電子能量方程

（9）

J=ene（u-ue）

代入式（9），并使用廣義歐姆定律消去電場強度，得到

（10）

1.3 電離模型

在高電流的工況條件下，內(nèi)部等離子體通過焦耳加熱迅速升溫，該溫度可能突破粒子的二次電離甚至更高電離需要的能量。因此，需要引入電離模型對MPDT內(nèi)部電離情況進行描述，考慮氬氣電離形成1～6價的離子，每種離子均滿足連續(xù)性方程，不考慮離子的遷移和激發(fā)態(tài)離子，得到離子的連續(xù)性方程如下

（11）

其中，ns為s價態(tài)離子數(shù)密度；源項ωs主要為電離與復(fù)合，其計算方式如下

（12）

其中，kf,i、kb,i（i=1，…，6）分別為反應(yīng)速率和復(fù)合速率，可由文獻[15]求得。

1.4 MPDT中等離子體的性質(zhì)

（1）霍爾參數(shù)

在MPDT中，由于強磁場的存在以及低壓電弧放電，使得等離子體中存在很強的霍爾效應(yīng)，離子和電子的霍爾效應(yīng)可以用離子霍爾系數(shù)βi和電子霍爾系數(shù)βe來評估

（13）

（14）

式中，ωci、ωce分別為離子和電子的回旋頻率；vii、vin、ven分別為離子-離子碰撞頻率、離子-原子碰撞頻率和電子-原子碰撞頻率；vjz為組分j和組分z的碰撞頻率；mj為組分j的粒子質(zhì)量；nz為組分z的數(shù)密度；Qjz為組分j和組分z的碰撞截面；Tj為組分j的溫度。當Te=Ti=3 eV，ne=1021/m3時，βi和βe分別等于1.5×10-2和4。這表明，離子的霍爾效應(yīng)相較于電子的霍爾效應(yīng)是可以忽略的。此外，若工質(zhì)數(shù)密度在某些區(qū)域顯著降低，需要考慮離子滑移的影響。由于在MPDT中，推進劑幾乎是完全電離的，因此離子滑移并不重要。

（2）碰撞頻率

帶電粒子之間的碰撞可由庫倫碰撞進行描述[16]。其中碰撞類型包括電子-電子碰撞、電子-離子碰撞、電子-原子碰撞、離子-離子碰撞以及離子-原子碰撞。其電子-離子碰撞截面Qei可由式（15）表示

（15）

其中，Zi為離子價態(tài)；Z為有效電離分數(shù)。對于氬氣，其電子-中性原子以及離子-中性原子的碰撞截面可認為是常數(shù)[17]，即Qeo≈4.0×10-20m2，Qio≈1.4×10-18m2。電子-電子碰撞可認為與電子-離子碰撞等價。離子-離子碰撞截面為

（16）

電子-電子以及原子和離子-電子碰撞頻率如下

（17）

（18）

其中，Qes、Qee分別為電子-離子碰撞截面和電子-電子碰撞截面；me為電子質(zhì)量。

離子和中性原子與離子的碰撞頻率可表示為

（19）

其中，Qis為離子與原子或離子的碰撞截面；mi為離子質(zhì)量。

（3）電阻率

根據(jù)經(jīng)典等離子體理論，等離子體中的電阻率可由碰撞頻率表示為

（20）

由于在MPDT中，交叉場引起的等離子體不穩(wěn)定性形成的湍流波動,會導(dǎo)致等離子體電阻率異常波動，從而影響有效電導(dǎo)率[18]，因此需要考慮反常碰撞引起的反常電阻率。根據(jù)Choueiri的研究表明[19]，當電子漂移速度與離子熱運動速度滿足以下關(guān)系式時,需要考慮反常碰撞的影響，即

（21）

其中，ude為電子漂移速度；vti為離子熱運動速度。

反常碰撞率van與經(jīng)典碰撞率vej的比值取決于經(jīng)典電子霍爾參數(shù)和離子與電子的溫度比值，其表示如下

（22）

因此，有效電導(dǎo)率可以由反常碰撞頻率和經(jīng)典碰撞頻率共同表示如下

（23）

（4）熱導(dǎo)率

根據(jù)文獻[16]，等離子體的重粒子熱導(dǎo)率λh和電子熱導(dǎo)率λe可以由碰撞頻率表示為如下形式。認為重粒子溫度與等離子體的黏性相關(guān)，而電子質(zhì)量較小，其導(dǎo)熱系數(shù)主要來源于粒子碰撞。

λh=μCv,

（24）

式中，μ是黏性系數(shù)；Cv是等離子體的定容比熱。

μ=

（25）

（26）

其中，U為氣體內(nèi)能，可由如下表達式求出，Zj為組分j的配分函數(shù)

（27）

2 數(shù)值方法

對于等離子體流動控制方程，使用有限體積法對二維柱坐標下的守恒型方程進行描述，其一般形式可表示為

（28）

其中，U為守恒型物理量；F、G為對流項的物理通量；S為源項。

在該研究中，使用TVD Lax-Friedrich格式對對流項進行離散求解，在空間上使用MUSCL差值方法獲得二階精度，在時間上使用預(yù)測校正方法獲得二階精度。對源項使用二階中心差分即可。

在時間n+1/2,空間網(wǎng)格（i,j）處構(gòu)建式（28）的有限差分格式，其中預(yù)測步為

（29）

其中，Δξ、Δη、Δt分別表示ξ方向空間步長、η方向空間步長以及時間步長。

在校正步中對通量進行重建，使用minmod限制器構(gòu)建左右守恒量通量（UL,UR）：

ξ方向的守恒量通量為

（30）

（31）

故校正步可表示為

（32）

TVD Lax-Friedrich形式

（33）

修正TVD Lax-Friedrich形式

（34）

其中，|aj+1/2|max為通量雅可比矩陣的最大特征值，該特征值為當?shù)厮俣燃由峡齑怕曀?。由于修正TVD Lax-Friedrich的耗散限制器效果較為出色，擬采用該修正形式。

同理，η方向采用相同的處理方式，得到該方向的數(shù)值通量，并代入校正步得到n+1時間步的守恒量。

使用交替方向隱式差分格式（Alternating Direction Implicit method，ADI）對磁感應(yīng)方程進行求解，該方法無條件穩(wěn)定并能較快收斂。

3 幾何構(gòu)型及邊界條件

針對陰陽極同軸的自身場推力器建立對稱軸形式的計算域，只考慮其中橫截面一半的部分，其結(jié)構(gòu)模型以及計算域如圖1所示，其中陽極內(nèi)徑55mm，陽極長度與陰極長度相同，均為20cm，陰極半徑10mm。為了更好地捕捉羽流中等離子體流動特性，將計算域的長擴展到45cm，徑向擴展到10cm。推進劑從陰陽極之間的背板流入推進器內(nèi)部，并假設(shè)推進劑以均勻的方式注入。在陰陽極高電流的作用下，推進劑注入推進器內(nèi)部后在很短的距離被迅速電離，在該區(qū)域內(nèi)準中性假設(shè)并不能嚴格成立，并且有理論認為該區(qū)域等離子體的不穩(wěn)定性產(chǎn)生了非麥克斯韋速度分布的電子[20]，該物理過程并不能很好地用流體理論描述，因此舍棄等離子體被迅速電離的過程以及距離，認為推進劑以完全電離的形式注入。

圖1 結(jié)構(gòu)模型及計算域Fig.1 Structural model and computational domain

4 仿真結(jié)果及討論

使用氬氣作為推力器的工質(zhì)，并以6g/s的質(zhì)量流量注入到推力器內(nèi)部，并且在陰陽極之間施加15kA的放電電流以產(chǎn)生較高磁場。通過求解非定常流場控制方程并耦合磁場，得到等離子體流場軸向速度分布如圖2所示。

其中，圖2所示為等離子體流場軸向速度分布，等離子體羽流最大速度達到了12000m/s，在徑向洛倫茲力的作用下，羽流整體速度分布更加偏向于對稱軸;并且由于膨脹以及軸向洛倫茲力的作用，在陰極表面等離子體逐漸被加速，在陰極出口處已經(jīng)被加速到11000m/s。但由于羽流區(qū)的磁場離放電電流較遠，磁場很難擴散到羽流，因此羽流區(qū)磁感應(yīng)強度較小，使得軸向洛倫茲力減小，加速效果將明顯減弱。

圖3 磁場等值線分布（與最大值的比值）Fig.3 Magnetic field contour distribution （ratio to maximum value）

圖3所示為磁場的空間分布狀況，并以最大值的比值展示出來，其中陰極根部的磁感應(yīng)強度最大，達到0.33T，并且可以很明顯地看到沿軸線方向磁場的線性分布。此外。在陰陽極出口平面，磁感應(yīng)強度已經(jīng)降到了入口磁場的5%，意味著在推力器內(nèi)部，磁場具有很大的梯度。羽流區(qū)的磁感應(yīng)強度較小，意味著洛倫茲力的減弱。

電子溫度隨空間的分布如圖4所示，陰極尖端展現(xiàn)出較強的加熱情況，這與等離子體流動特性相關(guān)，由于對稱軸的作用，等離子體在陰極尖端產(chǎn)生停滯作用，將動能轉(zhuǎn)化為內(nèi)能。此外，由于陰極尖端處于比較靠近對稱軸的地方，該處電流受尖端放電等影響產(chǎn)生比較高的電流密度，因此，陰極尖端處產(chǎn)生的焦耳加熱同樣導(dǎo)致溫度的升高，并向下游擴散。陽極離對稱軸較遠，其尖端放電電流較小，該處焦耳加熱并不是很明顯，因此并沒有出現(xiàn)較為明顯的加熱情況。停滯作用對于大質(zhì)量的離子來說比碰撞升溫更為明顯，因此在陰極尖端離子也會產(chǎn)生極高的溫度，高達12eV。

圖4 電子溫度分布Fig.4 Electron temperature distribution

圖5 霍爾參數(shù)分布Fig.5 Hall parameter distribution

電子霍爾效應(yīng)分布如圖5所示，其主要分布在陰極上表面以及陽極側(cè)表面。根據(jù)文獻[21]中對于霍爾參數(shù)的理論分析，其主要與磁感應(yīng)強度和溫度有關(guān)，磁場越強，溫度越高，電子霍爾參數(shù)越高，這在陰極表面得到了很好的驗證。但陽極側(cè)面附近的高電子霍爾參數(shù)主要與數(shù)密度有關(guān)，該處的數(shù)密度較低，碰撞頻率較低，同樣也會導(dǎo)致霍爾參數(shù)過高?；魻枀?shù)過高會使得電流線更加彎曲，壓降將會增加，不利于性能的提升。

圖6所示為二價氬離子分布云圖，由于陰極尖端的極高溫度，導(dǎo)致該處的電離急劇加強，二次電離明顯增強，但溫度沒有達到可以大規(guī)模產(chǎn)生三次電離的程度。因此，可認為在高電流的工況下，工質(zhì)在推力器內(nèi)部幾乎完全電離并且能夠產(chǎn)生高次電離。在羽流下半?yún)^(qū)域，溫度較低，數(shù)密度也偏低，使得該區(qū)域復(fù)合形成原子數(shù)密度并不高，因此，復(fù)合絕大部分產(chǎn)生于溫度較低以及數(shù)密度較高的陽極附近。

圖6 二價氬離子分布Fig.6 Ar++ distribution

5 結(jié) 論

本文通過使用流體耦合電磁場的方法,對特定構(gòu)型以及特定工況的自身場磁等離子體推力器內(nèi)部及羽流進行數(shù)值仿真，對流體控制方程使用計算量較小的TVD Lax-Friedrich格式，對電磁場使用ADI方法進行求解，得到了推力器等離子體流場空間分布特性。仿真結(jié)果與實驗較為吻合，能夠解釋SF-MPDT中出現(xiàn)的一些物理現(xiàn)象。

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

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