碘工質霍爾推力器內部過程數值模擬

王 尚，梁子軒，王平陽，徐宗琦，杭觀榮

（1.上海交通大學 機械與動力工程學院，上海 200240；2.上海空間推進研究所上海空間發動機工程技術研究中心，上海 201112）

0 引言

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

對于碘工質推力器放電過程的數值模擬，國內外學者均開展了相關研究。法國的GRONDEIN等采用流體動力學模型對碘工質離子推力器進行了數值模擬，獲取了離子、電子數密度和電子溫度等特征參數，將碘和氙工質情況下的推力、比沖和效率進行了比較；NASA 的MARIA采用混合單元粒子法模型對200 W 霍爾推力器進行了模擬，對比了碘和氙工質情況下羽流區的中性氣體和離子的數密度分布；法國的LUCKEN 等針對離子推力器進行數值模擬，研究了沿推力器軸向方向的磁場強度變化時，碘工質和氙工質情況下推力器的推力和效率變化；哈爾濱工業大學的牛翔等采用一維流體動力學模型對碘工質會切場推力器進行數值模擬，研究了碘工質會切場推力器內部的等離子體參數分布。相較于上述研究成果，目前對百瓦級碘工質霍爾推力器放電通道內部過程的數值模擬研究工作還鮮有報道。

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

1 數值仿真

1.1 計算基本模型

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

基于PIC 計算時采用了靜電求解模型，電場則通過差分求解泊松方程得到，對于二維軸對稱模型，在柱坐標系中，泊松方程可以簡化為

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

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

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

式中：為碰撞概率；為中性原子的數密度；為電子和中性原子的相對速度（可近似為電子速度）；為電子和中性原子的碰撞截面的總和，包括彈性碰撞截面、激發碰撞截面和電離碰撞截面，與電子能量相關；Δ為時間步長。

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

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

式中：e 為電子。

相關碘原子的彈性、激發、電離碰撞截面參數從文獻［14-16］中獲取。在模型中主要考慮的電荷交換碰撞（CEX 碰撞）反應如下：

上述CEX 碰撞截面表達式在文獻［17］中已經給出：

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

整體仿真流程可簡述為：初始化計算模型，對計算區域進行網格劃分，將磁場數據耦合加載到計算程序中，初始化物理參數，載入初始粒子。根據工質流量，在計算的每個時間步長內射入對應數目的中性碘分子，根據準中性條件射入對應數目的電子，依據粒子的現有速度計算其運動狀態。通過PIC 方法建立粒子間的電磁場作用模型；通過DSMC 方法求解原子、離子間的碰撞和動量交換、電荷交換；通過MCC 方法求解電子和原子間的碰撞；通過求解泊松方程得到電勢、電場分布，再根據加載的磁場數據計算粒子的加速度和速度。當流場中的粒子凈流量小于設置值時，視為計算進行至推力器完全穩定工作狀態，輸出計算結果；如不穩定則繼續進行計算，待計算至穩定狀態后，再統計相關參數。

1.2 計算區域

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

圖1 推力器結構及流場仿真區域Fig.1 Construction and simulated zone of the thruster

1.3 加速方法

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

2 數值模擬以及分析對比

2.1 磁場分布

推力器磁場考慮為靜態磁場，由ANSYS 商用軟件計算得到分布結果，以輸入文件形式將磁場數據加載至模擬計算程序。計算所得放電通道中軸線上磁感應強度分布如圖2 所示，從陽極到羽流近場，整體呈先上升再下降趨勢，最大值253.51×10T出現在通道出口附近，符合霍爾推力器內部磁場設計準則和分布規律。

圖2 通道軸線磁感應強度分布Fig.2 Distribution of the magnetic flux density along the axis

2.2 設計工況下內部參數分布及性能參數

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

圖3 設計工況下離子數密度分布Fig.3 Distribution of ion number density under design

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

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

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

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

2.3 碘工質與氙工質對比

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

圖7 氙工質離子數密度分布Fig.7 Distribution of number density of Xe+

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

3 結束語

本文基于PIC/DSMC/MCC 混合方法，數值模擬了設計工況下碘工質霍爾推力器放電通道內部的等離子體參數，并和氙工質進行對比，為未來200 W級碘工質霍爾推力器的設計優化和地面試驗提供了理論支撐，計算結果表明：

1）設計工況下，通道內的碘離子數密度最高達到4.12×10m，隨著通道內自洽電場對離子的加速作用，在出口位置下降至7.66×10m，最大出口速度約為16 860 m/s，通道內分為明顯的近陽極區、電離區和加速區。

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

3）相較于氙工質，碘工質還存在解離區，位于近陽極區之后、電離區之前，在解離區內，大量碘分子解離為碘原子。在此計算條件下，解離區寬度略大于2 mm。