考慮電荷交換的柵極區離子流數值模擬

劉 洋，李 娟，楚豫川，曹 勇

(1.哈爾濱工業大學深圳研究生院，廣東深圳518055;2.蘭州空間技術物理研究所，真空低溫技術與物理重點實驗室，甘肅蘭州730000)

1 引言

離子推力器是利用電能電離、加速推進劑離子并高速噴出后，獲得反推力的一種電推進裝置。其比沖是常規化學火箭的數倍到數十倍［1］。與常規化學火箭相比，電推進系統具有比沖高、壽命長等顯著特點，可有效減小推進劑的總質量，增加有效載荷，適應了空間飛行器的發展對推進系統提出的要求。因此，離子推力器系統越來越多地受到人們的關注，目前已成為我國航天領域的迫切需求［2］。

受電場等因素的限制，離子推力器的推力通常較小，因而，要完成一定的任務，離子推力器需要更長的工作時間。這就需要離子推力器有更長的工作壽命，一般要達到上萬小時。離子推力器光學系統是推力器的關鍵部件之一，它決定離子推力器的性能和壽命。離子光學系統由柵極及其支撐結構組成，本文中離子光學系統柵極包括屏柵極和加速柵極。其失效形式主要有3種:屏柵極和加速柵極短路，柵極結構失效和柵極系統下游等離子體中的電子回流［3］。柵極的腐蝕是影響柵極壽命的主要因素。其中加速柵極下游區域電荷交換碰撞形成的電荷交換離子對加速柵極的濺射腐蝕最嚴重［4］。柵極長時間工作后，加速柵極可能產生穿孔、加速柵極孔變大等現象，嚴重影響離子推力器的工作性能和使用壽命。因此研究光學系統的腐蝕情況對于離子推力器的壽命評估有重大意義。

目前研究光學系統腐蝕主要采用實驗和數值模擬2種方法。實驗方法可以直接準確的獲得真實的實驗數據，但需投入大量的資金和時間，消耗巨大。與實驗方法相比，數值方法結果直觀，實驗資金投入較小且大大降低了時間消耗，因此數值方法已被廣泛采用。作者基于嵌入式有限元－質點網格［5～7］與蒙特卡洛［8，9］(IFE－PIC－MCC)的三維等離子體離子模擬程序，對離子推力器光學系統中離子密度分布及加速柵極腐蝕情況進行了研究。

2 離子運動計算模型設置

2.1 計算域和邊界

離子推力器的柵極是很薄的電極板，柵極上開有數千個柵極孔，柵極孔呈六邊形分布，具有對稱結構。在計算模擬中，根據柵極的對稱性，從柵極區域中選取包含柵極所有信息的最小區域進行計算。圖1按照灰色區域將光學系統劃分為長方形區域。計算域中包含了2個1/4的柵極孔，這樣選取的計算域包含六邊形分布的柵極孔的所有信息，也就是將相同的計算域拼接在一起后會形成完整的柵極排列。這樣便可以認為所選取的計算域包含了整個柵極孔的信息。計算模型的參數都經過無量綱化。柵極幾何參數等如表1所列。

圖1 計算域二維簡圖

對于中性原子，只需考慮粒子運動邊界，無需考慮電場條件。當中性原子越過上游邊界或下游邊界時，將該原子刪除，如果中性原子到達4個側面中的任意一個或柵極，則將該原子鏡面反射回計算區域。

表1 柵極幾何參數

2.2 CEX 的產生與加速柵極腐蝕

由于加速柵極下游表面的腐蝕主要是由交換電荷離子造成的，在離子束引出的過程中，如果不考慮離子與原子之間的碰撞，就忽視了交換電荷離子，只能計算得到部分束流對加速柵極孔壁的腐蝕情況，而無法獲得加速柵極下游表面的腐蝕情況［10］。因此，要計算加速柵極下游表面的腐蝕，必須考慮離子之間的碰撞。離子之間的碰撞，主要包括離子與中性原子之間的電荷交換碰撞(CEX)，離子之間的庫侖碰撞，中性原子之間的動量交換碰撞等。其中，中性原子與離子之間的電荷交換碰撞是加速柵極電流產生及柵極腐蝕的主要原因［11］。因此，在計算加速柵極腐蝕過程中，本文只考慮中性原子與離子之間的CEX碰撞，CEX碰撞是經過電場加速的高速離子和熱運動的中性原子碰撞，產生慢速的電荷交換離子和快速的中性原子。即

作者用IFE－PIC方法模擬離子在離子光學系統中的運動，加入了蒙特卡洛［12］碰撞模塊，來模擬考慮CEX離子在離子光學系統中的產生和運動。空間任意一個帶電離子與空間分布的中性原子在每個步長Δt內的碰撞頻率［12］為

碰撞概率為

式中 nn(xi)為離子位置的中性原子的數密度;vi為帶電離子的速度;σT(vi)為離子電荷交換碰撞截面，采用經驗公式［13］

式中 k1=－1.303×10－10;k2=30×10－10;v為帶電離子與中性原子的相對速度，由于帶電離子的速度遠大于只有以熱運動為主的中性原子的速度，所以本文用帶電離子的速度來代替v。

在模擬計算中，離子在一個時間步長內運動的距離不超過最小網格長度的1/2，即

式中 vmax為計算域內模擬離子的最大速度;Δz為z方向的網格長度。

帶電離子與中性原子的碰撞概率很小，但需要計算所有離子的碰撞概率與碰撞截面，非常浪費時間，所以本文采用零碰撞技術［13］來簡化蒙特卡洛法的計算過程。首先，找出所有離子的最大碰撞頻率νm

則可得到離子的最大碰撞概率為

假設在計算域內總的帶電離子數為N，隨機選擇NPnull個離子，分別計算每個離子的碰撞頻率ν。選一個均勻分布的隨機數RAN，如果就認為該離子發生了電荷交換碰撞。離子發生電荷交換碰撞后，將該位置的中性原子的速度賦給帶電離子，產生的交換電荷離子將參與下一時間步的電場計算和離子運動過程。

離子撞擊到材料的表面，使材料脫離本體的效應稱為濺射效應，它是離子發動機加速柵極下游表面電荷交換碰撞腐蝕形成的物理原因。每個入射離子所轟擊出的材料原子的個數稱為濺射產額，用Y表示。濺射產額計算公式［14］如下

單個模擬離子的入射能量EXe計算公式為

入射離子的入射角θi對濺射產額的影響是通過上述公式得到的濺射產額進行修正而獲得:

3 計算流程

IFE－PIC－MCC三維等離子體離子模擬程序計算離子(包括CEX離子)和中性原子在離子光學系統中的運動。程序可分為3個模塊。首先在不考慮電荷交換碰撞的情況下，分別計算離子和中性原子在計算空間的分布。再利用離子和中性原子的空間分布，根據MCC碰撞模塊計算CEX離子的產生，最后得到包含CEX離子的所有離子在計算空間的分布。

(1)模擬不考慮電荷交換碰撞的離子在離子光學系統中的運動

在每個時間步長內都有一定數目的模擬粒子從計算域的左邊界(z=0平面)進入計算域，在此過程中，程序需要給定的參數包括，每個時間步長進入計算域的模擬粒子的數量，初始位置，初始速度等。進入計算域后，通過給定的邊界條件和粒子所帶的電量，用體積權重法將離子所帶電量向所在網格的節點進行插值，得到網格節點上的電荷密度;利用網格節點上的電荷密度計算出網格節點上的電勢和電場;由網格節點上的電場向單個粒子的位置進行插值，計算出該離子位置的電場強度進而計算出離子所受的力;用牛頓第二定律計算出粒子在下一時刻的位置。在下一個步長中，重復以上步驟，直到計算域中離子的數目達到穩定狀態。

(2)模擬中性原子在離子光學系統中的運動

中性原子不帶電，因此中性原子不受電場的影響，在離子光學系統中勻速運動。在一個時間步長內，通過給定的邊界條件，用體積權重法將粒子向所在網格的節點進行插值，計算出網格節點上的粒子密度。計算下一個步長中原子的位置，重復以上步驟，直到計算域中原子的數目達到穩定狀態。最后統計出中性原子的粒子密度。

(3)CEX在離子光學系統中的產生

在IFE－PIC中加入蒙特卡洛碰撞模塊模擬CEX離子的產生和運動。束流中部分經過電場加速的高速離子和以熱速度運動的中性原子之間會發生電荷交換碰撞，產生CEX離子。生成的CEX離子將參與下一個時間步長的離子運動。

4 模擬結果與分析

柵極孔中心線沿z軸方向的電勢如圖2，從圖中可以看出，屏柵極和加速柵極之間存在很高的電勢差，可以使帶電離子以很高的速度噴出。加速柵極的電勢要低于下游邊界的電勢。圖3為柵極孔中心線在z軸方向速度變化。離子在計算域內的初速度是以波爾茲曼分布給出的。經過屏柵極與加速柵極之間的強電場加速后，離子的速度急劇增加，最終以很高的速度噴出，最高速度能到幾萬米/秒。結合圖2的電勢分布可以看出，在屏柵極和加速柵極之間離子電勢急速下降，離子急劇加速，加速柵極下游區，電勢是上升的，離子在該區域稍有減速。

中性原子的粒子密度見圖4。中性原子不帶電，粒子在計算域中的運動不受電場的影響，做勻速運動。

圖2 柵極孔中心線在z軸方向的電勢變化(電勢以5 V無量綱)

圖3 柵極孔中心線上離子速度Z方向分量的變化(速度以1 917.4 m/s無量綱)

圖4 中性原子的粒子數密度(以1×1017m－3無量綱)

圖5 電荷交換離子的離子數密度(以1×1017m－3無量綱)

圖6 2種情況的離子密度比較

圖5 是初始離子數密度(即放電室等離子體密度)n0=1.5×1017m－3時的電荷交換離子的數密度，電荷交換離子的數密度比離子數密度小得多。圖6是2種情況對應工況下的離子數密度在計算域內的空間分布。上半部分是不考慮電荷交換離子時，帶電離子數密度的空間分布;下半部分為考慮電荷交換的離子數密度的空間分布。從圖中可以看出，由于柵極系統的聚焦效果，從屏柵極開始到加速柵極，主束流的半徑逐漸減小。在正常情況下，主束流半徑在加速柵極孔附近小于加速柵極孔半徑。由于柵極系統這種聚焦作用，使得主束流邊緣的離子數密度比軸線附近的離子數密度小得多。

在加速柵極上游區，2種情況下的離子束密度分布基本沒有差別。在加速柵極下游區，不考慮電荷交換情況下，離子在電場作用下聚焦通過加速柵極孔，沒有撞擊加速柵極?？紤]電荷交換情況下，有離子撞擊加速柵極。因此與加速柵極碰撞的不是主束流離子而是電荷交換離子;離子與加速柵極碰撞產生的加速柵極電流，主要是由電荷交換離子引起的。2種情況下的柵極電流如表2所列?？紤]電荷交換后加速柵極上的柵極電流由0增大到10.68。

表2 柵極電流對比(以8.489×10－8A無量綱)

圖7為加速柵極下游表面電流密度分布數值計算結果。電流密度參數體現的是與加速柵極碰撞的離子在柵極表面的分布情況。在孔中心連線形成的等邊三角形中心附近的電流密度最大。圖8為離子入射角度分布，單位為弧度。圖9為加速柵極下游表面濺射產額分布的數值計算結果。圖10為加速柵極下游表面的腐蝕深度計算結果。腐蝕深度的分布與電流密度的分布相似，與濺射產額分布趨勢相同。圖10中的腐蝕深度為發動機工作時間為1 000 h時的腐蝕深度。

圖7 加速柵極表面電流分布(A/m2)

圖8 離子的入射角分布(rad)

圖9 加速柵極表面濺射產額分布

圖10 加速柵極下游表面腐蝕深度(m)

圖11 為加速柵極孔壁電流密度分布數值計算結果。電流密度參數體現的是離子對柵極孔壁的沖擊情況。圖中顯示靠近加速柵極下游處的離子數密度較大，撞到孔壁的離子數也會較多，在靠近下游的孔壁上產生的電流密度就會增大。圖12為離子入射角度分布。離子的入射角度不同，會影響離子的濺射產額。圖13為加速柵極孔壁濺射產額分布的數值計算結果，濺射產額與入射能量的分布相近，圖中顯示濺射產額在孔壁下游部分會比較大。圖14為加速柵極孔壁的腐蝕深度計算結果。腐蝕深度的分布與電流密度的分布相似，與濺射產額分布趨勢相同。圖14中的腐蝕深度為離子發動機工作時間為1 000 h時的腐蝕深度，最大腐蝕深度靠近加速柵極下游，在孔壁可能會形成凹槽。

圖11 加速柵極孔壁電流分布(A/m2)

圖12 加速柵極孔壁離子的入射角分布(rad)

圖13 加速柵極孔壁濺射產額分布

圖14 加速柵極孔壁腐蝕深度(m)

5 結論

采用IFE－PIC－MCC三維等離子體離子模擬程序模擬了考慮電荷交換的帶電離子在離子光學系統中的運動。電荷交換離子可以使加速柵極沖擊電流增加，電荷交換離子是產生柵極電流和柵極腐蝕的主要因素。由于電荷交換離子對加速柵極的沖擊作用，加速柵極下游表面會形成“坑”形和“凹槽”形腐蝕，孔壁的腐蝕會造成孔徑增大。

［1］SUTTON R G，BIBLARZ O.Rocket Propulsion Elements［M］.John Wiley ＆ Sons，2001.

［2］吳漢基，蔣遠大，張志遠.電推進技術的應用與發展趨勢［J］.推進技術，2003，24(5):385～392.

［3］KAFAFY R.Immersed Finite Element Particle－In－Cell Simulations of Ion Propulsion［D］.Virginia Polytechnic Institute and State University，2005.

［4］PENG X，RUYTEN W，KEEFER D.Three dimensional particle simulation of grid erosion in ion thrusters［R］.IEPC 91－119，1991.

［5］ARAKAWA Y，NAKANO M.An efficient three dimensional optics code for ion thruster research［R］.AIAA 96－3198，1996.

［6］BOND R，LATHAM P.Ion thruster action extraction grid design and erosion modeling using computer simulation［R］.AIAA 95－2923，1995.

［7］HAYAKAWA Y.Three－dimensional numerical model of ion optics system［J］.Journal of Propulsion and Power，1992，8(1):110～117.

［8］GROTE D.Three dimensional simulations of space charge dominated heavy ion beams with applications to inertial fusion energy［D］.University of California，Davis，1994.

［9］SHIRAISH T，KUNINAKA H，SATORI S，et al.Numerical simulation of grid erosion for ion thruster［R］.IEPC 95－90，1995.

［10］WANG J，CAO Y，Kafafy R.Numerical and experimental investigations of crossover ion impingement for subscale ion optics［J］.Journal of Propulsion and Power，1994，24(3):562 ～570.

［11］鐘凌偉，劉宇，王海興，等.電荷交換離子對柵極系統束流影響的數值研究［J］.航空動力學報，2009，24(8):53～56.

［12］熊家貴，王德武.離子引出的二維PIC－MCC模擬［J］.物理學報，2000，49(12):2420～2426.

［13］BIRDSALL C.K，LANGDON A B.Plasma Physics via Computer Simulation［M］.Mcgraw－Hill，1985.

［14］ROSENBERG D，WEHNER G K.Sputtering Yield for Low Energy He+，Kr+，and Xe+Ion Bombardment［J］.Journal of Applied Physics，1962，33(5):1842 ～1845.

［15］李娟，顧左，江豪成，等.氙離子火箭發動機補償柵極設計［J］.真空與低溫，2005，11(1):29～33.