碘工質(zhì)霍爾推力器內(nèi)部過(guò)程數(shù)值模擬

王 尚，梁子軒，王平陽(yáng)，徐宗琦，杭觀榮

（1.上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，上海 200240；2.上?？臻g推進(jìn)研究所上?？臻g發(fā)動(dòng)機(jī)工程技術(shù)研究中心，上海 201112）

0 引言

霍爾推力器具有比沖高、推力精確可調(diào)、推力功率比大、壽命長(zhǎng)等優(yōu)勢(shì)，能夠給航天器帶來(lái)質(zhì)量、壽命、經(jīng)濟(jì)等增益，作為航天器的推進(jìn)裝置得到了廣泛應(yīng)用。隨著霍爾推進(jìn)技術(shù)的不斷發(fā)展，碘作為氙的替代工質(zhì)，因其自身升華溫度低、電離能低、相對(duì)原子質(zhì)量高、儲(chǔ)量豐富、價(jià)格低廉等特點(diǎn)，得到了廣泛研究。由于碘電離后，等離子體組分不同于氙，導(dǎo)致電離區(qū)位置和加速效應(yīng)均不同，為提高推力器的性能，需重點(diǎn)研究其內(nèi)部電離和加速過(guò)程。通過(guò)數(shù)值模擬方法可以得到推力器內(nèi)部等離子特性的詳細(xì)參數(shù)分布，得到放電通道內(nèi)部電勢(shì)、等離子體數(shù)密度、電子溫度等參數(shù)，以及推力器的性能參數(shù)，進(jìn)而深入分析放電通道內(nèi)部等離子體的生成機(jī)理。上述工作為推力器性能和壽命的優(yōu)化提供了詳細(xì)的理論支撐，具有重要參考價(jià)值。

對(duì)于碘工質(zhì)推力器放電過(guò)程的數(shù)值模擬，國(guó)內(nèi)外學(xué)者均開展了相關(guān)研究。法國(guó)的GRONDEIN等采用流體動(dòng)力學(xué)模型對(duì)碘工質(zhì)離子推力器進(jìn)行了數(shù)值模擬，獲取了離子、電子數(shù)密度和電子溫度等特征參數(shù)，將碘和氙工質(zhì)情況下的推力、比沖和效率進(jìn)行了比較；NASA 的MARIA采用混合單元粒子法模型對(duì)200 W 霍爾推力器進(jìn)行了模擬，對(duì)比了碘和氙工質(zhì)情況下羽流區(qū)的中性氣體和離子的數(shù)密度分布；法國(guó)的LUCKEN 等針對(duì)離子推力器進(jìn)行數(shù)值模擬，研究了沿推力器軸向方向的磁場(chǎng)強(qiáng)度變化時(shí)，碘工質(zhì)和氙工質(zhì)情況下推力器的推力和效率變化；哈爾濱工業(yè)大學(xué)的牛翔等采用一維流體動(dòng)力學(xué)模型對(duì)碘工質(zhì)會(huì)切場(chǎng)推力器進(jìn)行數(shù)值模擬，研究了碘工質(zhì)會(huì)切場(chǎng)推力器內(nèi)部的等離子體參數(shù)分布。相較于上述研究成果，目前對(duì)百瓦級(jí)碘工質(zhì)霍爾推力器放電通道內(nèi)部過(guò)程的數(shù)值模擬研究工作還鮮有報(bào)道。

為了配合200 W 級(jí)碘工質(zhì)霍爾推力器的設(shè)計(jì)優(yōu)化和實(shí)驗(yàn)研究工作，本文基于單元粒子法/直接模擬蒙特卡羅法/蒙特卡羅碰撞模型（Particel-in-Cell/Direct-Simulation-Monte-Carlo/Monte-Carlo-Collision，PIC/DSMC/MCC）混合方法，對(duì)一定結(jié)構(gòu)的200 W 級(jí)碘工質(zhì)霍爾推力器進(jìn)行了二維數(shù)值模擬仿真研究，得到了設(shè)計(jì)工況下推力器放電通道內(nèi)部等離子體的各項(xiàng)參數(shù)分布；根據(jù)數(shù)值模擬的結(jié)果分析推力器的內(nèi)部過(guò)程，和氙霍爾的等離子體參數(shù)分布進(jìn)行對(duì)比。此項(xiàng)研究工作可為百瓦級(jí)碘工質(zhì)霍爾推力器的性能、壽命優(yōu)化、地面實(shí)驗(yàn)的開展提供理論依據(jù)和技術(shù)參考。

1 數(shù)值仿真

1.1 計(jì)算基本模型

采用PIC/DSMC/MCC 方法對(duì)推力器的內(nèi)部過(guò)程進(jìn)行數(shù)值仿真，PIC 用于求解自洽電場(chǎng)和計(jì)算帶電粒子的加速度；DSMC 用于求解相對(duì)速度較小的中性分子-離子之間、中性分子-中性分子和離子-離子之間的碰撞，主要處理碘原子、碘離子之間的動(dòng)量交換和電荷交換碰撞；MCC 用于求解電子參與的碰撞過(guò)程，主要處理電子和碘分子、碘原子之間的碰撞，三者結(jié)合使用完成粒子間的碰撞輸運(yùn)和運(yùn)動(dòng)過(guò)程。文獻(xiàn)［10］已對(duì)此混合方法的可行性和正確性進(jìn)行了驗(yàn)證，在此基礎(chǔ)上繼續(xù)開展工作。

基于PIC 計(jì)算時(shí)采用了靜電求解模型，電場(chǎng)則通過(guò)差分求解泊松方程得到，對(duì)于二維軸對(duì)稱模型，在柱坐標(biāo)系中，泊松方程可以簡(jiǎn)化為

式中：為電勢(shì)；為電荷密度；為真空介電常量；為徑向距離；為縱向深度。

泊松方程的求解采用五點(diǎn)差分格式，使用超松弛迭代法，通過(guò)改變松弛因子以加速收斂。為了減少仿真計(jì)算量，采用虛擬陰極，位于推力器的中部。選擇準(zhǔn)中性模型作為虛擬陰極模型，在計(jì)算過(guò)程中的每個(gè)時(shí)間步長(zhǎng)內(nèi)，通過(guò)區(qū)域邊界向等離子體發(fā)射一定量的電子，以維持陰極邊界內(nèi)的電子和離子滿足準(zhǔn)中性條件。計(jì)算陰極邊界所有網(wǎng)格的總的凈離子數(shù)Δ，具體的處理方式如下：

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

式中：為碰撞概率；為中性原子的數(shù)密度；為電子和中性原子的相對(duì)速度（可近似為電子速度）；為電子和中性原子的碰撞截面的總和，包括彈性碰撞截面、激發(fā)碰撞截面和電離碰撞截面，與電子能量相關(guān)；Δ為時(shí)間步長(zhǎng)。

內(nèi)部碰撞過(guò)程主要為電子和原子、離子之間的碰撞，忽略中性粒子相互之間的彈性碰撞。離子相比原子有著較高的速度和能量，加速噴出的離子即為推力器推力來(lái)源。粒子的運(yùn)動(dòng)和加速通過(guò)解耦計(jì)算獲得，求解電場(chǎng)后將電場(chǎng)分配至粒子，得到粒子所受的電場(chǎng)力，通過(guò)靜態(tài)磁場(chǎng)分布得到粒子所受的磁場(chǎng)力，進(jìn)一步計(jì)算粒子的加速度，更新粒子的速度和位置。

和氙氣相比，碘蒸汽在放電通道內(nèi)部和陰極發(fā)射的電子相互作用，存在電離、解離、復(fù)合等多種不同的碰撞類型，產(chǎn)生的粒子種類更多，放電過(guò)程更為復(fù)雜。I和I的電離能分別為10.45 eV 和9.40 eV，Xe的電離能為12.10 eV。單個(gè)粒子的I和I的電離碰撞截面約為6.0×10cm和12.3×10cm，大于Xe 的4.8×10cm。I的解離能約為1.57 eV。在仿真中，主要實(shí)現(xiàn)了解離-電離過(guò)程，主要考慮的反應(yīng)如下：

式中：e 為電子。

相關(guān)碘原子的彈性、激發(fā)、電離碰撞截面參數(shù)從文獻(xiàn)［14-16］中獲取。在模型中主要考慮的電荷交換碰撞（CEX 碰撞）反應(yīng)如下：

上述CEX 碰撞截面表達(dá)式在文獻(xiàn)［17］中已經(jīng)給出：

式中：為CEX 碰撞截面；為相對(duì)碰撞速度；為碘的電離能，約10.45 eV；為氫的電離能，約13.60 eV；系 數(shù)為1.81×10；系 數(shù)為2.12×10。文獻(xiàn)［7］中也采用了相同的公式進(jìn)行模擬。

整體仿真流程可簡(jiǎn)述為：初始化計(jì)算模型，對(duì)計(jì)算區(qū)域進(jìn)行網(wǎng)格劃分，將磁場(chǎng)數(shù)據(jù)耦合加載到計(jì)算程序中，初始化物理參數(shù)，載入初始粒子。根據(jù)工質(zhì)流量，在計(jì)算的每個(gè)時(shí)間步長(zhǎng)內(nèi)射入對(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 方法求解原子、離子間的碰撞和動(dòng)量交換、電荷交換；通過(guò)MCC 方法求解電子和原子間的碰撞；通過(guò)求解泊松方程得到電勢(shì)、電場(chǎng)分布，再根據(jù)加載的磁場(chǎng)數(shù)據(jù)計(jì)算粒子的加速度和速度。當(dāng)流場(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ū)域

推力器在周向不同方位角各類參數(shù)基本分布均勻，為了簡(jiǎn)化計(jì)算，仿真模型采用二維軸對(duì)稱型結(jié)構(gòu)。數(shù)值模擬時(shí)所選取的計(jì)算區(qū)域主要包括推力器的放電室、出口和部分羽流區(qū)域，左側(cè)邊界為推力器陽(yáng)極邊界，以推力器中軸線為對(duì)稱軸，相關(guān)邊界條件設(shè)置如圖1 所示。碘原子和內(nèi)外壁面之間的作用采用漫反射模型。根據(jù)推力器的尺寸，放電通道的長(zhǎng)度=20.00 mm，內(nèi)壁面處于=13.38 mm 處，外壁面處于=25.00 mm 處，在加速通道外部選取半徑為38.00 mm、軸向長(zhǎng)度為16.00 mm 的圓柱形羽流場(chǎng)區(qū)域。

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

1.3 加速方法

碘離子和電子的質(zhì)量相差約5 個(gè)數(shù)量級(jí)，在相同的能量條件下，電子速度將是碘離子速度的百倍以上，如果按照電子的時(shí)間步長(zhǎng)同時(shí)推進(jìn)碘離子，需花費(fèi)大量時(shí)間和計(jì)算資源以達(dá)到穩(wěn)定，采取降低碘離子質(zhì)量和提高介電常數(shù)的辦法來(lái)提高計(jì)算速度。具體處理方式如下：將碘離子質(zhì)量降低到1/2 500，則其運(yùn)動(dòng)速度相應(yīng)提高50 倍；將介電常數(shù)提高100 倍，則德拜長(zhǎng)度相應(yīng)提高10 倍，可以降低達(dá)到穩(wěn)定時(shí)所需要的時(shí)間步，減少網(wǎng)格的總數(shù)量。

2 數(shù)值模擬以及分析對(duì)比

2.1 磁場(chǎng)分布

推力器磁場(chǎng)考慮為靜態(tài)磁場(chǎng)，由ANSYS 商用軟件計(jì)算得到分布結(jié)果，以輸入文件形式將磁場(chǎng)數(shù)據(jù)加載至模擬計(jì)算程序。計(jì)算所得放電通道中軸線上磁感應(yīng)強(qiáng)度分布如圖2 所示，從陽(yáng)極到羽流近場(chǎng)，整體呈先上升再下降趨勢(shì)，最大值253.51×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 axis

2.2 設(shè)計(jì)工況下內(nèi)部參數(shù)分布及性能參數(shù)

對(duì)設(shè)計(jì)工況下200 W 級(jí)推力器工作情況進(jìn)行仿真，得到放電通道內(nèi)的粒子數(shù)密度分布、電子溫度和通道內(nèi)離子軸向速度等結(jié)果。設(shè)計(jì)工況的輸入?yún)?shù)包括：功率為200 W、放電電壓為200 V、工質(zhì)流量為13.83 sccm。初始電場(chǎng)通過(guò)設(shè)置左側(cè)陽(yáng)極邊界為200 V、右側(cè)自由邊界為0 V 求解得到；穩(wěn)態(tài)電場(chǎng)通過(guò)使用文獻(xiàn)［19］中的方法，迭代求解方程組得到。離子數(shù)密度分布如圖3 所示，最大值出現(xiàn)在通道中間區(qū)域，最大值約為4.12×10m，隨著通道內(nèi)電磁場(chǎng)對(duì)碘離子加速，數(shù)密度迅速降低，在通道出口位置的碘離子數(shù)密度約為7.66×10m。沿通道軸向離子數(shù)密度先增加后降低，分布規(guī)律和HOFER 等的結(jié)果相似，通道出口處碘離子數(shù)密度和文獻(xiàn)［7］中的結(jié)果相似。推力器放電通道內(nèi)分為近陽(yáng)極區(qū)、電離區(qū)和加速區(qū)，滿足放電通道粒子的分布特征。通道內(nèi)的速度分布如圖4 所示，離子速度在通道后半段迅速增加，并在電場(chǎng)作用下保持加速效果至通道出口，最大速度約為16 860 m/s。

圖3 設(shè)計(jì)工況下離子數(shù)密度分布Fig.3 Distribution of ion number density under design

圖4 設(shè)計(jì)工況下離子軸向速度分布Fig.4 Distribution of the ion axial velocity under the design condition

通道內(nèi)碘原子數(shù)的密度先增大后減少分布情況如圖5 所示，在內(nèi)外壁面處積附了大量的碘原子，碘原子的數(shù)密度達(dá)到最高值；在電離區(qū)則由于大量碘原子在此處發(fā)生電離，碘原子的數(shù)密度較低。碘原子軸向的速度分布情況如圖6 所示，離子被電場(chǎng)加速，在加速區(qū)速度較高，和中性原子發(fā)生碰撞后，區(qū)域內(nèi)中性原子的速度也有所增加，達(dá)到千米每秒量級(jí)；其他區(qū)域內(nèi)碘原子的速度量級(jí)則和文獻(xiàn)［12］中一致。

圖5 碘原子數(shù)密度分布Fig.5 Distribution of the number density of iodine atoms

圖6 碘原子軸向速度分布Fig.6 Distribution of the axial velocity of iodine atoms

2.3 碘工質(zhì)與氙工質(zhì)對(duì)比

使用氙工質(zhì)代替碘工質(zhì)，在設(shè)計(jì)工況下進(jìn)行數(shù)值模擬，比較2 種工質(zhì)的性能。通道內(nèi)氙離子的數(shù)密度分布如圖7 所示，其分布趨勢(shì)和碘離子大致相同，最大值出現(xiàn)在通道中央，約為1.38×10m，出口處最大軸向速度達(dá)到了約15 150 m/s。相較于碘工質(zhì)，在通道內(nèi)沿通道中央的氙原子數(shù)密度，從陽(yáng)極開始就一直呈降低趨勢(shì)，由于近陽(yáng)極區(qū)電離率較低，下降趨勢(shì)較平緩，在電離區(qū)則快速降低。對(duì)于碘工質(zhì)，其沿通道中央的碘原子數(shù)密度呈先增大后降低的分布趨勢(shì)，在陽(yáng)極附近碘原子數(shù)密度較低，在軸向1～3 mm 處碘原子數(shù)密度明顯增加，再在電離區(qū)快速降低，如圖7 所示。說(shuō)明碘作為雙原子分子，其機(jī)理和氙原子存在差異。由于存在解離過(guò)程，并且解離閾值較低，在中性碘分子通過(guò)氣體分配器進(jìn)入到放電室后，需先和電子碰撞解離為碘原子，之后再進(jìn)行下一步反應(yīng)。碘分子數(shù)密度分布也證明了這一點(diǎn)，如圖8 所示。在陽(yáng)極附近處碘分子數(shù)密度較高，碘原子數(shù)密度較低，因?yàn)榈夥肿觿傔M(jìn)入通道，還未充分解離；然后碘分子數(shù)密度快速降低，同處的碘原子數(shù)密度則明顯上升，說(shuō)明解離過(guò)程在此處頻繁發(fā)生。對(duì)于碘工質(zhì)霍爾推力器而言，除了傳統(tǒng)的近陽(yáng)極區(qū)、電離區(qū)和加速區(qū)外，還存在解離區(qū)，在此區(qū)域內(nèi)，大量的碘分子解離成碘原子，其位置在近陽(yáng)極區(qū)之后、電離區(qū)之前。在此計(jì)算條件下，其寬度略大于2 mm。為了提高推力器的性能，應(yīng)設(shè)法提前解離區(qū)的位置，使碘原子充分電離，提高推力器效率。

圖7 氙工質(zhì)離子數(shù)密度分布Fig.7 Distribution of number density of Xe+

圖8 碘分子數(shù)密度分布Fig.8 Distribution of the number density of I2

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

本文基于PIC/DSMC/MCC 混合方法，數(shù)值模擬了設(shè)計(jì)工況下碘工質(zhì)霍爾推力器放電通道內(nèi)部的等離子體參數(shù)，并和氙工質(zhì)進(jìn)行對(duì)比，為未來(lái)200 W級(jí)碘工質(zhì)霍爾推力器的設(shè)計(jì)優(yōu)化和地面試驗(yàn)提供了理論支撐，計(jì)算結(jié)果表明：

1）設(shè)計(jì)工況下，通道內(nèi)的碘離子數(shù)密度最高達(dá)到4.12×10m，隨著通道內(nèi)自洽電場(chǎng)對(duì)離子的加速作用，在出口位置下降至7.66×10m，最大出口速度約為16 860 m/s，通道內(nèi)分為明顯的近陽(yáng)極區(qū)、電離區(qū)和加速區(qū)。

2）碘原子在通道內(nèi)沿通道中央的數(shù)密度先增大后減小，在內(nèi)外壁面處積附了大量的碘原子，在電離區(qū)碘原子的數(shù)密度迅速降低，在加速區(qū)碘原子軸向速度有所上升。

3）相較于氙工質(zhì)，碘工質(zhì)還存在解離區(qū)，位于近陽(yáng)極區(qū)之后、電離區(qū)之前，在解離區(qū)內(nèi)，大量碘分子解離為碘原子。在此計(jì)算條件下，解離區(qū)寬度略大于2 mm。