陳篤華，王平陽(yáng)，王 尚，劉 佳

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

0 引言

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

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

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

1 數(shù)值仿真

1.1 計(jì)算基本模型

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

基于PIC 計(jì)算時(shí)采用了靜電求解模型，電場(chǎng)的微分解析通過(guò)泊松方程給出，即

式中：

Φ

為電勢(shì)；

ρ

為電荷密度；

ε

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

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

如果Δ

N

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

N

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

在推力器運(yùn)行過(guò)程中原子速度為百米每秒量級(jí)，遠(yuǎn)遠(yuǎn)小于離子和電子的運(yùn)動(dòng)速度，可以認(rèn)為計(jì)算模型中的原子為近似背景場(chǎng)。原子的分布情況滿足通道等離子體擴(kuò)散準(zhǔn)則，考慮在模型中電子與原子碰撞只生成一價(jià)離子，電子與原子碰撞的概率為

式中：

P

為碰撞概率；

n

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

v

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

σ

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

t

為時(shí)間步長(zhǎng)。

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

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

式中：

q

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

q

為碰撞過(guò)程中離子動(dòng)能的變化量；Δ

q

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

e

為電荷量；Δ

Φ

為鞘層電勢(shì)。

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

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

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

L

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

R

=165 mm處，外壁面處于

R

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

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

1.3 加速計(jì)算說(shuō)明

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

1.4 推力及比沖計(jì)算

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

式中：

I

為比沖；

F

為推力；

m

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

g

為重力加速度。

2 標(biāo)準(zhǔn)工況數(shù)值模擬

2.1 磁場(chǎng)分布

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

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

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

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

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

圖4 標(biāo)準(zhǔn)工況下離子軸向速度分布Fig.4 Distribution of the ion axial velocity under the standard condition

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

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

3 多工況模擬

推力器實(shí)際運(yùn)行過(guò)程中根據(jù)空間任務(wù)需要，需在多種工作模式下穩(wěn)定運(yùn)行?？刂仆屏ζ鬏敵龅闹饕绞桨ㄕ{(diào)節(jié)工作電壓和工質(zhì)流量，針對(duì)上述工作情況，本章選擇了標(biāo)準(zhǔn)工況工質(zhì)流量、不同放電電壓和標(biāo)準(zhǔn)工況放電電壓、不同工質(zhì)流量的輸入條件，對(duì)多種工況下運(yùn)行的推力器進(jìn)行仿真研究。

3.1 不同放電電壓

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

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

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

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

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

圖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 推力及比沖

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

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

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

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

4 結(jié)束語(yǔ)

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

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

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

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