電子轟擊式離子推力器放電腔結(jié)構(gòu)對(duì)等離子體特性影響的全粒子仿真研究

孫安邦,閆涵,李昊霖

(西安交通大學(xué)電氣工程學(xué)院,710049,西安)

電推力器主要應(yīng)用于深空探測(cè)器、通信衛(wèi)星及空間科學(xué)實(shí)驗(yàn)等方面[1],其中離子推力器因比沖高、壽命長(zhǎng)、控制精度高等優(yōu)點(diǎn)受到重視,被廣泛應(yīng)用于航天工程中。放電腔是電子轟擊式離子推力器產(chǎn)生并引出等離子體的關(guān)鍵部件,其結(jié)構(gòu)是決定離子生成效率、均勻性、單一性的關(guān)鍵,并對(duì)推力器的推力、比沖、效率和壽命有著重要影響[2]。放電腔內(nèi)產(chǎn)生等離子體的過程是:空心陰極發(fā)射高能初始電子與推進(jìn)劑中性粒子發(fā)生碰撞電離,生成等離子體,等離子體中的離子通過電場(chǎng)加速向放電腔出口運(yùn)動(dòng)并被引出到柵極,其結(jié)構(gòu)如圖1所示。在放電腔設(shè)計(jì)中需要考慮到工質(zhì)氣體進(jìn)氣口位置以及空心陰極位置這些結(jié)構(gòu)參數(shù)。調(diào)整進(jìn)氣口位置可以改變工質(zhì)氣體在腔內(nèi)的運(yùn)動(dòng)路程,進(jìn)而改變主要發(fā)生碰撞電離的區(qū)域,影響電離率;空心陰極的長(zhǎng)度會(huì)影響腔內(nèi)電子分布以及離子的回流程度,改變引出離子的分布以及效率。離子光學(xué)系統(tǒng)的運(yùn)行對(duì)等離子體密度有著一定的要求,離子均勻性會(huì)顯著影響柵極壽命[3]。因此,需要進(jìn)一步開展相關(guān)研究來分析不同因素對(duì)放電腔等離子體特性的影響。

圖1 電子轟擊式離子推力器放電腔結(jié)構(gòu)示意圖

由于放電腔體的封閉結(jié)構(gòu)給實(shí)驗(yàn)測(cè)量帶來了一定的局限性和困難,為了研究放電腔內(nèi)的電離過程和引出機(jī)理,近年來科研人員開展了大量的數(shù)值模擬工作。Brophy等基于能量守恒方程建立了離子推力器放電腔數(shù)值模型,確定了產(chǎn)生等離子體所需的能量,研究了改善推力器性能的方法[4]。Goebel等建立了自洽等離子體放電模型,計(jì)算了放電腔內(nèi)的等離子體密度、能量和電勢(shì),發(fā)現(xiàn)磁場(chǎng)對(duì)壁面附近的等離子體密度影響較大[5]。Wirz等建立了二維軸對(duì)稱混合模型,對(duì)二次電子和離子采用流體描述,模擬得到了腔內(nèi)磁場(chǎng)、密度分布和電流密度等參數(shù)[6]。Mahalingam等開發(fā)了二維粒子模擬結(jié)合蒙特卡羅的全粒子模型[7]。該模型考慮了放電腔內(nèi)的初始電子、二次電子、單價(jià)和雙價(jià)離子以及中性粒子5種主要粒子,并采用人為增大介電常數(shù)的方法加速計(jì)算[8],主要研究了粒子權(quán)重、時(shí)間步長(zhǎng)和磁場(chǎng)大小等因素對(duì)等離子體分布的影響。Stueber建立了放電腔三維粒子模型,主要研究了磁場(chǎng)對(duì)初始電子的影響[9]。蘭州空間技術(shù)物理研究所陳娟娟等采用網(wǎng)格粒子和蒙特卡羅碰撞(PIC/MCC)方法對(duì)LIPS-200離子推力器放電腔原初電子動(dòng)力學(xué)行為進(jìn)行了模擬,研究了電磁場(chǎng)對(duì)原初電子的影響以及等離子體的運(yùn)動(dòng)特性[10]。

目前,對(duì)于離子推力器放電腔的研究,主要采用的是混合模擬和全粒子模擬?；旌夏M相比全粒子模擬計(jì)算量較小,但是由于對(duì)部分粒子做了流體近似,因此不能得到局部粒子的微觀特性。對(duì)于大多數(shù)采用全粒子模擬的研究,主要是對(duì)放電腔內(nèi)工作特性進(jìn)行描述,或只針對(duì)仿真中的計(jì)算參數(shù)如粒子權(quán)重、時(shí)間和空間步長(zhǎng)對(duì)等離子體特性的影響,以及磁極結(jié)構(gòu)對(duì)推力器某一特性的影響[5-10],而進(jìn)氣口位置和空心陰極長(zhǎng)度對(duì)離子推力器放電腔的優(yōu)化設(shè)計(jì)有著重要的參考意義,且相關(guān)研究不充分,有必要進(jìn)行深入討論。因此,本文基于PIC/MCC方法,采用全粒子的二維軸對(duì)稱模型,仿真研究了不同工質(zhì)氣體進(jìn)氣口位置以及陰極位置對(duì)放電腔等離子體特性的影響,并對(duì)仿真結(jié)果進(jìn)行分析。

1 仿真模型及參數(shù)設(shè)置

1.1 放電腔模型

鑒于放電腔的軸對(duì)稱性,可以在二維柱坐標(biāo)下對(duì)其進(jìn)行模擬,本文的仿真區(qū)域長(zhǎng)10 cm,半徑5 cm,柵極半徑為4.5 cm,放電腔外包圍著3個(gè)環(huán)狀釤鈷永磁體,如圖2所示。

圖2 放電腔計(jì)算結(jié)構(gòu)示意圖

本文仿真采用靜電模型,通過泊松方程求解電勢(shì)

(1)

式中:φ為電勢(shì);ε0為真空介電常數(shù);e為元電荷;ni和ne分別表示離子和電子密度。在二維軸對(duì)稱柱坐標(biāo)下,式(1)可展開為

(2)

式中z和r分別表示軸向和徑向坐標(biāo)。

本文采用超松弛迭代法來求解式(2),放電腔內(nèi)電場(chǎng)通過電勢(shì)梯度得到,即

E=-φ

(3)

帶電粒子的運(yùn)動(dòng)遵循牛頓運(yùn)動(dòng)定律[11]

(4)

式中m、q、v分別為帶電粒子的質(zhì)量、電荷以及速度。本文采用Boris算法[12]求v和x。Boris算法首先通過半個(gè)時(shí)間步長(zhǎng)的電場(chǎng)加速,速度由vn-1/2更新為v-;接著在一個(gè)時(shí)間步長(zhǎng)內(nèi),粒子的旋轉(zhuǎn)速度由v-更新為v+;最后再經(jīng)過半個(gè)時(shí)間步長(zhǎng)的電場(chǎng)加速,速度更新為vn+1/2。具體計(jì)算公式如下

(5)

式中:上標(biāo)n表示當(dāng)前時(shí)刻;t′和s′分別定義為

(6)

(7)

根據(jù)求解的速度,粒子位置可更新為

xn+1=xn+vn+1/2Δt

(8)

仿真中主要考慮了電子與中性粒子之間的彈性、激發(fā)和電離碰撞,以及離子與中性粒子之間的彈性碰撞和電荷交換碰撞。工質(zhì)氣體選擇氙氣,并在仿真的最開始通過直接蒙特卡羅碰撞(DSMC)的方法進(jìn)行模擬,得到其空間密度和速度分布特征,并在之后帶電粒子的計(jì)算中設(shè)置為背景氣體。

為了簡(jiǎn)化計(jì)算,本模型采用空碰撞的方法[11]使總碰撞率恒定。最大碰撞頻率為

Ri=max(nneu)max(σiv)

(9)

式中:nneu為背景粒子數(shù)密度;σi為目標(biāo)粒子與背景粒子的碰撞截面;v為目標(biāo)粒子速度。帶電粒子與氙氣的碰撞截面數(shù)據(jù)取自Lxcat網(wǎng)站[13]。

在Δt時(shí)間內(nèi),目標(biāo)粒子與背景粒子發(fā)生空碰撞的概率可以表示為

P=1-exp(RiΔt)

(10)

則每個(gè)時(shí)間步長(zhǎng)內(nèi)發(fā)生碰撞的最大次數(shù)Ncoll=NtotalP,其中Ntotal為模擬中的總粒子數(shù)。

通過空碰撞可以避免對(duì)每個(gè)粒子都進(jìn)行碰撞的判斷和處理,從而減小了計(jì)算開銷。對(duì)于可能發(fā)生碰撞的粒子,結(jié)合蒙特卡羅方法判斷其碰撞類型[14]。

由于全粒子模擬受德拜長(zhǎng)度的限制,網(wǎng)格數(shù)多且時(shí)間步長(zhǎng)小,需要采取一定的加速方法來減少計(jì)算量。本模型采用了Taccogna提出并已驗(yàn)證可行性的相似設(shè)計(jì)方法[15],核心思想是人為將放電腔尺寸乘以ξ-1,其中ξ是縮比系數(shù)且大于1,同時(shí)遵循相似原理,保持主要特征量不變,從而減小計(jì)算域中的網(wǎng)格數(shù)?？s放后各參數(shù)變化情況如表1所示。

表1 相似設(shè)計(jì)方法參數(shù)變化

此外,隨著時(shí)間的增加,計(jì)算域內(nèi)粒子數(shù)急劇增多,為了保證計(jì)算效率,引入宏粒子的概念,宏粒子權(quán)重表示真實(shí)粒子數(shù)。由于權(quán)重過大或過小都會(huì)影響仿真,所以模型中采用粒子自適應(yīng)權(quán)重(APM)的方法來動(dòng)態(tài)調(diào)整粒子權(quán)重[16-17],實(shí)現(xiàn)流程如下。

(1)通過定義每個(gè)網(wǎng)格內(nèi)的理想粒子數(shù)Nppc獲得第i個(gè)網(wǎng)格內(nèi)粒子的理想權(quán)重

(11)

式中:ni為第i個(gè)網(wǎng)格粒子數(shù)密度;Vcell為網(wǎng)格體積。

(2)當(dāng)權(quán)重小于(2/3)w時(shí)粒子被合并以增大權(quán)重,當(dāng)大于1.5w時(shí)粒子被分裂。

(3)對(duì)于需要合并的粒子,選擇與其距離最近的粒子合并。粒子間的距離d為

d2=(xi-xj)2+η2|vi-vj|2

(12)

式中η為距離與速度的比例因子,模擬中設(shè)置為1 ps。新形成的粒子速度與原始粒子有一個(gè)隨機(jī)的偏移量,權(quán)重為兩個(gè)原始粒子權(quán)重之和,位置可由下式得到

(13)

(4)需要分裂的粒子將被分裂為兩個(gè)新粒子,它們的速度和位置與原始粒子相同,權(quán)重為原始粒子的一半。整體仿真流程如圖3所示。

圖3 PIC/MCC仿真流程圖

1.2 參數(shù)設(shè)置

本模型的具體參數(shù)設(shè)置如表2所示。放電腔外加電勢(shì)一般為上千伏,但為了便于分析,將所有電勢(shì)均減去陰極電勢(shì)1 200 V,從而保證陰極電位為0。實(shí)驗(yàn)中測(cè)得空心陰極發(fā)射的初始電子能量大約在2～3 eV[18],仿真中取2 eV;中性粒子的柵極透過率設(shè)置為0.3,離子的柵極透過率為0.86。

表2 仿真參數(shù)設(shè)置

2 仿真結(jié)果分析

2.1 進(jìn)氣口位置對(duì)等離子體特性的影響

本文選擇了圖4所示的4種進(jìn)氣口位置,其中雙進(jìn)氣口由陰極進(jìn)氣口和主進(jìn)氣口構(gòu)成。

(a)位置1

因?yàn)楣べ|(zhì)氣體是背景氣體,所以在運(yùn)行等離子體程序前,需要采用DSMC方法預(yù)先模擬4種工況下中性粒子的運(yùn)動(dòng)情況。

中性粒子數(shù)密度分布如圖5所示:進(jìn)氣口處的粒子數(shù)密度最大,單進(jìn)氣口約為2.1×1019m-3,雙進(jìn)氣口在1×1020～4×1020m-3范圍內(nèi),高出單進(jìn)氣口一個(gè)數(shù)量級(jí);雙進(jìn)氣口工況下放電腔出口處中性粒子數(shù)密度最小,在2.3×1019～3×1019m-3區(qū)間內(nèi),單進(jìn)氣口工況下側(cè)壁面中性粒子數(shù)密度最小,為1.2×1018m-3。

(a)位置1

以上述中性粒子分布為背景,計(jì)算得到離子數(shù)密度在放電腔內(nèi)的分布,如圖6所示。陽極壁面的離子數(shù)密度最低,約為1×1012m-3,這是由于磁場(chǎng)對(duì)電子的約束作用減少了電子對(duì)陽極壁面的撞擊,碰撞電離主要發(fā)生在腔體中部,因此離子也主要聚集在此區(qū)域。由圖5a知,進(jìn)氣口位置1對(duì)應(yīng)的腔體上游中性粒子密度低,因此電離產(chǎn)生的離子主要集中在放電腔下游,隨著主進(jìn)氣口位置靠近下游,離子高密度區(qū)也逐漸向下游出口處移動(dòng)。

放電腔出口處離子數(shù)密度分布如圖7所示。由圖可以看到:雙進(jìn)氣口的離子數(shù)密度峰值較單進(jìn)氣口偏離了對(duì)稱軸,這是雙進(jìn)氣口工況的中性粒子高密度區(qū)偏向陽極壁面所導(dǎo)致的;進(jìn)氣口位置4整體離子數(shù)密度最高,進(jìn)氣口位置2和3的密度相差不大,進(jìn)氣口位置1最小;4種工況下最大離子數(shù)密度為7.9×1017m-3,平均離子數(shù)密度在2.6×1017～4.1×1017m-3之間,與蘭州空間技術(shù)物理研究所結(jié)果的數(shù)量級(jí)相吻合[19];在徑向位置增大至4.5 cm左右時(shí),由于陽極后壁面的存在,離子數(shù)密度趨近于0,出口處離子密度變化趨勢(shì)與文獻(xiàn)[3]基本一致。

圖7 不同進(jìn)氣口位置下腔體出口處離子數(shù)密度隨徑向位置的變化

為了保證離子推力器引出束流的穩(wěn)定性和準(zhǔn)直性,需要保證腔體內(nèi)尤其是柵極附近等離子體參數(shù)的均勻性。用一個(gè)區(qū)域網(wǎng)格點(diǎn)上的密度與平均密度比值的標(biāo)準(zhǔn)差來衡量該區(qū)域粒子分布的均勻度,如下式所示

(14)

通過計(jì)算得到腔內(nèi)的平均電離率和出口處離子均勻度,如圖8所示。雙進(jìn)氣口工況的平均電離率高于單進(jìn)氣口,進(jìn)氣口位置1的電離率為9.67×1023m-3,進(jìn)氣口位置2最大,為1.1×1024m-3;雙進(jìn)氣口工況的離子均勻性差于單進(jìn)氣口,進(jìn)氣口位置1的S=0.65,進(jìn)氣口位置3的離子均勻性最差,S=0.81。

圖8 不同進(jìn)氣口位置的平均電離率和離子均勻度

表3對(duì)比了推力器的宏觀性能,發(fā)現(xiàn)雙進(jìn)氣口工況的性能優(yōu)于單進(jìn)氣口,其中進(jìn)氣口位置4最優(yōu)。4種工況下計(jì)算的推力在0.094～0.112 N之間,在文獻(xiàn)[20]中所報(bào)道的實(shí)驗(yàn)數(shù)據(jù)范圍內(nèi)。

表3 不同進(jìn)氣口位置下推力器的宏觀性能

綜合來看,進(jìn)氣口位置4既能保證較高的電離率1.07×1024m-3,得到較均勻的離子分布,S=0.69,又能擁有最優(yōu)的推力器宏觀性能。

2.2 陰極位置對(duì)等離子體特性的影響

放電腔內(nèi)的初始電勢(shì)可以空心陰極前端為界分為兩部分,一部分是由于陰極低電勢(shì)使電場(chǎng)指向陰極,另一部分是由于屏柵低電勢(shì)使電場(chǎng)指向屏柵,這就可能導(dǎo)致電離產(chǎn)生的離子回流到陰極并被還原為中性粒子,影響推力器效率,因此在設(shè)計(jì)推力器時(shí)就需要考慮陰極位置對(duì)等離子體特性的影響。

基于2.1小節(jié)的結(jié)果,將進(jìn)氣口位置4設(shè)定為本節(jié)背景氣體初始條件,并選擇陰極長(zhǎng)度lc為3、4、5、6和7 cm這5個(gè)工況進(jìn)行研究。

不同的陰極位置會(huì)直接影響電子分布,進(jìn)而影響出口離子均勻性以及電離率,電子數(shù)密度分布如圖9所示。由圖可以看出:在磁場(chǎng)的約束作用下,只有少數(shù)電子撞擊到永磁體附近的陽極壁面,從而減弱了對(duì)壁面的腐蝕,電子數(shù)密度在3.5×1014～1.9×1015m-3之間;隨著陰極長(zhǎng)度的增加,放電腔上游電子數(shù)密度逐漸減小,下游電子數(shù)密度逐漸增大。

(a)lc=3 cm

表4為不同陰極位置下腔內(nèi)的平均電子密度和均勻度。增加陰極長(zhǎng)度會(huì)降低平均電子密度,極長(zhǎng)度過短,又會(huì)影響電子在腔內(nèi)的均勻性,lc在4～6 cm之間是相對(duì)較合適的選擇。

表4 不同陰極位置下腔內(nèi)平均電子密度和均勻度

但陰

圖10為不同陰極長(zhǎng)度下腔內(nèi)離子數(shù)密度分布。由于離子分布主要由電子分布決定,因此其隨陰極長(zhǎng)度的變化趨勢(shì)與電子基本一致。由圖可以看到,陰極長(zhǎng)度為3、4和5 cm的工況中陰極上方區(qū)域、陽極前壁面、上游側(cè)壁面離子數(shù)密度較高,回流現(xiàn)象較明顯。

(a)lc=3 cm

為了對(duì)回流現(xiàn)象進(jìn)行更精確的研究,選取陰極上方區(qū)域的離子數(shù)密度進(jìn)行分析,如圖11所示。由于陰極長(zhǎng)度不同,對(duì)軸向坐標(biāo)進(jìn)行歸一化處理,即

(15)

式中:z*為相對(duì)軸向坐標(biāo);z為實(shí)際軸向坐標(biāo)。

從圖11可以看出,隨著軸向位置的增加,陰極上方區(qū)域離子數(shù)密度整體呈上升趨勢(shì),而陰極長(zhǎng)度的增加會(huì)減小陰極上方區(qū)域的離子數(shù)密度,但當(dāng)lc增加到6 cm之后,離子數(shù)密度下降趨勢(shì)不再明顯。

圖11 不同陰極位置下陰極上方區(qū)域離子數(shù)密度隨軸向位置的變化

圖12為放電腔出口處離子數(shù)密度分布,可以看到:隨著徑向位置的增加,出口處離子密度先保持在一定的范圍上下波動(dòng),而后當(dāng)徑向坐標(biāo)接近4.5 cm時(shí),離子數(shù)密度逐漸減小并趨近于0;5種工況下離子平均密度在2.9×1017～8.9×1017m-3范圍內(nèi),與文獻(xiàn)[19]的數(shù)據(jù)基本吻合;隨著陰極長(zhǎng)度的增加,出口處離子數(shù)密度逐漸增大,同樣當(dāng)lc增加到6 cm之后,離子數(shù)密度基本不再增大。

圖12 不同陰極位置下腔體出口處離子數(shù)密度隨徑向位置的變化

圖13為5種工況下腔內(nèi)平均電離率和出口處離子均勻度,可知lc=7 cm的平均電離率最小,為8.74×1023m-3,lc=3 cm的平均電離率最大,為2.54×1024m-3;lc=5 cm的離子均勻性最好,S=0.63,lc=3 cm的離子均勻性最差,S=0.8。

圖13 不同陰極位置的平均電離率和離子均勻度

通過對(duì)比表5給出的推力器宏觀性能,發(fā)現(xiàn)陰極長(zhǎng)度的增大可以提升推力器的宏觀性能。

表5 不同陰極位置下推力器的宏觀性能

綜合來看:雖然lc=3 cm工況的電離率高出其他工況近1倍,但是回流現(xiàn)象最嚴(yán)重,出口處離子均勻性最差,推力器宏觀性能最差;lc=7 cm工況的推力器宏觀性能最優(yōu),但是其電離率最小,離子分布均勻性較差;lc=4～5 cm的工況內(nèi),在保證較好離子均勻性的同時(shí),平均電離率相對(duì)較高,且回流現(xiàn)象不嚴(yán)重,推力器宏觀性能較好。

3 結(jié) 論

本文基于PIC/MCC方法,建立了電子轟擊式離子推力器放電腔的全粒子二維軸對(duì)稱模型,仿真發(fā)現(xiàn)進(jìn)氣口位置和陰極位置會(huì)對(duì)等離子體分布、推力器的推力、比沖和放電損耗產(chǎn)生影響,可為放電腔的實(shí)際設(shè)計(jì)提供一定的參考。通過分析,得出以下結(jié)論。

(1)雙進(jìn)氣口工況的平均電離率和推力器性能優(yōu)于單進(jìn)氣口,但放電腔出口處離子均勻性相對(duì)較差,綜合等離子體分布和推力器宏觀性能,進(jìn)氣口位置4工況的性能最優(yōu)。

(2)陰極長(zhǎng)度的增加可以減弱離子回流現(xiàn)象,推力器宏觀性能參數(shù)也得到提升,但陰極過長(zhǎng)又會(huì)降低電離率;陰極長(zhǎng)度在4～5 cm區(qū)間內(nèi)時(shí),出口處離子分布均勻性最好,電離率較高,但推力器的宏觀性能會(huì)有所降低。