感應式脈沖推力器中等離子體加速數值研究?

成玉國 夏廣慶

1)(中國人民解放軍91550部隊91分隊,大連 116023)

2)(大連理工大學,工業裝備結構分析國家重點實驗室,大連 116024)

1 引 言

目前預研的具有較大應用前景的無電極燒蝕、無羽流中和器的感應式脈沖等離子體推力器中,一類以FARADAY為代表[1],其將無電極等離子體源(如螺旋波源等[2,3])作為預電離級,在源室下游進一步利用脈沖電流圈加速;另一類以PIT為代表[4]的推力器則直接利用特殊構型的線圈,如圖1(a)所示,電離其表面的中性氣體在等離子體中形成環形電流片,與線圈排斥產生推力,實現中性氣體的電離和加速.這種推力器中,中性氣體通過脈沖流量閥的引導,均勻地分布在線圈表面,電容器能量經線圈激勵產生電磁場,傳遞至氣體,使其電離和加速(圖1(b)).后一類脈沖工作方式具有推力器結構簡單、方便調節、比沖高等優點,對其工質加速過程和推力產生機理的研究仍處于不斷深化和發展之中.

圖1 (a)PIT實物圖;(b)PIT結構簡圖[5]Fig.1.(a)The picture of the actual object of the PIT;(b)the Simpli fied con figuration of the PIT[5].

Jahn[6]提出上述第二類推力器推進應用中脈沖氣體電離和加速的兩個主要問題:1)氣體擊穿特征時間與脈沖電場回路放電時間的匹配問題;2)如何在氣體運動遠離線圈表面之前實現絕大部分脈沖能量的注入,以提高電離和加速效率的問題,并進一步提煉出線圈等離子體解耦距離這一影響加速的主要參數.

針對上述問題,研究人員分別發展了線圈等離子體脈沖感應耦合機電模型和等離子體加速的磁流體力學模型.機電模型實現了線圈回路和等離子體負載回路通過感應線圈的耦合[5].Polzin和Choueiri[7]總結出動態匹配系數α,將回路放電和等離子體加速過程聯系,分析了等離子體的解耦距離.在近似正弦的脈沖電流中,等離子體在前1/4周期內達到解耦距離時,等離子體可實現高效率的電離和加速.

二維磁流體數值研究中,Mikellides等[8]利用MACH2計算發現了在氣體質量＞3 mg時,推力器的效率保持穩定.這一結果與實驗中觀測到的臨界質量現象一致.在一些實驗中,測量結果表明PIT推進效率隨著能量的變化而改變[9],且Mikellide等[10]將MACH2與電路模型結合的計算進一步表明了超過臨界質量之后效率的下降是由于內能的沉積所致.這些矛盾的結論表明,感應式脈沖加速推力效率這一問題值得進一步的探索.已有的模擬中,對于影響推力產生的等離子體運動過程、運動參數的空間分布以及磁場強度變化對等離子體加速和電離過程的影響研究較少.

本文研究感應式脈沖等離子體的電離和加速過程,計算在單脈沖、短放電周期(10μs)條件下等離子體流場參數分布及其特性,并著重分析脈沖強度變化時等離子體的加速特性和推進性能的變化情況.

2 感應式脈沖等離子體流動控制方程

2.1 物理模型

在μs脈沖作用下,等離子體中的感應電磁場變化較快,電子與重粒子之間的耦合強烈,形成整體的等離子體團運動,采用文獻[11—13]處理等離子體組分的方法,對等離子體在脈沖電磁場激勵下的運動采用單流體磁流體力學描述,認為電子和離子處于熱力學平衡狀態.其連續、動量和能量方程分別為:

連續方程

式中,ρ為等離子體密度,V為等離子體速度;動量方程

式中,B為等離子體流場中的磁場強度,p為等離子體壓力,J為等離子體電流,J=?×B/μ0,μ0為真空磁導率;

能量方程

式中,et為總比能,,q為等離子體熱流,σe為等離子體電導率,f(p,ρ)為真實氣體狀態函數,其計算方法在2.4節中討論,Qrad,tn為輻射能量損失,采用文獻[14]的處理方法.

對等離子體流場中的電磁場運動采用參考文獻[15]的處理方法,用如下形式的磁擴散方程關聯磁場和速度場:

式中,?×(V×B)=?·(V B?BV).

為了解決以上方程組中的無黏項Jacobi矩陣不滿秩和磁場散度由于離散造成的不為零的問題,本文采用文獻[16,17]給出的雙曲型散度清除代數方程,其已經在經典算例中得到了驗證.方程如下:

式中,ψ為廣義拉格朗日乘子;ch,cb分別為控制?·B誤差的對流和耗散系數.

2.2 廣義拉格朗日乘子形式的磁流體動力學(MHD)方程組

由以上方程,得到廣義拉格朗日乘子形式的MHD方程組的二維軸對稱形式為

式中,μ為黏性系數.

以上方程組由等離子體的狀態方程(11)封閉,

式中,ne為電子數密度;na為中性粒子數密度;j=1,2···,表示離子電離階數;T為等離子體溫度,與電子溫度Te相同;kB為玻爾茲曼常數.

2.3 定解條件

2.3.1 氣動邊界條件

圖1所示推力器的計算簡化結構如圖2所示,其主要由電路和等離子體部分組成,電容器經線圈激勵在其表面產生高頻振蕩電磁場,并逐步滲透至等離子體內部.文中計算主要關注等離子體部分,電容器放電對等離子體的影響在計算中體現為:在等離子體邊緣(線圈表面,AD邊)施加時變磁場.氣體流經閥門,在線圈表面沉積,在脈沖電流線圈的作用下,加熱并電離,向近似真空的環境中膨脹加速.等離子體的入流/出流、壁面以及軸對稱邊界等三類邊界條件根據文獻[18]的討論給出.

圖2 (網刊彩色)感應式脈沖等離子體計算域構型Fig.2.(color online)Calculation con figuration of the inductive pulsed plasma.

2.3.2 磁場邊界條件

入口處磁場的條件由電容器的放電頻率和強度決定,假設其為正弦變化,線圈電流變化為J(t),徑向上的磁感應強度Br隨時間變化為

即Br(t)=B0sin(ωt),B0為峰值強度,ω=2πf,f為放電頻率,文中設置放電周期為10μs,router,rinner為分別為線圈外徑和內徑.線圈表面的磁場主要在徑向,在角向為零、徑向均勻的情況下,忽略軸向磁場,此時表面的磁場散度為零.

文中計算的放電波形為正弦變化,根據COMSOL軟件的計算,對于2.3節中所述的線圈的幾何構型,單周期正弦脈沖能量和峰值磁場強度B0的關系如圖3所示.

圖3 單脈沖正弦能量Esingle-pulse與峰值磁場強度B0關系曲線Fig.3. The relationship between single sinusoidal pulse energy Esingle-pulseand peak field intensity B0.

壁面處的磁場采用無穿透條件.流場下游區域足夠大,認為在脈沖作用的時間內,尚未影響下游的計算區域,邊界處的軸向磁場Bz(t)采用外推條件.

實際計算中,發現在脈沖作用周期內,在電流圈斥力的作用下,等離子體軸向運動距離在1—2個線圈半徑的長度.因此,為了突出主要運動區域,研究其中的流動機理,文中徑向起算點為A,B點,軸向將AB線段設置為壁面條件,暫不考慮拐角處熱膨脹對流動的影響.

對于文中的非定常MHD流動,空間格式應用M-AUSMPW+[15,19],采用結構網格離散,在時間上采用具有真實垂直深度性質的三階Runge-Kutta法.文中計算的構型中,線圈內徑為0.05 m,外徑為0.2 m,軸向方向取5倍內徑長度.文中對單脈沖周期中等離子體的電離、加熱以及加速過程開展研究.

2.4 熱力學性質

計算采用Ar氣,高溫情況下,其將會出現三階以上電離.本文考慮氬氣的五階電離.文中采用的簡并度及其能級數值由美國標準與技術研究院數據庫(National Institute of Standards and Technology,NIST)給出.

混合組分的比焓表達式為[20]

式中,zj,Ej為各能級配分函數和電離能.

比焓與比內能的關系為

由(13)式,等離子體的比焓由重粒子分量、電子分量、電子激發能分量以及電離分量構成.比焓計算中對文獻中報道較少的溫度范圍(2.5—5.0 eV)進行了拓展,如圖4所示,在低溫范圍內(0—2.5 eV),結果與文獻符合[20].

圖4 100.0 Pa條件下Ar等離子體比焓隨溫度變化Fig.4.The speci fic enthalpy of Argon plasma varies with Tefor P=100.0 Pa.

3 計算結果與分析

3.1 脈沖作用周期內等離子體的運動分析

圖5(a)—(c)中給出了在峰值磁場強度為0.1 T情況下,單脈沖周期內等離子體流場參數隨時間的變化.可以看出,隨著時間的推進,等離子體的能量和速度逐漸增加,且在全部脈沖作用時間內等離子體的速度始終在增加.

由等離子體在各時刻的狀態參數分布可以看出,在脈沖電流的作用下其運動圖景大致描述如下.

放電初始階段,在感應線圈產生的時變磁場作用下,等離子體中產生感應電磁場,電場強度達到氣體電離閾值后,中性氣體電子被剝離、加速,并進一步撞擊中性氣體.線圈表面等離子體迅速被加熱、升溫并向外膨脹,在瞬變磁場的持續作用下,表面的高速、相對高密度氣體向前方運動(圖5(a),(b)),同時,這一過程伴隨著線圈的磁場能量在等離子體中的逐步滲透,且流場中的磁場強度及其影響范圍也在逐漸增大.

隨著放電和加速過程的持續進行,達到1/2周期以后,線圈電流的方向發生改變.氣壓足夠大或者峰值磁場強度較小時,反向作用周期內表面殘余氣體可獲得進一步加熱,與前方等離子體團的電流片異向,等離子體中前后異向的電流片相互吸引,從而降低前向等離子體團的速度,并減弱等離子體的總體加速效果.

圖5(c)給出了等離子體的磁場隨時間的變化.脈沖電流超過1/4周期后,電流方向不變,但其變化率與前1/4周期相反.這段時間內,中性氣體在脈沖電磁場作用下進一步電離,表面的高溫等離子體向前緣移動,但氣體密度較高的情況下,線圈表面將滯留帶電粒子和中性氣體的混合物.放電脈沖進入后半周期后,這部分氣體受負變化率的影響產生反向等離子體電流環,這種新產生的電流環一方面對前緣電流環具有吸引作用,另一方面阻止了激勵磁場進一步向前緣的滲透.

3.2 不同單脈沖能量下等離子體加速過程研究

圖6給出了線圈峰值磁場強度為0.45 T時等離子體速度峰值時刻流場的參數分布(密度、軸向速度、徑向磁場、溫度、電子數密度和電流密度).在強脈沖電磁場的作用下,將有更多的電磁場能量用于加速等離子體.從圖6的結果可以看到,此時流場運動與壁面黏滯阻力特征時間差異較大,尤以密度和電流的分布較為明顯,表明在較大的磁場強度下,壁面黏性對等離子體造成的影響減小.

圖6中,線圈表面等離子體密度與背景氣體幾乎相同,從電子數密度的分布圖中可以看出,此時表面附近的電子數密度較前緣峰值降低了約1—2個量級(圖6中的密度分布),但由于這一區域獲得的電磁場能量最大,等離子體溫度為流場的高值區域,如圖6中的溫度分布所示.在較高磁場強度下,放電一段時間后,表面的等離子體粒子數密度是極低的.因此,雖然溫度較高,但基本不具備進一步產生可觀沖量的價值.

圖5 (網刊彩色)0.1 T時不同時刻等離子體流場參數變化 (a)密度隨時間變化的分布;(b)速度隨時間變化的分布;(c)等離子體磁感應強度隨時間變化的分布Fig.5.(color online)The variations of the flow parameters at different time for 0.1 T:(a)The density distribution as the time is varied;(b)the velocity distribution as the time is varied;(c)the plasma inductive field intensity distribution as the time is varied.

圖6 (網刊彩色)峰值磁場強度為0.45 T等離子體流場參數分布Fig.6.(color online)The plasma parameter distribution for peak field intensity of 0.45 T.

圖7 (網刊彩色)峰值磁場強度為0.45 T時不同時刻等離子體速度的分布Fig.7.(color online)The plasma velocity distribution varies with the time for the peak field intensity of 0.45 T.

從圖6中磁場和等離子體電流片的分布可看出,磁場的有效作用范圍受到等離子體密度分布影響,加速過程中,遠離線圈表面時,等離子體中磁場逐漸減弱.

圖7中給出了不同時刻流場中等離子體速度隨時間的變化(從零加速至峰值以及之后的減速).從圖中可以看出,速度相對于3.1小節中的低強度磁場提升了約1—2個量級,其最大值出現于放電的1/4與1/2周期之間.隨著放電時間的推移,等離子體團速度逐漸增加,但受到磁場與流體運動相位滯后以及慣性的影響,越過電流峰值之后,等離子體速度在1/4周期后達到最大.

從圖7可以看出,在較高的脈沖能量情況下,由推力器幾何構型決定的線圈等離子體解耦距離對等離子體加速影響顯著,表明了高脈沖輸入能量可獲得較高的加速效率,但加速過早地達到解耦距離,也限制了能量的進一步有效利用.

圖8給出了根據電感計算理論得到的線圈結構與解耦距離的關系[21,22],對于文中內徑0.05 m、外徑0.2 m的線圈,其解耦距離約為0.04 m,這與0.45 T峰值磁場情況下等離子體在1/4周期內運動距離相近,表明在這一脈沖能量下,等離子體的運動開始達到解耦距離.

圖8 (網刊彩色)解耦距離與線圈內外徑的關系Fig.8.(color online)The relationship between the decouple distance and the coil outer-inner diameter.

圖9 不同時刻等離子體放電參數沿軸向的分布(0.45 T)Fig.9.The axial plasma discharge parameters vary with the time for peak intensity of 0.45 T.

圖10 不同時刻等離子體密度沿軸向的分布(0.45 T,徑向截面平均)Fig.10.The axial distribution of the plasma density as the time is varied(0.45 T,radial cross section averaged).

圖11 (網刊彩色)不同磁場強度的峰值速度時刻電流密度分布Fig.11.(color online)The current density distribution of different field intensity at the time of peak velocity.

由上文中峰值強度0.1 T和0.45 T的結果對比可以看出,磁場強度較低的情況下,等離子體速度低,加速距離短,脈沖作用時間內可始終獲得加速;但電磁場能量低,等離子體的能量以及感應強度有限,加速性能提高有限.磁場強度增大后,一部分能量增加等離子體內能,另一部分增加等離子體團的動能,于是等離子體中誘導電流片吸收的能量增加[23],從而導致等離子體團推進速度增加.主等離子體的運動距離超出解耦距離后,受脈沖電流及磁場的影響減弱,同時部分能量傳遞通過傳導以及輻射等方式傳遞給中性氣體,且隨著磁場強度的增強,誘導電流片速度進一步增加,與電流線圈解耦的時間縮短.

圖9給出了不同時刻等離子體放電參數沿軸向的分布(參數對徑向進行了平均).隨時間變化的曲線表明,等離子體不斷向前推進過程中,受影響的氣體區域不斷擴大.可以看出磁場的影響區域明顯大于等離子體流場的擾動區域,這表明磁場在等離子體中有效滲透,有利于等離子體的電離.

圖10給出了不同時刻等離子體數密度沿軸向的分布,密度對徑向進行了平均.結果表明,隨著時間的發展,等離子體數密度逐漸增加,高階電離粒子數密度組分增大,線圈表面附近等離子體幾乎完全電離.由上文的結論,密度的這一分布結果主要由于表面是加熱的主要區域造成.

圖11和圖12給出了不同峰值磁場強度下最大速度時刻等離子體流場中電流密度和流體密度的分布.從圖中可以看到,磁場強度越高,等離子體獲得的能量越多,其加速越顯著,磁場強度提高的極限是在放電經過峰值時.進入負半周期之后,由于沒有中性氣體的補充,等離子體團遠離耦合距離,這段放電時間的能量幾乎得不到合理的利用.因此,提高磁場強度可能會造成能量利用效率增長的緩慢.從圖13速度的變化趨勢可看出,磁場強度增大后,速度幅值增量逐漸減小.

圖12 (網刊彩色)不同磁場強度的峰值速度時刻密度分布Fig.12.(color online)The plasma density distribution of different field intensity at the time of peak velocity.

圖13 不同峰值磁場強度的最大激勵速度變化曲線Fig.13.The maximum velocity varies as the peak filed strength is changed.

3.3 比沖和推力效率分析

圖14給出了感應式脈沖等離子體比沖和推力效率隨脈沖能量的變化,結果表明,在前半周期運動能夠達到解耦距離的情況下,等離子體的比沖和效率將有較大的提升.比沖從低脈沖(0.1 T)情況下約700 s提高到高脈沖能量(0.55 T)的約為2000 s.感應式加速中能量的提高對推進性能的提升顯著,顯示出感應式脈沖加速的優勢,尤其是在四分之一周期電離充分的情況下,等離子體的加速效果明顯.

圖14 比沖和推力效率隨單脈沖能量的變化Fig.14. The speci fic impulse and thrust efficiency varies as the pulse energy is changed.

4 結果與討論

本文建立了二維軸對稱的感應式脈沖等離子體作用的磁流體力學模型,針對計算中出現的高溫情況,對Ar等離子體的熱力學性質進行了計算,著重分析了不同脈沖電流能量情況下等離子體的加速特性、磁場形態的變化.

計算中發現在峰值較低的磁場強度(～0.1 T)情況下,電容器的放電能量主要用于中性氣體的電離以及氣動加熱,造成感應電磁場對粒子的加速效果有限,等離子體吸收的能量不足以支持高階電離,從而不能形成強度較高的感應電流片,與線圈之間的斥力較弱,因此,等離子體在低強度脈沖電流作用周期內,速度始終較低.

當峰值磁場強度增大(＞0.45 T)時,推進工質以及壓力不變的情況下,更高比例的輸入脈沖能量將被用于加速等離子體,劇烈變化并且幅值較大的電磁場有利于等離子體重粒子和電子之間的耦合,提高了等離子體感應電流片的電流密度和強度,可在等離子體-線圈之間產生較大的斥力,增加等離子體的速度,使得脈沖加速的效果有較大的提升.

磁場強度較高的情況下,等離子體速度可加速至20 km/s以上,從文中的計算結果可以看出,等離子體運動達到解耦距離的時間縮短,同時伴隨著等離子體高階電離、能量輻射、磁場影響迅速擴大等非線性過程,且隨著強度的不斷增大,等離子體速度提高,這一時間逐漸減小.獲得這段時間內各種物理過程對能量吸收和加速的影響,以及等離子體的解耦距離與磁場之間的關系,是下一步提高感應式脈沖等離子體推進性能的重要研究方向.

[1]Choueiri E Y,Polzin K A 2006J.Prop.Power223

[2]Cheng Y G,Cheng M S,Wang M G,Li X K 2014Acta Phys.Sin.63035203(in Chinese)[成玉國,程謀森,王墨戈,李小康2014物理學報63035203]

[3]Cheng Y G,Cheng M S,Wang M G,Li X K 2014Chin.Phys.B23105202

[4]Mikellides P G,Neilly C 2007J.Prop.Power231

[5]Polzin K A 2011J.Prop.Power273

[6]Jahn R G 1968Physics of Electric Propulsion(New York:McGraw-Hill)p268

[7]Polzin K A,Choueiri E Y 2006IEEE Trans.Plasma Sci.343

[8]Mikellides P G,Villarreal J K 2007J.Appl.Phys.10210

[9]Dailey C L,Lovberg R H 1993NASACR-1993-191155

[10]Mikellides P G,Turchi P J,Roderick N F 2000J.Prop.Power165

[11]Xie Z H 2013Ph.D.Dissertation(Changsha:National University of Defense Technology)(in Chinese)[謝澤華2015博士學位論文(長沙:國防科學技術大學)]

[12]Polzin K A,Sankaran K,Ritchie A G,Reneau J P 2013J.Phys.D:Appl.Phys.46475201

[13]Tian Z Y 2008Ph.D.Dissertation(Changsha:National University of Defense Technology)(in Chinese)[田正雨2008博士學位論文(長沙:國防科學技術大學)]

[14]Li X K 2011Ph.D.Dissertation(Changsha:National University of Defense Technology)(in Chinese)[李小康2011博士學位論文(長沙:國防科學技術大學)]

[15]Kim S,Soogab L,Kyu H K 2008J.Comput.Phys.2278

[16]Dedner A,Kemm F,Kr?ner D,Munz C D,Schnitzer T,Wesenberg M 2002J.Comput.Phys.1752

[17]Dedner A 2003Ph.D.Dissertation(Freiburg im Breisgau:Alert Ludwigs-Universit?t Freiburg)

[18]Yan C 2006Computational Fluid Dynamic(Beijing:Beihang University Press)pp254–258(in Chinese)[閻超2006計算流體力學方法及其應用(北京:北京航空航天大學出版社)第254—258頁]

[19]Kim K H,Kim C 2005J.Comput.Phys.2082

[20]Chen X 2009Thermal Plasma Heat Transfer and Flow(Beijing:Science Press)pp48–59(in Chinese)[陳熙 2009熱等離子體傳熱與流動(北京:科學出版社)第48—59頁]

[21]Kalantarov P L,Neitlqy L A(translated by Chen T M,Liu B A,Luo Y L,Zhang Y H)1992 Inductance Calculation Handbook(Beijing:Mechanism Industry Press)pp1–9(in Chinese)[卡蘭塔羅夫P L,采伊特林L A 著 (陳湯銘,劉保安,羅應力,張奕黃 譯)1992電感計算手冊(北京:機械工業出版社)第1—9頁]

[22]Che B X 2015M.S.Thesis(Changsha:National University of Defense Technology)(in Chinese)[車碧軒2015碩士學位論文(長沙:國防科學技術大學)]

[23]Polzin K A,Choueiri E Y 2006IEEE Trans.Plasma Sci.343