考慮電荷交換的柵極區(qū)離子流數(shù)值模擬

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

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

1 引言

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

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

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

2 離子運動計算模型設(shè)置

2.1 計算域和邊界

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

圖1 計算域二維簡圖

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

表1 柵極幾何參數(shù)

2.2 CEX 的產(chǎn)生與加速柵極腐蝕

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

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

碰撞概率為

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

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

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

式中 vmax為計算域內(nèi)模擬離子的最大速度;Δz為z方向的網(wǎng)格長度。

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

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

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

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

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

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

3 計算流程

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

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

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

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

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

(3)CEX在離子光學系統(tǒng)中的產(chǎn)生

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

4 模擬結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

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

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

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

圖9 加速柵極表面濺射產(chǎn)額分布

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

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

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

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

圖13 加速柵極孔壁濺射產(chǎn)額分布

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

5 結(jié)論

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

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

［2］吳漢基，蔣遠大，張志遠.電推進技術(shù)的應(yīng)用與發(fā)展趨勢［J］.推進技術(shù)，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］鐘凌偉，劉宇，王海興，等.電荷交換離子對柵極系統(tǒng)束流影響的數(shù)值研究［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］李娟，顧左，江豪成，等.氙離子火箭發(fā)動機補償柵極設(shè)計［J］.真空與低溫，2005，11(1):29～33.