叢云天，陳 健

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

*總裝預研基金資助項目(9140A20050315HT05002).

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

叢云天，陳 健

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

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

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

0 引 言

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

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

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

1 物理數學模型

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

1.1系統(tǒng)方程

1.1.1 電子擴散

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

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

(1)

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

(2)

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

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

存在磁場情況下，真實的電子遷移率的表達式不能用一個簡單的形式來表達，而應該采用張量形式的電子遷移率

(3)

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

電子的擴散系數、能量遷移率、能量擴散系數分別為

(4)

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

1.1.2 重物質擴散

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

(5)

式中，Jk為擴散通量矢量，Rk為k種離子的比率表達式，u是中性粒子速度矢量，ρ是混合成分的密度，wk是第k種離子的質量分數.

擴散通量矢量Jk為

Jk=ρwkVk

(6)

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

(7)

(8)

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

(9)

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

1.1.3 磁場

J=σV×B+Je

(10)

(11)

其中，A為磁勢矢量，M為磁化強度矢量

1.1.4 N-S方程

(12)

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

1.2化學反應

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

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

1.3邊界條件

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

1.3.1 電子壁邊界條件

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

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

(13)

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

(14)

1.3.2 層流邊界條件

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

1.4物理模型

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

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

2 計算結果分析

2.1性能參數

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

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

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

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

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

2.2電離過程分析

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

3 優(yōu)化設計

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

3.1磁聚焦

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

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

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

3.2磁屏設計

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

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

4 結 論

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

[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] 于達仁，劉輝，丁永杰,等. 空間電推進原理[M].哈爾濱：哈爾濱工業(yè)大學出版社，2014.

[9] 潘海林，丁鳳林，李永等. 空間推進[M].西安：西北工業(yè)大學出版社，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—)，男，助理工程師，研究方向為電推進技術；陳健(1969—)，男，研究員，研究方向為推進技術.