霍爾推力器工作性能數值模擬研究

龍建飛，張天平，孫明明，高 俊，賈連軍

(蘭州物理研究所 真空低溫技術與物理重點試驗室，蘭州 730000)

?

霍爾推力器工作性能數值模擬研究

龍建飛，張天平，孫明明，高 俊，賈連軍

(蘭州物理研究所 真空低溫技術與物理重點試驗室，蘭州 730000)

為了模擬霍爾推力器推力、比沖等工作性能，采用粒子網格單元與蒙特卡洛相結合(PIC/MCC)方法，建立了霍爾推力器二維軸對稱模型。模型中電子和離子均采用粒子描述，中性原子為背景氣體，自洽電勢通過求解泊松方程獲得。跟蹤推力器出口處離開的離子數量、軸向速度等信息，通過統計計算得到推力器的推力、比沖。以LHT100推力器為研究對象，針對不同的工作參數(陽極流量在4.6～5.4 mg/s，陽極電壓在280～320 V之間)共進行了9種工況數值模擬，并進行試驗對比驗證，模擬結果與試驗測試結果均較好吻合，最大誤差小于10%。

霍爾推力器；工作性能；蒙特卡洛；粒子單元法；數值仿真

0 引言

霍爾推力器具有高可靠、高功推比等特點，作為空間動力裝置而廣泛應用于衛星的位置保持和姿態控制等空間任務中[1-2]。推力和比沖作為霍爾推力器的主要性能而備受關注，是衡量該推力器能否空間應用的重要指標。對霍爾推力器的推力和比沖等性能進行仿真研究，可明晰該推力器工作物理過程，以及為性能優化提供理論支撐，具有重要意義。

霍爾推力器的推力及比沖等工作性能的理論計算主要是通過經驗模型[3-4]獲得，建立的經驗模型一般是推力器工作參數和結構參數的函數。在經驗模型建立過程中，涉及電勢分布、等離子體密度、工質利用率等相關敏感參數，而這些敏感參數很難通過試驗精確測出。因此，均是通過經驗估值給定，這便使得經驗模型存在一定偏差。同時，推力器在不同工況下，敏感參數的不一致性，甚至會隨工作參數存在較大變化，這便加大了經驗模型誤差，進而使得模型不能滿足廣泛工作條件。關于霍爾推力器放電室數值模擬研究已經開展了大量研究，Szabo[5]建立了放電室的全粒子模型，Fife[6]和Hofer等[7]建立了混合模型等，均對霍爾推力器放電室工作過程進行了數值模擬研究，得到了放電室內電勢分布、等離子體密度和能量分布等重要信息。國內劉輝等[8]采用PIC/MCC方法，對霍爾推力器放電室內電子運動行為進行了數值模擬研究。而關于霍爾推力器工作性能的數值仿真，特別是廣泛工作條件下的仿真計算，目前還鮮有文章報道。

本文將采用全粒子PIC/MCC方法，建立霍爾推力器放電室二維軸對稱模型，對放電室等離子體進行仿真，跟蹤推力器出口處離開的離子信息(主要是離子數量和軸向速度)，通過建立的性能模型計算，進而模擬出霍爾推力器的推力和比沖性能，針對廣泛的工作參數(陽極流量在4.6～5.4 mg/s，陽極電壓在280～320 V之間)共進行了多工況數值模擬，并開展試驗驗證研究。

1 數值模型建立

1.1 模擬區域及網格劃分

霍爾推力器具有軸對稱的幾何結構，且推力器的參數在周向的變化相對較小。因此，可采用2維柱坐標進行模擬，具體包括軸向(z)和徑向(r)的二維空間和3維速度(Vz，Vr，Vθ)，θ表示軸向方向。圖1為所選取的計算區域，計算區域主要為推力器放電室加速通道，左側邊界為推力器陽極。根據LHT100推力器尺寸，加速通道長度取26 mm，內壁面位于r=35 mm位置處，外壁面位于r=50 mm位置處。采用正交等距劃分網格，長度為0.5 mm。因此，通道網格劃分為52×30的等距網格。

圖1 計算區域選取示意圖Fig.1 Schematic of calculation area

1.2 電磁場及粒子運動

在采用PIC粒子模擬方法建模時，均采用靜電模型。因此，電場采用泊松方程求解[9]。

(1)

式中Φ為電勢；ε0為真空介電常數；ni為離子密度；ne為電子密度。

泊松方程求解采用有限差分法，所使用的格式為中心差分五點格式。為了提高電場計算收斂的速度，文中采用逐次超松弛(SOR)方法。

采用三維磁場測試儀，對LHT100推力器的真實工況磁場進行測試，并將結果進行導入。磁場在整個計算過程一直保持不變。由于在霍爾推力器中，帶電粒子自身產生的磁場較小，對計算結果影響不大。因此，為了簡化計算，忽略帶電粒子自身產生的磁場，即只考慮線圈產生的磁場。

霍爾推力器的陰極位于推力器外部，模擬真實的陰極需要較大計算量，并增加了程序設計的難度。因此，本文沒有模擬真實的陰極，而是在推力器出口采用虛擬的陰極向通道內噴射電子。每個時間步長噴射進入通道的電子數量為

(2)

式中Id為放電電流；Δt為時間步長；w為仿真粒子的權重；q為電子電荷；Δnc為上一個時間步長從推力器出口處離開的粒子數。

入射電子滿足半麥克斯韋分布，其能量為15 eV(與陰極觸持極電壓保持一致)。

由于原子的速度遠小于離子和電子，因此在模型中，認為原子的分布是背景場。原子數密度分布近似地可表示為

(3)

式中na為原子密度；m為陽極流量；s為推力器出口端面積；v0為原子速度，原子速度在放電室內處處相等。

在模型中，模擬粒子的加速采用牛頓-洛倫茲力求解[10]

(4)

模型中認為粒子之間只有電子一原子碰撞才能生成離子，進一步認為只產生單電荷離子。電子和原子的碰撞概率表示為[11]

p=1-exp(-naveσtΔt)

(5)

式中ve為電子速度；σt為電子與原子間的碰撞截面，其值與電子能量有關。

由于在每個時間步長里都要計算所有電子的動能及碰撞概率，計算量很大。因此，本文采取Birdsall和Vahedi提出的空碰撞的方法[12]。

1.3 模擬計算簡化

在目前的計算機硬件條件下，使用普通計算機直接模擬離子與電子的運動是非常困難的，因為離子比電子重，其質量比電子質量高3個數量級，這將導致離子的收斂很慢，計算時間過長。為了能夠同時模擬電子和離子的運動，引入人工質量比[8]，將離子的質量人為的減小2 500倍數，以達到加速收斂的作用。此外，采用增大介電常數[8]1 600倍減少網格數量，進一步減小計算收斂時間。

1.4 工作性能模型

霍爾推力器是通過電能轉化將工質氣體(氙氣)電離，并通過自洽電場將離子加速噴出，進而產生推力。這一過程可采用動量守恒方程進行描述，推力產生示意圖見圖2。圖2中，Δv、m分別為推力器在Δt時間后推力器所獲得的速度增量和噴出的粒子質量；u為噴出粒子所獲得的平均軸向速度(相對于推力器)。

圖2 霍爾推力器推力產生示意圖Fig.2 Schematic diagram of thrust form in Hall thruster

假設推力器受到的外力的合力為∑F，根據動量定理，該過程滿足方程：

(M-m)(v+Δv)+m(v-u)-Mv=(∑F)Δt

(6)

對式(6)做進一步分析：(1)忽略上式中mΔv高階項；(2)推力器在空間工作過程中，只受到噴出束流反作用力的作用，即外力為零；(3)推力器出口噴出的粒子主要考慮離子。由于電子質量遠小于離子，同時原子不受電場加速(速度過小)。因此，推力主要由噴出的離子產生。根據作用力與反作用力，則可得到推力器受到的推力，以及推力器獲得的比沖。

(7)

從式(7)中可得出，推力器獲得的推力等于單位時間內噴出束流離子的總動量。而比沖為噴出離子的平均軸向速度。為了獲得推力和比沖性能，通過穩態工況下的數值仿真，統計出口處離開離子的數量和軸向速度，由于單個離子質量已知(Xe+)，則可求出單位時間噴出的離子總動量，對應為霍爾推力器的推力。根據比沖定義公式，可進一步求出比沖。

2 放電室等離子體仿真

對額定工況(陽極流量：5.0 mg/s，陽極電壓為300 V)下LHT100推力器工作過程進行了仿真，分別得到了放電室內的電勢分布、等離子體密度分布和離子軸向速度分布等重要信息。

圖3為推力器放電室內電勢分布。仿真結果顯示，電勢沿徑向具有較好的一致性，沿軸向電勢降主要集中在放電室末端，離陽極約20 mm以后開始大幅下降，仿真結果與Fife等[6]研究結果基本一致。放電室內電勢分布主要由電子的電導率分布決定的。為了保證有足夠的電子到達陽極以維持放電的穩定性，近陽極區域的電導率相對需要較大。而在電離區，為了保證電子具有足夠的能量向陽極遷移以維持放電，使得電離區電導率要根據電子溫度的變化向陽極方向逐漸增大。在推力器出口末端很窄一段區域內，為了讓離子獲得較高的加速電壓，同時要使電子從陰極穿越加速區到達電離區的過程中得以“升溫”，確保電場提供給電子的能量足以彌補電子遷移過程中因碰撞而造成的損失，即該區域的電導率較小。因此，電導率的特征分布就決定了電勢分布。

圖3 等離子體電勢分布Fig.3 Plasma potential distribution in discharge room

圖4為放電室離子密度分布。仿真結果顯示，通道內部最大離子數密度量級為1018mm-3，與Fife等[6]研究結果基本一致。離子密度分布與放電室內部的磁場分布是密切相關的。放電室徑向磁場沿軸向呈“正梯度”磁場分布，即出口處的磁感應強最強，而陽極附近最小。磁場強度決定了放電室內對電子約束能力，也間接地決定了原子的電離區域。仿真結果顯示，等離子體主要分布在中下游，對應于放電室的電離區，與放電室的特征區域[8]是相匹配的。

圖4 放電室離子密度分布Fig.4 Pion density distribution in discharge room

圖5為放電室內離子速度分布。從圖5可看出，離子軸向速度沿軸向逐漸增大，在出口處有最大軸向速度，約16 000 m/s，與文獻[6]基本一致。離子的軸向速度分布與放電室內加速電場強度分布有密切關系，而電場分布決定于放電室的電勢分布。通過前面電勢分布可看出，在近陽極區電勢分布變化不大，因而電場強度很小。在電離區特別是加速區，電勢降沿軸向明顯增大，即加速電場增大，使得離子獲得速度增量增大，即在出口處獲得最大速度。

圖5 放電室離子軸向速度分布Fig.5 Ion axial velocity distribution in discharge room

3 性能仿真及試驗驗證

3.1 仿真工況及測試方法

霍爾推力器的推力和比沖等性能主要受陽極流量和陽極電壓等工作參數的影響。在空間應用中，面臨復雜的工作環境，工作參數一般具有額定工作點±5%的振蕩偏差。因此，為了系統驗證該模型的準確性，驗證條件盡量涵蓋了空間應用中各種變化條件。同時考慮計算量，本文共選取9中工況進行推力、比沖進行仿真和試驗驗證研究，具體工況見表1。表1中，*表示額定工況。

霍爾推力器地面試驗中，推力采用微小推力測試方法獲得，比沖根據試驗測試推力除以試驗時流量得到。微小推力測試裝置示意圖及推力器測試工作情況見圖6。

表1 工作參數選擇Table1 The list of working parameter

3.2 推力驗證

圖7為推力性能的仿真與試驗驗證，其中共仿真了9種工況。從推力性能對比結果來看，仿真結果相比試驗結果偏低，在第9種工況下有最大誤差約9%(仿真結果為79.6 mN，測試值87.2 mN)。分析認為，主要存在3點原因：(1)仿真推力計算中，只是統計了離子，忽略了電子和原子的微推力效應，特別是原子經過放電室加熱后具有一定的速度，從推力器出口流出也會產生一定的推力；(2)研究表明，出口到陰極之間仍存在電勢差[13]，這部分電勢差仍將對離子加速，對應仿真推力比實際測試值偏小；(3)推力測試試驗中本身存在誤差，特別是本地壓強的存在(約1×10-3Pa)，會導致推力測試偏大，而仿真中并未考慮本底壓強因素。

(a)推力測試裝置

(b)推力測試過程

3.3 比沖驗證

圖8為比沖性能的仿真與試驗驗證，其中共仿真了9種工況。從比沖性能對比結果來看，在陽極流量變化情況下(前5種工況)，仿真結果比較平穩，而測試結果隨陽極流量增大而較快增長，以致兩者差異增大。而在陽極電壓變化情況下，仿真結果與測試結果之間的差異較為穩定。其中工況5情況下，兩者有最大誤差約7.5%(仿真結果為1 587 s，測試結果為1 715 s)。分析認為，本模型中忽視了二價離子的存在，而這部分離子可獲得較高速度，使得仿真出來的比沖比測試值偏低。同時，從出口到陰極之間仍占有一定的電勢差，這其中有一部分電勢降仍對出口后離子加速。所以，根據出口離子速度統計得到的比沖要比實際比沖偏小。

(a)推力隨陽極流量變化關系 (b)推力隨陽極電壓變化關系

(a)比沖隨陽極流量變化關系 (b)比沖隨陽極電壓變化關系

4 結論

(1)建立了霍爾推力器新的性能模型。模型中，推力、比沖是等離子體微觀參數的相關函數。

(2)采用(PIC/MCC)數值模擬方法跟蹤放電室內離子的運動軌跡，通過統計得到了離開離子的速度和數量等微觀信息，進而得到了推力和比沖等性能。

(3)對放電室等離子體進行了數值仿真研究，得到了放電室電勢分布、等離子體密度分布及離子軸向速度分布等信息，仿真結果與國外文獻能較好吻合。

(4)在陽極流量4.6～5.4 mg/s、陽極電壓280～320 V范圍內，共仿真計算了9種工況下推力器推力、比沖性能。將仿真結果與試驗測試值進行對比驗證，推力最大誤差9%，比沖最大誤差7.5%。

[2] 張天平. 國外離子和霍爾電推進技術最新進展[J]. 真空與低溫，2006, 12(4)：187-190.

[2] 趙杰，唐德禮，耿少飛. 圓柱形霍爾推力器內等離子體數值模擬[J]. 固體火箭技術，2010, 33(1)：54-58.

[3] Kim V. Main physical features and processes determining the performance of stationary plasma thrusters[J]. Journal of Propulsion and Power, 1998, 32(14)：5736-5743.

[4] Brophy J, Barnett J. Performance of the stationary plasma thruster: SPT-100[R]. AIAA 92-3155.

[5] Szabo J. Fully kinetic numerical modeling of a plasma thruster[D]. Massachusetts: Massachusetts Institute of Technology, 2001.

[6] Fife J. Hybrid-PIC modeling and electrostatic probe survey of hall thrusters[D]. Massachusetts: Massachusetts Institute of Technology, 1998.

[7] Hofer R R, Mikellides I G, Katz I．Wall sheath and electron mobility modeling in Hybrid-PIC Hall thrusters simulations[R]. AIAA 2007-5267.

[8] 劉輝. 霍爾推力器電子運動行為的數值模擬[D].哈爾濱: 哈爾濱工業大學，2009.

[9] 賀碧蛟，張建華，蔡國彪. 穩態等離子體推進器羽流場數值模擬[J]. 北京航空航天大學學報，2005, 31(9): 145-148.

[10] 孫安邦，毛根旺，陳茂林, 等. 霍爾效應推力器羽流的PIC/MCC模擬[J]. 固體火箭技術，2009, 32(6): 638-640.

[11] Garrigues L, Bareilles J, Boeuf J. Modeling of the plasma jet a stationary thruster[J]. Journal of Applied Physics, 2002, 91(12): 9521-9524.

[12] Birdsall C K，Langdon A B．Plasma physics via computer simulation[M]．New York：McGraw Hill Higher Education，1984.

[13] Yu D , Li Y, Song S. Ion sputtering erosion of channel wall corners in hall thrusters[J]. Journal of Physics D: Applied Physics, 2006, 39(4): 2205-2211.

(編輯：崔賢彬)

Numerical simulation of performance in Hall thruster

LONG Jian-fei, ZHANG Tian-ping, SUN Ming-ming, GAO Jun, JIA Lian-jun

(Science and Technology on Vacuum Cryogenics Technology and Physics Laboratory, Lanzhou Institute of Physics, Lanzhou 730000, China)

To simulate the thrust and specific impulse of a Hall thruster, a two-dimensional axisymmetric model is developed for the discharge channel of the thruster, the model is based on a particle description of electrons and ionic components, the neutral atoms are set as background, and the self-consistent electric potential is calculated from Poisson equation. In the process of tracking, the number and energy of ions ejecting from exit plane were recorded, then the thrust and specific impulse were obtained by statistical calculation. Treat LHT100 thruster as research object, according to different working parameters (anode flow rate is 4.6～5.4 mg/s, anode voltage is 280～320 V), numerical simulation were conducted in 9 kinds of conditions, respectively. By experiment validation, the simulation results and experiment results are consistent, the maximum error is no more than 10%.

Hall thruster；performance；Monte Carlo；particle-in-cell；numerical simulation

2014-08-25；

：2014-12-15。

真空低溫技術與物理重點實驗室基金(9140c550206130c5503)。

龍建飛(1984—)，男，博士，研究領域為空間電推進技術。E-mail:ljf510@163.com

V493

A

1006-2793(2015)04-0514-05

10.7673/j.issn.1006-2793.2015.04.012