2 cm 電子回旋共振離子推力器離子源中磁場對等離子體特性與壁面電流影響的數值模擬*

夏旭 楊涓? 付瑜亮 吳先明 耿海 胡展

1) (西北工業大學 航天學院 西安 710072)

2) (蘭州空間技術物理研究所 蘭州 730000)

電子回旋共振離子推力器(electron cyclotron resonance ion thruster, ECRIT)離子源內等離子體分布會影響束流引出, 而磁場結構決定的ECR 區與天線的相對位置共同影響了等離子體分布.在鞘層作用下, 等離子體中的離子或電子被加速對壁面產生濺射, 形成壁面離子或電子電流, 造成壁面磨損和等離子體損失, 因此研究壁面電流與等離子體特征十分重要.為此本文建立2 cm ECRIT 的粒子PIC/MCC (particle-in-cell with Monte Carlo collision)仿真模型, 數值模擬研究磁場結構對離子源內等離子體與壁面電流特性的影響.計算表明, 當ECR 區位于天線上游時, 等離子體集中在天線上游和內外磁環間, 柵極前離子密度最低, 故離子源引出束流、磁環端面電流和天線壁面電流較低.ECR 區位于天線下游時, 天線和柵極上游附近的等離子體密度較高, 故離子源引出束流、天線壁面電流和磁環端面電流較高.腔體壁面等離子體分布與電流受磁場影響最小.

1 引 言

ECRIT 具有無陰極燒蝕特點, 除了用于空間飛行器控制[1,2], 還能用于材料處理等方面[3,4].目前國外研究出10 cm 和2 cm ECRIT, 并分別應用于HAYABUSA[2,5]和50 與70 kg 微小衛星[6,7].目前國內外熱衷于微小衛星推進裝置的研究, 因此開展2 cm ECRIT 研究具有重要意義.離子源是2 cm ECRIT 重要的部件, 其內部物理過程包括電離、碰撞以及鞘層加速帶電粒子對壁面轟擊, 這些過程決定了等離子體和壁面電流特征.等離子體特征影響推力器引出束流大小, 同時推力器內壁面濺射電流是約束工作壽命和引起等離子體損失的重要因素,而且等離子體和壁面電流存在密切關系并同時受磁場結構和工作參數等的影響, 因此開展不同磁場結構對等離子體和壁面電流特征的影響研究十分重要.由于實驗研究[8?10]很難觀測到離子源內瞬態物理過程變化, 需要借助數值模擬手段研究離子源內部等離子體過程.雖然PIC 模擬已經揭示了ECRIT 的等離子體特征[11?14], 但2 cm ECRIT 離子源壁面電流數值模擬少有報道.為此本文采用PIC/MCC 方法研究磁場結構對離子源內等離子體與壁面電流特性的影響, 為推力器的性能預估提供參考.

2 物理和計算模型

2.1 物理模型

如圖1 所示, 2 cm ECRIT 離子源由圓柱腔體、環形天線、磁軛、內外環形永磁體和雙柵極組成, 其中兩磁環分布在離子源底部, 其間環向均布8 個進氣孔; 4.2 GHz 的微波能量通過環形天線饋入放電室, 在0.15 T 的磁場區域形成ECR 區[15].在ECR 區, 電子繞磁力線的回旋運動和微波電場變化同步產生諧振, 由此微波電場加熱電子形成高能電子, 高能電子和原子發生碰撞激發和電離, 從而形成ECR 等離子體, 等離子體中的離子經由雙柵極被高速引出.圖1 中的H1,W1,H2和W2分別表示外磁環高度、外磁環寬度、內磁環高度和內磁環寬度.L1為天線環形段下表面與屏柵上表面之間的距離,L2為內磁環下表面與天線環形段上表面之間的距離, 本文中不同磁場結構中L1和L2的值相同.定義A, B, C 和D 分別表示天線環形段的內表面, 下表面, 外表面和上表面; E, F 分別表示內磁環端面、外磁環端面; G 表示腔體內表面.

圖1 2 cm ECRIT 離子源結構示意圖Fig.1.Schematic diagram of the 2 cm ECRIT ion source internal structure.

2.2 計算模型

圖2 2 cm ECRIT 離子源的放電形貌Fig.2.Discharge image of 2 cm ECRIT ion source.

圖2 為2 cm ECRIT 離子源的氣體放電圖像.從圖2 中可以看出,z軸周圍有一個圓形的等離子體區域, 其亮度較均勻, 環形天線的陰影區域是可見的.由圖可見, 離子源內的等離子體分布近似軸對稱分布.文獻[9]也表明2 cm ECRIT 離子源中靜磁場軸對稱, 高頻場接近于軸對稱性.基于以上分析, 計算域等離子體參數分布具有軸對稱性, 因此將放電室內的三維的環形電離區簡化為二維軸對稱區域.

2.3 全粒子PIC/MCC 模型

2.3.1 電磁場

本文將ECRIT 離子源內部的電磁場分解為3 個模塊: 靜電場、靜磁場和微波高頻電場, 這里忽略高頻磁場, 因為其比靜磁場小得多.靜電場通過離散泊松方程?·?φ=?ρ/ε0直接求解得到, 其中φ,ρ和ε0分別為電勢, 電荷密度和真空介電常數.除了對稱邊界條件為其余邊界電勢均為0 V; 靜磁場作為背景場, 由COMSOL 多物理場軟件計算, 靜磁場的計算參數設置: 永磁體(釤鈷合金)的相對磁導率為1.36; 磁軛(軟鐵)的相對磁導率為4000; 空氣的相對磁導率為1.靜磁場的求解域為物理模型的所有區域, 邊界條件僅需設置外邊界為零磁標量勢邊界.微波高頻電場采用COMSOL 計算, 計算時不考慮等離子體電流, 這是由于微波功率較低時等離子體電流對微波電場影響甚微[16].微波電場的計算參數設置: 壁面(不銹鋼)的相對介電常數為1, 電導率為 4×106S/m ;天線(銅)的相對介電常數為1, 電導率為 6×107S/m ;空氣的相對介電常數為1, 電導率為0 S/m.微波電場的求解域為放電室區域, 邊界條件僅需設置微波輸入端口邊界, 輸入端口類型設置為同軸, 輸入功率設為2 W.PIC 程序計算時, 每個時間步上高頻場的瞬時值與靜電場疊加共同推進電子運動; 而離子質量較大, 故高頻場與離子作用弱, 認為離子不受高頻場的影響, 只有靜電場推動離子運動.

2.3.2 粒子碰撞

離子源以氙氣為工質, 假設中性原子為背景粒子, 均勻分布在計算域中.考慮的碰撞類型有: 電子與原子的彈性碰撞、激發碰撞和一級電離碰撞,以及帶電粒子間的庫倫碰撞.電子與中性原子之間的碰撞是通過蒙特卡洛(MCC)方法處理, 碰撞截面采用Szabo[17]的截面數據; 采用直接蒙特卡洛(DSMC)方法來描述帶電粒子間的庫倫碰撞[18].帶電粒子的運動采用Boris“蛙跳”算法求解[17].粒子運動的邊界條件: 對稱邊界時, 粒子在該邊界上反射, 其余邊界上, 粒子完全吸收.

圖3 不同磁場結構下磁場分布: (a) 1 號源; (b) 2 號源; (c) 3 號源Fig.3.Distribution of magnetic flux density inside of the discharge chamber: (a) 1; (b) 2; (c) 3.

2.3.3 磁場結構

如表1 所示, 計算中模擬了3 種不同磁場結構.通過COMSOL 軟件計算離子源中磁場強度,其ECR 區如圖3 所示.可以看出, 1, 2 和3 號源的ECR 區分別位于天線上游、與天線相交、天線下游.隨著外磁環的高度和寬度增大, ECR 區逐漸向屏柵移動并遠離天線.

表1 磁場結構參數Table 1.Magnetic circuit structure parameters.

2.4 計算條件

計算條件為氙氣流量0.3 sccm (1 sccm =98.3 μg/s), 微波功率2 W, 微波頻率4.2 GHz.中性粒子密度當作背景氣體代入PIC 程序, 0.3 sccm流量時, 估算得到的中性原子密度范圍為1.5×1020m?3[19].

根據文獻[12]給出的離子源內等離子體密度為 2×1017m?3, 再假設平均電子溫度為10 eV, 則德拜長度約為0.053 mm.保持真實的介電常數下,為避免數值過熱, 將網格步長設定為0.025 mm,即使等離子體密度達到5 × 1017m–3也能保證計算精度.

3 結果和討論

3.1 計算結果的收斂性

計算開始時, 設初始電子和離子溫度分別為2,0.5 eV, 初始等離子體密度為 1×1016m?3.圖4 記錄了不同磁場結構離子源的帶電粒子數隨計算時間的變化.如圖4 所示, 隨著計算時間的進行,在50 萬步以后, 電子和離子數都逐漸達到了線性穩態.

圖4 帶電粒子數隨計算步的變化Fig.4.Quantities of charged particles versus computation step.

圖5 不同磁場結構下等離子體分布: (a)電子密度; (b)離子密度Fig.5.Plasma distribution of ion source with different magnetic circuit: (a) Electron density; (b) ion density.

3.2 等離子體特性分布

不同磁場結構下離子源內部的電子和離子密度分布計算結果如圖5 所示.從圖中可以看出不同離子源內的電子和離子密度峰值都出現在ECR 區附近, 且電子密度和離子密度峰值分別是3.5×1017m?3/3.9×1017m?3,2.6×1017m?3/3.0×1017m?3, 4.1×1017m?3/4.2×1017m?3.由此可見, 3 種離子源內部的等離子體密度峰值均大于4.2 GHz微波頻率的截止密度ncut-off=ε0meω2/e2≈2.1×1017m?3[20,21], 其中ω為微波頻率, 說明了ECR等離子體的高密度特征.不同離子源內的等離子體密度峰值區面積大小為: 3 號源 > 1 號源 > 2 號源.2 號源的密度峰值區域最小, 是因為ECR 區被天線切斷, 從而影響磁鏡中的電子被ECR 區持續加熱.1 號源的等離子體密度峰值位于天線上游,并且ECR 區位于較小的磁鏡區, 所以天線上游的內外磁環間的等離子體密度較高而靠近屏柵區域的等離子體密度較低.這是因為天線環形段阻礙了等離子體擴散.3 號源的ECR 區面積大于1, 2 號源, 而且ECR 區位于較大的磁鏡區, 所以其等離子體密度最高而且靠近柵極前的等離子體高于1,2 號源.柵極前的等離子體密度越高越有利于離子束流的引出, 符合文獻[15]的實驗結果, 即ECR區位于天線下游的離子源引出效果優于ECR 區位于天線左側的離子源.

圖6 為不同磁場結構下離子源的電子溫度分布.1, 2 和3 號源中電子在ECR 區被加熱后, 沿磁力線擴散到ECR 區附近的磁鏡中, 電子溫度峰值分別為13.7 eV, 13.3 eV 和12.7 eV.1 號源中電子溫度最高, 因為電子ECR 加熱區靠近微波電場幅值高的天線上游[15].與圖6 對比, 1 號源的電子溫度峰值出現在內外磁環間的ECR 區和天線上表面附近的ECR 區附近, 但離子密度峰值出現在天線上游附近的ECR 區.這是因為天線上表面ECR 區的高能電子未與原子碰撞電離就損失在天線上表面, 而內外磁環間ECR 區的高能電子會與原子發生電離碰撞.2 號源的電子溫度峰值靠近在天線下表面ECR 區, 而高密度區偏移了ECR 區.這也是因為未發生電離碰撞的高能電子損失在天線下表面.3 號源中電子被ECR 區加熱后, 沿磁力線擴散到ECR 區和柵極上游附近的磁鏡中.柵極上游磁鏡中的高能電子多于ECR 區附近的高能電子, 這是因為柵極上游的磁鏡范圍更大, 而且ECR區附近的高能電子可能損失在天線環形段表面或與原子發生電離碰撞.3 號源的高密度區位于ECR區和柵極上游之間的磁鏡中.這是ECR 區加熱產生的高能電子向柵極擴散中, 部分高能電子與原子發生電離碰撞, 從而ECR 區和柵極上游之間的電子溫度較低.部分高能電子和未與原子發生電離碰撞的電子可能在磁鏡的約束下又往復穿越ECR區, 被ECR 區加熱并約束在柵極上游的磁鏡區.2 和3 號源在靠近柵極表面附近有少量高能電子(大于10 eV), 這是因為壁面附近鞘層的歐姆加熱效應.

圖6 不同磁場結構下離子源的電子溫度Fig.6.Electron temperature of ion source with different magnetic circuit.

3.3 壁面電流

圖7 不同磁場結構下離子源內各壁面的電流密度統計Fig.7.Current density statistics on each surface of ion source with different magnetic circuit.

4 結 論

本文采用全粒子PIC/MCC 方法模擬了2 cm ECRIT 離子源內部的等離子體與壁面電流.結論如下:

(1)電子被磁鏡約束在天線ECR 區附近, 使得等離子體也分布在ECR 區附近.并且ECR 區與天線的相對位置會共同影響離子源的引出束流.當ECR 區位于天線下游時, 更有利于離子引出.從等離子體分布結果中發現柵極前等離子體主要集中在外磁環與軸線之間, 推測當保持孔數不變時, 而減少壁面附近的柵極孔分布范圍, 將柵極孔集中分布在軸線中心附近, 能提高引出束流以及推進劑利用效率.

(2) ECR 區位于天線下游時, 等離子體主要分布在天線下游而且柵極前等離子體密度高, 但離子對天線濺射嚴重.這樣的磁場設計雖然能夠提高引出的離子束流, 但會使天線工作壽命縮減.不同壁面的電流分布規律為: 天線壁面 > 磁體端面 > 腔體壁面.由此可見離子源天線是壁面電流的主要來源, 有可能是影響工作壽命的主要因素.未來設計推力器時, 應權衡引出性能和壽命來選取合適的磁場結構.