叢云天，陳 健

(北京控制工程研究所，北京 100190)

*總裝預(yù)研基金資助項(xiàng)目(9140A20050315HT05002).

霍爾推力器仿真與優(yōu)化*

叢云天，陳 健

(北京控制工程研究所，北京 100190)

針對霍爾推力器通道內(nèi)的放電過程，建立一種基于COMSOL軟件的仿真模型.該模型以等離子體內(nèi)部電子和離子的漂移擴(kuò)散為核心，結(jié)合電磁場、氣體流動以及等離子體內(nèi)部的碰撞反應(yīng)，通過合理選取系統(tǒng)方程、邊界條件以及求解器，有效估算霍爾推力器的性能參數(shù)以及各物理量在通道內(nèi)的分布.將仿真結(jié)果與理論相比較，驗(yàn)證該模型用于霍爾推力器數(shù)值仿真的有效性并以此對推力器進(jìn)行磁場優(yōu)化設(shè)計.

霍爾推力器；等離子體；漂移擴(kuò)散；有限元

0 引 言

霍爾推力器(Hall effect thrust，HET)自20世紀(jì)70年代在前蘇聯(lián)流星號衛(wèi)星上應(yīng)用以來，經(jīng)過幾十年的發(fā)展日趨成熟，出現(xiàn)了多款經(jīng)典的飛行產(chǎn)品，如SPT-100、PPS-1350G、BPT-4000等.霍爾推力器由空心陰極、放電室、磁場生成器(內(nèi)/外磁線圈、磁極)以及陽極/氣體分配器等結(jié)構(gòu)組成.霍爾推力器的工作原理如圖1所示，內(nèi)外磁線圈在環(huán)形通道內(nèi)產(chǎn)生徑向磁場，陽陰極之間的放電等離子體在通道內(nèi)將產(chǎn)生自洽的軸向電場.陰極發(fā)射的部分電子在流向陽極過程中遇到正交的徑向磁場與軸向的電場將形成E×B漂移(霍爾漂移).推進(jìn)劑通常采用惰性氣體氙(Xe)或氪(Kr)，從陽極/氣體分配器注入推力器通道，同做漂移運(yùn)動的電子發(fā)生碰撞并被電離，離子在軸向電場的作用下加速噴出產(chǎn)生反沖推力，同時電子通過傳導(dǎo)機(jī)制到達(dá)陽極，在通道內(nèi)實(shí)現(xiàn)穩(wěn)定的等離子放電過程，提供持續(xù)穩(wěn)定的推力[1].

國外針對霍爾推力器的仿真建模已開展了大量研究，結(jié)合網(wǎng)格粒子仿真(particle in cell，PIC)、蒙特卡洛(MCC)或直接蒙特卡洛碰撞法(DSMC)等方法對放電室工作過程、壁面濺射腐蝕和羽流發(fā)散進(jìn)行仿真，并構(gòu)建了較為成熟的仿真代碼或程序，如HPHall(Hybrid-PIC Hall)、PICPluS3D等[2-3].

本文采用有限元分析軟件COMSOL Multiphysics對霍爾推力器放電通道建模仿真，對比分析1.5 kW 級霍爾推力器的數(shù)值仿真結(jié)果，驗(yàn)證軟件和數(shù)值模擬方法的可行性.在此基礎(chǔ)上，進(jìn)一步提出推力器的優(yōu)化設(shè)計.

1 物理數(shù)學(xué)模型

考慮到霍爾推力器放電室工作機(jī)理，對推力器內(nèi)部等離子體性質(zhì)的求解主要基于COMSOL軟件的直流放電模塊，該模塊引用一種新型的等離子體模擬混合模型，既具有流體模型的較快計算速度同時又適用于低壓等離子體模擬.此外，COMSOL適用于對等離子體動力學(xué)、推進(jìn)劑、電磁場進(jìn)行耦合計算.根據(jù)推力器放電過程的主要工作機(jī)理，采用合理的系統(tǒng)方程或控制方程描述等離子體內(nèi)部反應(yīng)機(jī)制以及各粒子的遷移輸運(yùn)性質(zhì).

1.1系統(tǒng)方程

1.1.1 電子擴(kuò)散

對于電子的運(yùn)動，近似采用漂移擴(kuò)散模型由電子數(shù)密度和平均電子能量來描述，通過求解電子連續(xù)性方程和能量守恒方程得到[4].

電子連續(xù)性方程為：

(1)

電子能量守恒方程計算電子能量密度：

(2)

式中，nε為電子能量密度，Rε為由于非彈性碰撞引起的能量損失或增加項(xiàng)，με為張量形式的電子能量遷移率，Dε為電子能量遷移率.

考慮到漂移擴(kuò)散方程固有的高度非線性，電子數(shù)密度在陽極附近等離子體鞘層區(qū)域內(nèi)很短的距離上就有10個數(shù)量級的跨度；在鞘層區(qū)域內(nèi)離子同電子的遷移率和擴(kuò)散性存在較大的差異，這種差異將形成一個與空間帶電隔離區(qū)，在此區(qū)域內(nèi)形成了一個能夠明顯提高電子能量的強(qiáng)電場區(qū)域.為解決這些問題，仿真中采用直接求解電子數(shù)密度和電子能量密度的對數(shù)[5].

存在磁場情況下，真實(shí)的電子遷移率的表達(dá)式不能用一個簡單的形式來表達(dá)，而應(yīng)該采用張量形式的電子遷移率

(3)

式中μdc為無外加磁場時的電子遷移率.

電子的擴(kuò)散系數(shù)、能量遷移率、能量擴(kuò)散系數(shù)分別為

(4)

式中，Te為電子溫度，定義Te=2nε/3ne.

1.1.2 重物質(zhì)擴(kuò)散

對于非電子粒子(重物質(zhì))，采用擴(kuò)散輸運(yùn)方程，假設(shè)計算中有k(k=1,…，Q)種類型物質(zhì)的流體以及j(j=1,…，N)種反應(yīng)類型，那么輸運(yùn)方程可表示為[5]：

(5)

式中，Jk為擴(kuò)散通量矢量，Rk為k種離子的比率表達(dá)式，u是中性粒子速度矢量，ρ是混合成分的密度，wk是第k種離子的質(zhì)量分?jǐn)?shù).

擴(kuò)散通量矢量Jk為

Jk=ρwkVk

(6)

Vk是k種成分總的擴(kuò)散速度矢量，采用混合平均模型計算，在Maxwell-Stefan公式基礎(chǔ)上進(jìn)行簡化，混合平均公式保留了質(zhì)量守恒定律且計算量較小，對擴(kuò)散系數(shù)Dk的簡化表達(dá)便于計算，但計算精度下降[6].Vk可表示為

(7)

(8)

而混合平均遷移率可通過愛因斯坦關(guān)系式計算

(9)

式中，q是單位電荷，kB為玻爾茲曼常數(shù).

1.1.3 磁場

J=σV×B+Je

(10)

(11)

其中，A為磁勢矢量，M為磁化強(qiáng)度矢量

1.1.4 N-S方程

(12)

式中，ρ為密度，u為速度矢量，P為壓強(qiáng)，F(xiàn)是體積力矢量，τ是粘性應(yīng)力張量，S是應(yīng)變率張量，μ是動態(tài)粘度.

1.2化學(xué)反應(yīng)

霍爾推力器通道內(nèi)部等離子體的粒子組成成分包括電子、離子和中性粒子(激發(fā)態(tài)和基態(tài))，氣體放電的本質(zhì)是帶電粒子與中性原子、激發(fā)態(tài)原子的相互碰撞，以及這些粒子與壁面之間的相互碰撞.電離主要發(fā)生在靠近推力器放電通道出口的電離區(qū)，中性氣體到達(dá)電離區(qū)后與高速運(yùn)動的電子發(fā)生碰撞，當(dāng)電子的能量高于中性氣體的電離能時，中性氣體可能發(fā)生電離.因?yàn)橥屏ζ髦须娮幽芰枯^低，高能電子碰撞電離很少發(fā)生，而低能電子與中性氣體碰撞后所生成的離子在加速電場的作用下噴出，所以二次電離很少發(fā)生，通常認(rèn)為推進(jìn)劑電離后主要生成一價離子.

上述全部反應(yīng)可歸為兩大類：①彈性碰撞，即粒子碰撞只改變運(yùn)動方向而粒子的總動量和動能均守恒，粒子本身能量未改變且無新的粒子或光子產(chǎn)生；②非彈性碰撞，碰撞過程引起了粒子內(nèi)能變換，或伴隨著新的粒子以及光子的產(chǎn)生[6].碰撞反應(yīng)的反應(yīng)速率由相應(yīng)的碰撞截面與電子能量分布函數(shù)(EEDF)進(jìn)行積分得到.霍爾推力器等離子體瞬態(tài)放電過程近似為平衡態(tài)系統(tǒng)的放電，則EEDF可取為麥克斯韋(Maxwell)分布函數(shù).而氙碰撞截面可根據(jù)Szabo[7]以經(jīng)驗(yàn)數(shù)據(jù)給出的分段函數(shù)計算.

1.3邊界條件

在霍爾推力器仿真過程中，COMSOL軟件在直流放電模塊和層流模塊中的邊界條件設(shè)置尤為重要，對仿真結(jié)果具有較大影響，需根據(jù)實(shí)際工況和理論進(jìn)行合理設(shè)置和選擇.

1.3.1 電子壁邊界條件

等離子體區(qū)域的壁面條件包括二次電子發(fā)射和熱發(fā)射.為表征等離子體中非電子類粒子的傳輸特性，將放電空間內(nèi)的固體邊界進(jìn)行邊界反應(yīng)設(shè)置.其主要的反應(yīng)包括在所有的等離子體放電空間固體邊界處均發(fā)生激發(fā)態(tài)氙原子和氙離子復(fù)原為基態(tài)原子的反應(yīng).而后者正是產(chǎn)生二次電子的反應(yīng)，該電子對于推力器通道內(nèi)放電過程尤為重要，正是由于轟擊產(chǎn)生的二次電子不斷地補(bǔ)充放電通道內(nèi)的電子損失，維持整個放電過程，因而應(yīng)設(shè)置邊界處的二次電子發(fā)射系數(shù)[5].

在電子的漂移擴(kuò)散過程中，對于放電空間內(nèi)的固體邊界要外加“壁”邊界條件.該種邊界條件主要用來描述電子在遇到固體邊界時產(chǎn)生電子密度變化的各種情況，其采用的電子交換機(jī)理是電子的損失取決于擴(kuò)散到邊界處的等離子體中的凈電子流量、在電子平均自由程內(nèi)與邊界條件的碰撞；電子密度的增加則取決于邊界處由于正離子的轟擊而引起的二次電子發(fā)射或邊界的熱發(fā)射(陰極).在壁上電子通量的法向分量方程以及電子能量密度的方程為：

(13)

其中，r是反射系數(shù)(一般情況下為0)，ve,th為熱運(yùn)動速度，γp是由離子引起的二次電子發(fā)射系數(shù)，Γp是離子通量，Γt是熱發(fā)射通量，εp是發(fā)射電子的平均能量，ε為平均熱離子能，n是向外的法線方向.熱運(yùn)動速度ve,th定義為：

(14)

1.3.2 層流邊界條件

等離子體內(nèi)部流動采用可壓縮層流計算，相對于穩(wěn)定的固體壁，層流在邊界上的速度均為u=0.因此，模型計算中采用無滑移壁.入口處應(yīng)用質(zhì)量流邊界條件，出口處一般通過法向應(yīng)力條件規(guī)定，給定出口的壓力邊界條件求解內(nèi)部速度.

1.4物理模型

針對1.5 kW級磁聚焦型(ATON)霍爾推力器進(jìn)行二維軸對稱建模，其簡化后的幾何尺寸結(jié)構(gòu)參見圖2.完整的霍爾推力器包括陽極、陰極、氣體分配器、通道套筒及磁路等，這里陰極簡化為出口平面處的邊界，并給定熱發(fā)射電子通量.在模型中設(shè)置各部分的材料，陽極采用1Cr18Ni9Ti，磁路選擇非線性磁性材料，勵磁線圈采用Cu，陶瓷套筒采用六方氮化硼(BN).

對于物理場的選擇，考慮層流模塊、磁場模塊以及直流放電模塊，根據(jù)霍爾推力器內(nèi)部反應(yīng)的機(jī)理選擇相應(yīng)幾何區(qū)域，并設(shè)置相應(yīng)的系統(tǒng)方程、化學(xué)反應(yīng)和邊界條件及其相關(guān)參數(shù).之后對仿真區(qū)域進(jìn)行網(wǎng)格劃分，在等離子體放電區(qū)、近陽極區(qū)和出口處采用較為細(xì)化的網(wǎng)格，其他僅有磁場的區(qū)域可適當(dāng)劃分.最終計算區(qū)域的網(wǎng)格劃分如圖3所示.

2 計算結(jié)果分析

2.1性能參數(shù)

霍爾推力器仿真采用COMSOL 5.2.0.166版本，仿真電腦采用Intel(R) Core(TM) i5-4570 CPU，主頻3.2 GHz，安裝內(nèi)存(RAM)16 GB，64位操作系統(tǒng).

對于1.5 kW級霍爾推力器，選取放電電壓和電流分別為368 V和3.7 A，陽極氣體標(biāo)準(zhǔn)流率為4.2 mg/s，計算初始條件為電子密度為1×1014m-3，初始平均電子能量為4 V.理論計算結(jié)果為：功率1 361.6 W、推力95.6 mN、比沖2 317 s.通過采用之前所描述的模型，可以模擬獲得的磁場、粒子數(shù)密度、電勢等主要參數(shù)的空間分布.

霍爾推力器等離子通道內(nèi)磁場線分布及通道中軸處磁場強(qiáng)度沿軸向分布見圖4.由圖可見，中軸線上磁場有正負(fù)梯度的變化，并且在軸線兩側(cè)的磁場逐漸增強(qiáng)，形成一個不均勻的鞍形磁場分布，這與ATON型霍爾推力器的凸向陽極的磁場形狀設(shè)計相符.中軸線上，磁場強(qiáng)度最大值位于磁極端面處，約為248.42 Gauss；最小值位于z=42.2 mm處，約為25.53 Gauss，使得在近陽極區(qū)磁場強(qiáng)度較弱，有利于減小陽極加速電壓損失.

霍爾推力器通道內(nèi)電子的傳導(dǎo)決定了電勢的分布，進(jìn)而影響工質(zhì)的電離和加速.從圖5中可以看出，推力器內(nèi)電勢降分布主要集中在推力器通道出口附近，在緩沖區(qū)及和緩沖區(qū)相接的區(qū)域電勢變化不明顯.而從通道內(nèi)軸線上的電勢分布可得，電勢在

z=54 mm之后逐漸下降，即可以得出放電通道加速區(qū)的長度.

2.2電離過程分析

通過COMSOL對放電通道的仿真，可記錄0～1 s 內(nèi)通道內(nèi)部電子或離子數(shù)密度的變化情況(見圖7)，即電離過程.圖8為不同時刻的推力器通道中軸線上離子數(shù)密度分布.在放電初期(t=10-6s)，放電集中于通道出口處，而隨著電子與原子發(fā)生碰撞電離，離子數(shù)密度逐漸增加.同時，電子發(fā)生霍爾漂移，向陽極方向傳導(dǎo)維持等離子體放電，因此離子數(shù)密度繼續(xù)升高，在近陽極區(qū)域由碰撞產(chǎn)生的離子數(shù)也逐漸增大.由于離子在加速段向推力器出口運(yùn)動，故在近出口處離子數(shù)密度最大.放電至10-4s后保持穩(wěn)定，此后推力器內(nèi)部離子數(shù)密度分布基本不變.目前離子數(shù)密度的仿真計算結(jié)果為：在氙氣氣量為4.2 mg/s、電壓為368 V條件下，推力器通道內(nèi)離子數(shù)密度最高為2.53×1019m-3，達(dá)到預(yù)期的電離效果.

3 優(yōu)化設(shè)計

在驗(yàn)證已有推力器主要性能參數(shù)與模型計算結(jié)果相符的基礎(chǔ)上，進(jìn)行優(yōu)化設(shè)計分析.對于給定功率與放電通道寬度的推力器，主要通過優(yōu)化磁場位形、改善磁聚焦特性等方法.

3.1磁聚焦

在ATON型霍爾推力器中，附加磁線圈的引入產(chǎn)生凸向陽極的磁場形狀，使離子流聚集到通道的中心，減少了離子壁面損失和離子束的羽流發(fā)散損失.若彎曲程度過小可能無法達(dá)到聚焦離子的效果，而彎曲程度過大則有可能出現(xiàn)過渡聚焦現(xiàn)象，即離子可能從一個壁面沿徑向遷移而運(yùn)動到另一個壁面，同樣會造成壁面腐蝕、性能下降[8].因而，磁力線的彎曲程度需要進(jìn)一步優(yōu)化設(shè)計，使得聚焦效果最好.

在保持通道中軸線處磁場強(qiáng)度幾乎不變的條件下，調(diào)整磁場磁力線與通道中心線的關(guān)系，設(shè)計3種典型的磁場位形(見圖9)，比較推力器的性能.在氙氣氣量為4.2 mg/s、電壓為368 V條件下，比較3種磁場位形下推力器通道內(nèi)的參數(shù)(見圖10)：聚焦型磁場的性能最佳，通道內(nèi)離子數(shù)密度最高為2.63×1019m-3，出口離子速度為22 371 m/s，相應(yīng)比沖為2 283 s，接近理論計算值.此外，電離區(qū)的離子分布顯示出聚焦磁場不僅可提高出口速度，還可提高電離率以及中心區(qū)域的離子密度.

因此，在霍爾推力器設(shè)計時，應(yīng)設(shè)計聚焦型磁場，有利于減小等離子體的發(fā)散和對內(nèi)外壁面的侵蝕作用，以及提高比沖.

3.2磁屏設(shè)計

ATON將推力器中心磁鐵分為兩部分設(shè)計，降低了放電通道前部的磁場強(qiáng)度，即形成靠近陽極的零磁場區(qū)，使得磁場向通道后部聚集，形成出口處的大梯度磁場區(qū).為了降低通道前半部分的磁場強(qiáng)度及調(diào)節(jié)通道內(nèi)的磁場形狀，還須進(jìn)行磁屏設(shè)計.外磁線圈產(chǎn)生的部分磁場通過與外磁屏及外圍空氣形成的局部回路而損失掉，同理，內(nèi)磁屏也消耗了部分內(nèi)磁線圈產(chǎn)生的磁場[9].在建立的模型中考慮在內(nèi)線圈靠近進(jìn)氣區(qū)域放置磁屏，采用高飽和、耐腐蝕的鐵鋁軟磁合金.

加入磁屏后放電通道中軸線上磁場強(qiáng)度分布如圖11所示.磁屏的加入使得通道前半部分的磁場強(qiáng)度明顯降低，近陽極區(qū)域的磁場強(qiáng)度降至6.83 Gauss，有利于抑制電子對陽極電極的侵蝕；同時，磁極端面處最大場強(qiáng)相對提高.此時通道內(nèi)離子數(shù)密度的分布如圖12所示，等離子體整體電離度改變較小，但受磁場影響更易于約束在放電中心區(qū)域內(nèi).相比無磁屏設(shè)計，等離子中心區(qū)域的分布比例較高，靠近壁面的數(shù)量較少.同時，磁聚焦中心區(qū)域離子受電場加速作用，出口離子速度進(jìn)一步增大，達(dá)到25 167 m/s，此時比沖為2 568 s.

4 結(jié) 論

對于霍爾推力器放電過程的二維仿真計算，本文利用COMSOL有限元分析軟件建立模型進(jìn)行仿真，可對推力器的磁場、電場以及等離子體的分布與性能參數(shù)進(jìn)行合理預(yù)測.在此基礎(chǔ)上，通過改進(jìn)磁場聚焦性能和磁屏設(shè)計，可改進(jìn)原有推力器的性能.

[1] GOEBELD M, KATZ I. Fundamentals of electric propulsion: ion and Hall thrusters[M]. John Wiley & Sons, 2008.

[2] GAMERO-CASTANO M, KATZ I. Estimation of Hall thruster erosion using HPHall[M]. Pasadena, CA: Jet Propulsion Laboratory, National Aeronautics and Space Administration, 2005.

[3] VICINI A, PASSARO A, BIAGIONI L. Hall thruster 3D plume modeling and comparison with SMART-1 flight data[J]. European Space Agency, 2004：555.

[4] HAGELAARG J M， PITCHFORD L C. Solving the boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models[J]. Plasma Sources Science and Technology,2005(14):722-733.

[5] CHRISTOU A G. Modeling of an electrical propulsion system: towards a hybrid system[R]. Royal Military College of Canada, 2015.

[6] DAVIDB. Go. Gaseous ionization and ion transport: an introduction to gas discharges[D]. Cambridge: University of Notre Dame, 2012

[7] SZABO Jr J J. Fully kinetic numerical modeling of a plasma thruster[D]. Cambridge: Massachusetts Institute of Technology, 2001.

[8] 于達(dá)仁，劉輝，丁永杰,等. 空間電推進(jìn)原理[M].哈爾濱：哈爾濱工業(yè)大學(xué)出版社，2014.

[9] 潘海林，丁鳳林，李永等. 空間推進(jìn)[M].西安：西北工業(yè)大學(xué)出版社，2016.

SimulationandOptimizationofHallEffectThruster

CONG Yuntian, CHEN Jian

(BeijingInstituteofControlEngineering，Beijing100190，China)

Based on COMSOL soft, a method is established to model the discharge processing in Hall effect thruster (HET) channel. The drift diffusion of electron and ion of the plasma is taken as the kernel of the model. Combined with the electromagnetic field, flow and impact reactions inside plasma， the characteristics of HET and distribution of parameters in the channel can be evaluated with proper selection of equations, boundary conditions and solvers. Compared with the results in theory, numerical simulation results demonstrate the effectiveness of the proposed approach， and then the magnetic field optimization design is applied.

HET； plasma； drift diffusion； finity

2016-12-08

V439

A

1674-1579(2017)05-0061-07

10.3969/j.issn.1674-1579.2017.05.010

叢云天(1993—)，男，助理工程師，研究方向?yàn)殡娡七M(jìn)技術(shù)；陳健(1969—)，男，研究員，研究方向?yàn)橥七M(jìn)技術(shù).