50 kW 級高功率霍爾推力器放電通道數(shù)值模擬研究

陳篤華，王平陽，王 尚，劉 佳

（1.上海交通大學(xué) 機械與動力工程學(xué)院，上海 200240；2.上海空間推進研究所，上海 201112）

0 引言

高功率霍爾推力器具有高比沖、大推力、長壽命以及參數(shù)可調(diào)范圍廣等優(yōu)勢，通過調(diào)節(jié)工質(zhì)流量和放電電壓可滿足空間任務(wù)在不同階段對推力器的性能需求，并大幅度節(jié)約燃料。同時，高功率霍爾推力器工作時也有著因放電通道內(nèi)產(chǎn)生較高密度等離子體導(dǎo)致的熱負荷以及溫度過高引起的磁路元件性能下降等問題。數(shù)值仿真可得到推力器詳細參數(shù)分布，為推力和比沖的計算提供放電通道內(nèi)部電勢、等離子體數(shù)密度、電子溫度等具體參數(shù)，同時也能通過計算結(jié)果分析放電通道內(nèi)等離子體生成機理和推力器熱負荷。上述工作為推力器性能分析及優(yōu)化設(shè)計提供詳細理論支撐，具有重要參考價值。

對于霍爾推力器的數(shù)值模擬研究由來已久，俄羅斯的MOROZOV 等采用一維流體動力學(xué)模型對霍爾推力器進行數(shù)值模擬分析，得到了一維條件下的粒子參數(shù)穩(wěn)態(tài)解析解和時變動態(tài)解；麻省理工學(xué)院FIFE提出了二維高精度的單元粒子法（Particle-in-Cell，PIC）計算模型，用以計算從通道內(nèi)到羽流場的霍爾推力器粒子運動過程及電子溫度、電場強度、磁感應(yīng)強度、粒子數(shù)密度分布等參數(shù)；加州理工學(xué)院MIKELLIDES 等采用二維軸對稱條件下的磁場網(wǎng)格一致性等離子體求解程序?qū)ASA-300 M 型霍爾推力器進行仿真，得到了標準工況下（放電電壓300 V，放電電流40 A）推力器通道內(nèi)及附近羽流場的電勢分布、電子溫度、電子數(shù)密度等參數(shù)；上海交通大學(xué)嚴立等模擬了SPT-70 加速通道的粒子參數(shù)分布以及粒子碰撞類型對流場分布的影響，并考慮了工作過程中生成的等離子體對加速通道壁面的相互作用所導(dǎo)致的熱效應(yīng)。相比于中低功率推力器研究進展，目前對于多工況下的高功率霍爾推力器性能仿真研究工作國內(nèi)鮮有報道。

為配合50 kW 級高功率霍爾推力器的樣機設(shè)計與實驗研究工作，本文基于PIC/蒙特卡羅碰撞模型（Monte Carlo Collision，MCC）/直接模擬蒙特卡羅碰撞模型（Direct Simulation Monte Carlo，DSMC）方法針對特定結(jié)構(gòu)的50 kW 級高功率霍爾推力器進行了二維數(shù)值模擬仿真研究，得到多種穩(wěn)定運行工況下推力器通道內(nèi)部的粒子參數(shù)分布；并依據(jù)仿真結(jié)果，分析推力器內(nèi)部工作過程以及推力器的熱損失功率情況，并計算了推力器推力及對應(yīng)比沖。此項研究工作可為改進推力器構(gòu)型、提高推力器工質(zhì)效率和開展高功率霍爾推力器地面試驗提供理論依據(jù)和技術(shù)參考。

1 數(shù)值仿真

1.1 計算基本模型

本文采用了PIC/MCC/DSMC 方法對推力器的工作過程進行了數(shù)值模擬工作：PIC 主要用于粒子運動過程跟蹤求解和通道內(nèi)自洽電場求解；MCC 主要用于求解相對速度較大的電子和中性粒子之間的碰撞過程；DSMC 主要用于求解相對速度較小的中性粒子和離子相互之間的碰撞過程。三者綜合使用以完整模擬推力器工作時的等離子體流動。

基于PIC 計算時采用了靜電求解模型，電場的微分解析通過泊松方程給出，即

式中：

Φ

為電勢；

ρ

為電荷密度；

ε

為真空介電常數(shù)。

泊松方程的求解考慮五點差分形式，為了提高計算的收斂速度，采用超松弛迭代法。計算模型霍爾推力器的陰極中置于推力器中部，如按照真實位置對陰極及電子運動軌跡進行模擬，則需要極大的計算量并難以捕捉電子的真實運動軌跡，故此選擇虛擬陰極。考慮以準中性入射模型建立虛擬陰極，在計算的每一時間步長內(nèi)通過區(qū)域邊界向等離子體入射一定量的電子以維持靠近陰極邊界區(qū)域的準中性條件。具體的處理方式為計算陰極邊界所有網(wǎng)格總的凈離子數(shù)目，計算式為

如果Δ

N

>0，則總的離子數(shù)目大于電子數(shù)目，偏離電中性，需要在邊界打入相應(yīng)的電子，以麥克斯韋分布方式打入；如果Δ

N

<0，則總電子數(shù)目大于總離子數(shù)目，不需要打入電子。

在推力器運行過程中原子速度為百米每秒量級，遠遠小于離子和電子的運動速度，可以認為計算模型中的原子為近似背景場。原子的分布情況滿足通道等離子體擴散準則，考慮在模型中電子與原子碰撞只生成一價離子，電子與原子碰撞的概率為

式中：

P

為碰撞概率；

n

為中性原子的數(shù)密度；

v

為電子和中性原子的相對速度，可近似為電子速度；

σ

為電子和中性原子碰撞截面的總和，包括彈性碰撞截面、激發(fā)碰撞截面和電離碰撞截面，與電子能量分布相關(guān)；Δ

t

為時間步長。

內(nèi)部的碰撞過程主要為電子與原子、離子之間的碰撞，忽略中性粒子相互之間的彈性碰撞。對于3 種粒子（中性粒子、氙離子、離子電荷）相互間的碰撞情況，對于中性粒子而言，碰撞集中發(fā)生在緩沖區(qū)和電離區(qū)，數(shù)密度相對較大。對于離子而言，相比原子有著較高的速度和能量，加速噴出的離子即為推力器推力來源。粒子的運動和加速通過解耦計算獲得，通過計算電場得到粒子所受電場力，通過磁場分布得到粒子所受磁場力，再進一步計算粒子的加速度和速度。

計算達到穩(wěn)定后，以粒子為基本單元統(tǒng)計其在陽極表面和放電通道內(nèi)外壁面的能量沉積所引起的推力器熱損失情況。推力器加速通道內(nèi)中性粒子能量較低，可視為與壁面處于熱平衡狀態(tài)；產(chǎn)生的電子質(zhì)量較小，在鞘層電勢影響下與壁面的碰撞過程受抑制，因此，能量沉積計算中忽略電子的熱流作用，只考慮離子在加速通道陶瓷壁面的能量沉積。等離子體沉積到加速通道壁面上的能量為

式中：

q

為等離子體沉積在壁面的總能量；Δ

q

為碰撞過程中離子動能的變化量；Δ

q

為離子內(nèi)能的變化量；

e

為電荷量；Δ

Φ

為鞘層電勢。

整體計算流程可描述為：初始化計算模型，對計算域進行網(wǎng)格劃分，將磁場數(shù)據(jù)耦合加載到計算程序，對各類物理參數(shù)初始化和粒子進行初始分布。根據(jù)工質(zhì)流量在計算的每個時間步長入射對應(yīng)數(shù)目中性氙原子，依據(jù)準中性入射模型入射對應(yīng)數(shù)目電子，根據(jù)粒子的現(xiàn)有速度計算其運動狀態(tài)。通過PIC 建立粒子間的電磁場作用模型，通過DSMC 計算原子、離子間的碰撞和動量交換、電荷交換，通過MCC 計算電子和原子間的碰撞。通過求解泊松方程得到電勢、電場分布，再依據(jù)載入的磁場數(shù)據(jù)計算粒子的加速度與速度。通過所編寫程序判定流場中的粒子凈流量小于某初設(shè)值時，則視為計算進行至推力器完全穩(wěn)定工作狀態(tài)，輸出計算結(jié)果，若不穩(wěn)定則繼續(xù)進行計算。計算至穩(wěn)定狀態(tài)后，統(tǒng)計相關(guān)參數(shù)。

1.2 計算區(qū)域及網(wǎng)格劃分

模擬計算的50 kW 級霍爾推力器為軸對稱型結(jié)構(gòu)，該推力器在周向上各類參數(shù)分布均勻，因此采用二維對稱計算模型可模擬推力器運行基本情況。數(shù)值計算時所選取的計算區(qū)域如圖1 所示。主要包括推力器的放電通道、出口及部分羽流區(qū)域，其中左側(cè)邊界為推力器陽極，以推力器中心軸線為對稱軸。根據(jù)推力器的固有尺寸，放電通道的長度

L

=95 mm，其中內(nèi)壁面處于

R

=165 mm處，外壁面處于

R

=235 mm 處，在加速通道外部選取半徑為400 mm、軸向長度95 mm 的圓柱形羽流場區(qū)域。采用結(jié)構(gòu)化均勻網(wǎng)格，網(wǎng)格單元長度選取為2.5 mm。

圖1 推力器結(jié)構(gòu)及流場仿真區(qū)域Fig.1 Construction and simulation zone of the thruster

1.3 加速計算說明

工質(zhì)氙原子和電子的質(zhì)量相差5 個數(shù)量級，在相同的能量條件下，氙原子通過加速通道的時間約為電子的500 倍，若按照電子的時間步長推進離子，將花費大量時間和計算資源達到穩(wěn)定，故此計算中采用降低氙原子質(zhì)量和提高介電常數(shù)的辦法以提高計算速度。具體表現(xiàn)為將氙原子質(zhì)量降低2 500 倍，則運動速度相應(yīng)提高50 倍；將介電常數(shù)提高100 倍，德拜長度提高10 倍，總體上降低時間步長，提高空間步長。

1.4 推力及比沖計算

推力器工作時放電通道內(nèi)生成的離子在電場作用下加速噴出產(chǎn)生推力，在單位時間內(nèi)推力器獲得推力等于從通道內(nèi)噴出離子的總動量。對推力器標準工況進行數(shù)值模擬，統(tǒng)計推力器穩(wěn)定工作時通道出口的離子數(shù)目及軸向速度。根據(jù)已知的氙離子質(zhì)量（Xe），可得單位時間內(nèi)的離子總動量，進而計算推力。再根據(jù)比沖定義，即可計算比沖為

式中：

I

為比沖；

F

為推力；

m

為單位時間內(nèi)的工質(zhì)流量；

g

為重力加速度。

2 標準工況數(shù)值模擬

2.1 磁場分布

推力器磁場考慮為靜態(tài)磁場，由ANSYSMAXWELL 商用軟件計算得到分布結(jié)果，以輸入文件形式將磁場數(shù)據(jù)加載至模型計算程序。計算所得放電通道中軸線上磁感應(yīng)強度分布如圖2 所示，從陽極出口到羽流近場，整體呈先上升再下降趨勢，最大值220.12×10T 出現(xiàn)在通道出口附近，符合霍爾推力器內(nèi)部磁場設(shè)計準則及分布規(guī)律。

圖2 通道軸線磁感應(yīng)強度分布Fig.2 Distribution of the magnetic flux density along the channel axis

2.2 標準工況內(nèi)部參數(shù)分布及性能參數(shù)

對標準工況下的50 kW 級推力器工作情況仿真分析，得到了放電通道內(nèi)的粒子數(shù)密度分布、電子溫度、離子通道軸向速度等結(jié)果。標準工況的輸入?yún)?shù)包括功率50 kW，放電電壓500 V，工質(zhì)流量為86.4 mg/s。離子數(shù)密度分布如圖3 所示，最大值出現(xiàn)在通道中間區(qū)域，最大值為1.635×10m，隨著通道內(nèi)電磁場對氙離子加速，數(shù)密度迅速降低，在通道出口位置的氙離子數(shù)密度約為1.679×10m。沿通道軸向離子數(shù)密度先增加后降低，分布情況和HOFER 等的結(jié)果相似，推力器放電通道內(nèi)分為明顯的緩沖區(qū)、電離區(qū)和加速區(qū)，滿足放電通道粒子分布特征。通道內(nèi)的速度分布如圖4所示，離子速度在通道后半段迅速增加，并在電場作用下保持加速效果至通道出口，最大出口速度達到了28 150 m/s。

圖3 標準工況下離子數(shù)密度分布Fig.3 Distribution of the ion number density under the standard condition

圖4 標準工況下離子軸向速度分布Fig.4 Distribution of the ion axial velocity under the standard condition

統(tǒng)計可知，在考慮壁面典型平均溫度600 K 條件下，標準工況下的推力器放電通道內(nèi)壁面的熱損失功率為3 068.97 W，外壁面的熱損失功率為4 191.84 W，陽極表面（溫度750 K）的熱損失功率為710.91 W，推力器整體熱損失功率占總功率比重為15.94%，與文獻［10］的計算結(jié)果相比偏小，推力器在穩(wěn)定運行時存在一定的熱負荷問題。

利用標準工況計算得到的流場參數(shù)，結(jié)合1.4 節(jié)的方法和式（4）分別獲得了推力為2.2 N，比沖為2 598 s，與文獻［16］中同類50 kW 級霍爾推力器實驗結(jié)果比較，推力和比沖的誤差分別為5.18%和3.35%。

3 多工況模擬

推力器實際運行過程中根據(jù)空間任務(wù)需要，需在多種工作模式下穩(wěn)定運行。控制推力器輸出的主要方式包括調(diào)節(jié)工作電壓和工質(zhì)流量，針對上述工作情況，本章選擇了標準工況工質(zhì)流量、不同放電電壓和標準工況放電電壓、不同工質(zhì)流量的輸入條件，對多種工況下運行的推力器進行仿真研究。

3.1 不同放電電壓

文中計算模型推力器在設(shè)計最佳運行條件下的放電電壓為500 V，考慮實際需求偏離標準放電電壓正負20%，選取放電通道中軸線上的計算結(jié)果輸出。

電子溫度隨放電電壓的變化如圖5 所示，在不同的放電工作電壓下，電子溫度呈現(xiàn)相同的變化趨勢，分布趨勢與文獻［17］中結(jié)果一致，電子在磁場中繞磁力線作拉莫爾回旋運動，磁場強度決定了放電通道約束電子能力的強弱。電子溫度從陽極到通道口的軸向分布趨勢和磁場分布趨勢近乎相同，可以說明磁場對電子有明顯束縛作用。不同差異點在于放電電壓不同導(dǎo)致的電子溫度最大值不同，偏差幅度與工作電壓的變化幅度相近。

圖5 不同放電電壓下電子溫度分布Fig.5 Distribution of the electron temperature at different discharge voltages

離子數(shù)密度的變化如圖6 所示，不同電壓下的數(shù)密度分布與標準工況下分布趨勢基本一致，數(shù)值上放電電壓越高，電子和離子的運動運動速度越大，碰撞和電離過程更劇烈，在通道中部的最大離子數(shù)密度越高。離子運動速度變化如圖7 所示，工作電壓的升高使離子在加速階段獲得更大的能量，電壓越高，出口處離子運動速度越高。

圖6 不同放電電壓下離子數(shù)密度分布Fig.6 Distribution of ion number density at different discharge voltages

圖7 不同放電電壓下離子軸向速度分布Fig.7 Distribution of the ion axial velocity at different discharge voltages

3.2 不同質(zhì)量流量

推力器設(shè)計標準工況下質(zhì)量流量為86.4 mg/s，考慮實際需求條件偏離標準流量正負20%，統(tǒng)計結(jié)果顯示，在69.12 mg/s 到103.68 mg/s 的輸入流量下，推力器保持近乎相同的流場參數(shù)分布趨勢。電子溫度和離子數(shù)密度隨質(zhì)量流量的變化如圖8～圖9 所示，變化趨勢相同，數(shù)值變化范圍超過了流量變化的20%。在低流量時粒子數(shù)密度降低，粒子相互間的碰撞作用減弱，電子溫度和生成離子的數(shù)密度進一步降低；質(zhì)量流量越高，中間粒子的相互作用過程越劇烈，電子溫度和離子數(shù)密度越高。離子運動速度在不同質(zhì)量流量條件下分布近乎相同，如圖10 所示，說明在此范圍內(nèi)的輸入流量變化幾乎不會影響通道內(nèi)的粒子加速過程。

圖8 不同質(zhì)量流量下電子溫度分布Fig.8 Distribution of the electron temperature at different mass flow rates

圖9 不同質(zhì)量流量下離子數(shù)密度分布Fig.9 Distribution of the ion number density at different mass flow rates

圖10 不同質(zhì)量流量下離子軸向速度分布Fig.10 Distribution of axis ion velocity at different mass flow rates

3.3 推力及比沖

推力器工作時推力和比沖直接體現(xiàn)了推力器工作性能，通過前面對放電通道內(nèi)部過程的數(shù)值模擬，可得到出口所有粒子的速度和質(zhì)量流量，進而求解推力和比沖。多工況下的推力和比沖如圖11～12所示。

圖11 多工況下推力Fig.11 Thrust under multiple conditions

圖12 多工況下比沖Fig.12 Specific impulse under multiple conditions

在相同流量下隨著電壓升高，出口離子速度越大，推力器比沖和推力越大；在相同工作電壓下考慮工質(zhì)流量變化，出口離子速度基本相同，比沖幾乎不變，推力隨流量增大而增大。

4 結(jié)束語

基于PIC/MCC/DSMC 混合數(shù)值模擬方法，計算了標準工況和其他工況下放電通道內(nèi)部的等離子體參數(shù)，為后續(xù)50 kW 級高功率霍爾推力器的樣機設(shè)計與地面試驗提供了理論支撐，計算結(jié)果表明：

1）標準工況下，仿真結(jié)果與文獻［16］中同類推力器實驗結(jié)果比較，推力和比沖誤差分別為5.18%和3.35%；通道內(nèi)的離子數(shù)密度最大達到1.635×10m，隨著通道內(nèi)自洽電場對離子的加速作用，數(shù)密度在出口位置下降至1.679×10m，通道內(nèi)分為明顯的緩沖區(qū)、電離區(qū)和加速區(qū)。

2）標準工況下，陽極表面平均溫度為750 K 和放電通道內(nèi)外壁面在典型平均溫度600 K 條件下的熱損失占總功率的15.94 %，推力器存在一定的熱負荷問題。

3）在工作電壓400～600 V、質(zhì)量流量69.12～103.68 mg/s 范圍內(nèi)考察了流場的分布情況。結(jié)果顯示：工作電壓提升會增加電子溫度和離子出口速度；調(diào)節(jié)質(zhì)量流量，電子溫度和離子數(shù)密度分布趨勢相同，離子速度基本不變。在此類工況條件下增大工作電壓，推力比沖均隨之增大；增大質(zhì)量流量，推力隨之增大，比沖幾乎不變。