車碧軒 李小康 程謀森 郭大偉 楊雄

(國防科技大學(xué)航天科學(xué)與工程學(xué)院,長沙 410073)

1 引 言

脈沖感應(yīng)推力器(pulsed inductive thruster,PIT)是一種通過脈沖感應(yīng)電磁場電離和加速等離子體從而產(chǎn)生沖量的空間電推進(jìn)裝置,具有工質(zhì)適用范圍廣、比沖效率較高、負(fù)荷大功率、推力變比高等諸多優(yōu)點(diǎn),是未來行星際空間貨物運(yùn)輸或載人深空探測任務(wù)的優(yōu)選推力器之一,近年來引起了廣泛關(guān)注[1?4].PIT尺寸較大,工作電壓達(dá)數(shù)十千伏,需要在高真空環(huán)境下開展實(shí)驗(yàn),實(shí)驗(yàn)系統(tǒng)十分復(fù)雜;其所涉及的等離子體過程具有微秒尺度的強(qiáng)瞬態(tài)特性,諸多等離子體診斷手段并不適用.因此,數(shù)值仿真成為研究PIT中等離子體瞬態(tài)參數(shù)特性、預(yù)測其推進(jìn)性能的重要手段.

PIT的單個脈沖工作過程中存在驅(qū)動電路放電和等離子體感應(yīng)加速兩類基本物理過程,二者通過激勵線圈產(chǎn)生強(qiáng)烈的感應(yīng)耦合作用.Lovberg和Dailey[4?6]最早提出機(jī)電模型,將電路和等離子體依次等效為變壓器的主次級,并將等離子體視作質(zhì)量不斷增長的“雪耙”電流片,建立集總參數(shù)的電流片運(yùn)動及電路方程組求解;Polzin等利用機(jī)電模型研究了推力器性能優(yōu)化的動態(tài)參數(shù)匹配問題[7,8],并引入對Ar等離子體物性參數(shù)的計算[9];Martin等[10]將其適用范圍由平板型線圈拓展至錐形.機(jī)電模型能較好地反映電路與等離子體之間的耦合關(guān)系,但忽略了復(fù)雜的等離子體參數(shù)分布,計算結(jié)果不夠準(zhǔn)確;還存在初始電流片厚度、電流片-激勵線圈互感等諸多關(guān)鍵參數(shù)需要根據(jù)實(shí)驗(yàn)甚至經(jīng)驗(yàn)給出,其計算結(jié)果的正確性難以判斷[11].針對這一缺陷,近年來開始采用二維軸對稱的磁流體力學(xué)(magnetohydrodynamic,MHD)模型計算等離子體感應(yīng)加速過程:Mikellide等[12?14]利用MACH2代碼對NH3工質(zhì)的推力器MK-Va進(jìn)行了數(shù)值仿真,其依據(jù)實(shí)驗(yàn)所得電路電流曲線給出MHD中的時變磁場邊界條件,計算所得推力器單脈沖沖量及放電前1/4個周期內(nèi)的等離子體瞬態(tài)磁場分布與實(shí)驗(yàn)數(shù)據(jù)均能較好地匹配;成玉國等[15,16]采用M-AUSMPW格式計算了正弦時變磁場激勵下Ar等離子體的瞬態(tài)參數(shù)分布,證實(shí)在感應(yīng)加速過程中等離子體形成了“雪耙”型的電流片結(jié)構(gòu).相較于集總參數(shù)的機(jī)電模型,MHD模型側(cè)重于對等離子體流場形態(tài)進(jìn)行計算,考慮了復(fù)雜的各類等離子體物理過程,計算結(jié)果更加可靠.然而,已有各類MHD模型均采用外部給定的瞬時磁場邊界條件驅(qū)動等離子體,沒有考慮等離子體對電路的反向耦合作用,隨著計算時間的推進(jìn),計算結(jié)果會逐漸偏離實(shí)際[14],因此計算結(jié)果的準(zhǔn)確性有待提升.綜上,對于PIT這樣一個涉及強(qiáng)瞬態(tài)等離子體流動過程且電路-等離子體高度耦合的復(fù)雜系統(tǒng),數(shù)值仿真研究還有待進(jìn)一步拓展.

本文結(jié)合兩類數(shù)值仿真模型的優(yōu)缺點(diǎn),在COMSOL平臺中建立了一種耦合外部電路的二維瞬態(tài)軸對稱MHD模型并開展數(shù)值仿真研究.本文第二部分系統(tǒng)介紹了PIT的工作原理以及數(shù)值模擬所采用的物理模型;第三部分對比了模型計算結(jié)果與文獻(xiàn)實(shí)驗(yàn)數(shù)據(jù),驗(yàn)證了計算結(jié)果的有效性;第四部分對計算結(jié)果進(jìn)行了分析與討論,合理解釋了PIT的工作物理圖景,并對電路與等離子體之間的雙向耦合作用進(jìn)行了分析;第五部分為總結(jié).本文建立的數(shù)值仿真模型實(shí)現(xiàn)了對驅(qū)動電路放電過程和等離子體二維流場形態(tài)演化過程的同步耦合求解,計算結(jié)果更加可靠、準(zhǔn)確.模型可用于研究PIT的工作原理,預(yù)測其推進(jìn)性能,為推力器結(jié)構(gòu)設(shè)計、參數(shù)優(yōu)化提供了技術(shù)支持和可行的計算方法.

2 模型簡介

2.1 工作原理

圖1 美國TRW公司MK-1推力器 (a)推力器實(shí)物照片[26];(b)激勵線圈構(gòu)型及線圈面板Fig.1.American TRW Company’s MK-1 thruster:(a)Photograph of the thruster assembly;(b)drive-coil con fi guration and coil-faceplate.

圖2 PIT工作原理示意圖 (a)氣體噴注;(b)感應(yīng)加速Fig.2.Schematic plot of PIT’s working principle:(a)Gas injection;(b)inductive acceleration.

以美國TRW公司1 m直徑的推力器MK-1作為數(shù)值模擬對象,其基本結(jié)構(gòu)如圖1(a)所示,由噴注塔、激勵線圈、電容組C、開關(guān)S等組成.其中,激勵線圈如圖1(b)所示,一般由多組螺旋線形導(dǎo)線按軸對稱方式并聯(lián)排布而成,再封裝于絕緣材料制成的平板型線圈面板中以避免和等離子體直接接觸.推力器的工作原理如下:如圖2(a),首先噴注塔以脈沖方式將數(shù)毫克氣體噴注至線圈面板表面使其達(dá)到較為均勻的徑向分布;緊接著如圖2(b),在氣體噴注完成瞬間,S接通C放電形成迅速上升的脈沖強(qiáng)電流I1,電流通過激勵線圈在空間中激發(fā)徑向發(fā)散型的時變磁場;根據(jù)“法拉第”電磁感應(yīng)定律,這一時變磁場將感生出周向閉合的渦旋電場;渦旋電場擊穿氣體并建立起周向閉合的等離子體電流;等離子體電流(垂直紙面向外)與線圈電流(垂直紙面向內(nèi))相當(dāng)于兩個異號的電流片,其洛倫茲力相互排斥,從而對等離子體實(shí)現(xiàn)軸向加速;等離子體電流片在向前運(yùn)動的過程中不斷電離和帶走下游氣體,最終噴射進(jìn)入真空環(huán)境,完成一個工作脈沖.線圈對等離子體的有效耦合距離約數(shù)厘米,為了約束氣體在面板外徑處設(shè)置一圈與耦合距離等高的圍壩,圍壩所圍成的中空柱形區(qū)域即是等離子體的加速通道[12].

2.2 物理模型

2.2.1 模型總體方案

根據(jù)PIT工作原理及國內(nèi)外已有研究成果,本文建立模型時所考慮的各類物理過程所具備的主要特征以及相應(yīng)的假設(shè)和近似處理總結(jié)如下:1)假設(shè)等離子體電流只存在周向分量,且主要為電子電流,忽略離子電流;2)載流的電子受洛倫茲力作用與離子產(chǎn)生軸向空間電荷分離,離子的加速本質(zhì)上由分離電場實(shí)現(xiàn),但由于電子數(shù)密度較高,其分離尺度相較于電流片尺度可以忽略,電子和離子視作單一流體被加速,滿足等離子體的宏觀電中性條件[5];3)圍壩、線圈面板、噴注塔等表面不存在明顯燒蝕,鞘層現(xiàn)象被證明可以忽略;4)電子和離子幾乎處于熱力學(xué)平衡狀態(tài),等離子體溫度可達(dá)2 eV以上,存在顯著的多級電離和等離子體輻射[12].

本文借助多物理場仿真平臺COMSOL MultiphysicsTM建立和求解計算模型,采用如圖3所示的二維軸對稱計算區(qū)域,假設(shè)等離子體只存在于激勵線圈前側(cè)的矩形加速通道域D1內(nèi)(對應(yīng)圖2(b)虛線內(nèi));激勵線圈則等效為長條形區(qū)域D2,具有垂直于紙面方向的周向電流密度分布jC,θ;噴注塔、圍壩、線圈面板及推力器工作時所處的真空環(huán)境等因其電導(dǎo)率均為0而統(tǒng)一作為真空域D3;為了準(zhǔn)確計算推力器磁場在空間中的分布,D3尺寸遠(yuǎn)大于D1,D2.

圖3 計算模型區(qū)域劃分與邊界Fig.3.Computational domains and boundaries.

圖4 計算模型中所采用的各物理場耦合關(guān)系及模型求解流程Fig.4.Coupling relationship between physical fi elds in the computational model and its solving fl ow chart.

對應(yīng)的模型多物理場耦合關(guān)系及求解流程如圖4:第一步,考慮如圖2所示的驅(qū)動電路,設(shè)激勵線圈端口間的電勢降為VP,建立集總參數(shù)的電路方程組計算線圈總電流I1,其中,VP根據(jù)D2的周向感應(yīng)電場強(qiáng)度Eθ分布及激勵線圈幾何構(gòu)型計算;第二步,根據(jù)激勵線圈幾何構(gòu)型及I1,計算D2中的等效線圈電流密度分布jC,θ,將jC,θ作為外部電流密度施加至D2,采用磁場模塊同時計算D1,D2,D3組成的半圓形區(qū)域磁場分布,再單獨(dú)計算D1的歐姆熱功率密度j2/σ和洛倫茲力密度j×B,以及D2中的周向感應(yīng)電場強(qiáng)度Eθ;第三步,將j2/σ和j×B分別作為體積熱源項和體積力項引入Navier-Stokes(N-S)方程組中的能量方程和動量方程,采用高馬赫數(shù)流動模塊計算D1的溫度T、壓力p等狀態(tài)參數(shù)分布及速度分布u,從而實(shí)現(xiàn)對MHD過程的求解.其中,A○為磁場邊界,B○,C○依次為流場的壁面和出口邊界.

特別地,D1區(qū)域所涉及的等離子體其電離組分不斷變化,等離子體的熱力學(xué)性質(zhì)、輸運(yùn)性質(zhì)等隨等離子體狀態(tài)不同而有很大差異.假設(shè)其處于局部熱力學(xué)平衡狀態(tài)(local thermal equilibrium,LTE),電子和重粒子組分具有相同的溫度T,給定狀態(tài)參數(shù)(p,T)即可惟一確定等離子體的組分及各項物性參數(shù).本文基于LTE模型預(yù)先計算出Ar等離子體在p=1—106Pa,T=298—105K參數(shù)范圍內(nèi)的電離組分及物性參數(shù),將其制作為按(p,T)索引的二維插值數(shù)表,在每一時間步更新.

2.2.2 控制方程組

1)磁場控制方程

磁場控制方程應(yīng)用于圖3所示的所有計算域,其基本形式如下:

式中,磁場通量密度B采用矢量磁勢A(B=?×A)形式求解,從而避免了在MHD數(shù)值計算中因磁場散度清除條件?·B=0限制帶來的數(shù)值誤差[17];σ,μ0表示電導(dǎo)率、真空磁導(dǎo)率;u為速度矢量;u×(?×A)表示由導(dǎo)體宏觀定向運(yùn)動引起的感生電場,該項只在D1中存在;jex表示外部電流密度分布矢量,該項只在D2中存在且只含周向分量,等于線圈電流密度分布jC,θ.

2)流場控制方程組

采用可壓縮黏性流體控制方程組,只在D1域求解,具體形式如下:

連續(xù)方程

動量方程

式中,ρ為密度,I為單位對角張量,τ為黏性應(yīng)力張量,B為磁感應(yīng)強(qiáng)度j×B為洛倫茲力體積力;能量方程[18]

式中,cp表示等離子體定壓比熱,k為等離子體熱導(dǎo)率,S為流體微元應(yīng)變張量;等號右邊五項依次為熱傳導(dǎo)項、黏性耗散項、壓力耗散項、等離子體歐姆熱源項、等離子體輻射損失項qrad.

3)狀態(tài)方程

在COMSOL中采用如下形式的氣體狀態(tài)方程:

其中,在常溫下Rs一般表示理想氣體常數(shù);對于本文涉及的高溫電離氣體,因其存在多級電離組分,Rs不再是常數(shù),而是等離子體狀態(tài)參數(shù)(p,T)的函數(shù)[19],Rs計算方法將在后文闡述.

4)電路控制方程組

考慮圖2所示的驅(qū)動電路建立集總參數(shù)的電路方程組:

其中,I1為驅(qū)動電路總電流,IP為等離子體總電流,R0為電路寄生電阻,L0為電路寄生電感,C為電容組總?cè)葜?VC為電容電壓,VP為激勵線圈端口間的電勢降.

5)電路-等離子體雙向耦合參數(shù)計算方程

驅(qū)動電路與等離子體通過激勵線圈相互耦合,耦合參數(shù)jC,θ,VP分別表示電路對的等離子體的激勵與反饋?zhàn)饔?與激勵線圈的幾何構(gòu)型緊密相關(guān).對MK-1的激勵線圈,其單支螺旋線可表示為如圖5所示的等速螺旋線,導(dǎo)線由線圈外側(cè)r2處回繞至r1處,再由絕緣面板的背側(cè)連出.其兩端IO通過同軸電纜與電容及開關(guān)相連,螺旋線線形可表示為

圖5 激勵線圈中的單支螺旋線導(dǎo)線幾何形狀Fig.5.The geometry of one individual spiral conductor in drive-coil.

激勵線圈域D2中任意位置處的等效周向電流密度jC,θ可表示為

式中,δc為D2區(qū)域的高度,等于激勵線圈的導(dǎo)體直徑.

激勵線圈端口電勢降VP則等于周向電場強(qiáng)度Eθ沿螺旋線導(dǎo)線的積分,

對上式做積分變換可得適用于二維軸對稱的表達(dá)式:

其中,Eθ可根據(jù)磁場分布情況由式Eθ=(??A/?t)θ計算.

2.2.3 初邊值條件

1)初值條件

電路控制方程組初值條件:初始電路電流I1=0,初始電容電壓V=V0.

D1,D2,D3域磁場初值條件:A=0.

D1域流場初值條件:本文不考慮中性氣體的擊穿過程,假設(shè)在初始時刻加速通道內(nèi)的工質(zhì)已被輕度電離,其等離子體溫度T=0.5 eV;根據(jù)室溫下給定的中性氣體密度分布,計算其T升高至0.5 eV時的p分布,以之作為p分布的初值條件.

2)邊值條件

對于磁場邊界條件,由于三部分計算子域磁場是連續(xù)的,故只需給出圖示邊界A○的磁場邊界條件,這里采用磁絕緣邊界條件;對于流場邊界條件,線圈面板、噴注塔及圍壩組成的壁面B○采用速度無滑移、絕熱壁面邊界條件;開口邊界C○采用遠(yuǎn)場邊界條件.

2.2.4 等離子體物性參數(shù)

1)電離組分

計算等離子物性參數(shù),首先需要獲得其電離組分,針對Ar工質(zhì)在LTE假設(shè)下考慮6級平衡電離復(fù)合反應(yīng),聯(lián)立沙哈平衡方程組(12)、電荷守恒方程(13)及混合組分氣體狀態(tài)方程(14)求解[9]:

式中,ne表示電子數(shù)密度;ns表示s級(s=0—6)電離離子數(shù)密度,特別地,當(dāng)s=0時n0表示中性原子數(shù)密度;me為電子質(zhì)量;κB為玻爾茲曼常數(shù);h為普朗克常數(shù);εsl為s級電離離子的第l能級,gsl為相應(yīng)的統(tǒng)計權(quán)重.

2)熱力學(xué)參數(shù)

根據(jù)電離組分的數(shù)密度計算等離子體總ρ,得到其Rs(p,T):

cp通過比焓ht計算[20]:

對于單原子分子等離子體,ht由平動能貢獻(xiàn)htran,重粒子內(nèi)部電子的激發(fā)能貢獻(xiàn)hex及電離離子的電離能貢獻(xiàn)hi組成:

各部分貢獻(xiàn)的計算式依次為:

式中,Es表示s級電離離子的電離能,zs表示s級電離離子的內(nèi)部配分函數(shù).

3)輸運(yùn)系數(shù)

在等離子體的輸運(yùn)性質(zhì)中,電導(dǎo)率是極其重要的一種,本文所涉及的等離子體可能處于完全電離或弱電離之間的任意電離狀態(tài),需要同時考慮長程碰撞和短程碰撞的作用.忽略磁場對電導(dǎo)率的影響,假設(shè)電導(dǎo)率為各向同性的標(biāo)量σ,采用Kantrowitz[21]提出的電導(dǎo)率并聯(lián)模型計算復(fù)合電導(dǎo)率:

其中,Spitzer[22]給出了完全電離等離子體的電導(dǎo)率σc計算公式:

對于弱電離等離子體電導(dǎo)率σw,采用經(jīng)典碰撞理論計算:

式中e為電子電荷量,ves為能量加權(quán)平均的電子和s級電離離子之間進(jìn)行動量交換的碰撞頻率[23,24].

對于等離子體輻射,其光學(xué)厚部分作為等離子體熱傳導(dǎo)的一種增強(qiáng)機(jī)制,依據(jù)文獻(xiàn)[25—27]計算了等離子體復(fù)合熱導(dǎo)率;其光學(xué)薄部分作為體積熱源項qrad引入能量方程.qrad及黏性系數(shù)μ依據(jù)文獻(xiàn)[25]計算.

3 數(shù)值模擬與實(shí)驗(yàn)數(shù)據(jù)對比

數(shù)值仿真所采用的初始ρ分布如圖6,取自文獻(xiàn)[26]給出的MK-1推力器實(shí)測值,對應(yīng)的脈沖氣體質(zhì)量mbit=15 mg,電容初始電壓V0=20 kV.為驗(yàn)證計算結(jié)果的有效性,本文對比了不同時刻下徑向磁場強(qiáng)度Br沿內(nèi)外徑中線處的軸向分布(圖7)、軸向洛倫茲力密度(j×B)z的二維分布(圖8)、推力器推力-時間曲線(圖9)以及V0=20—24 kV下的比沖Isp、效率η等性能參數(shù)(圖10)的數(shù)值仿真結(jié)果與實(shí)驗(yàn)測量數(shù)據(jù).

圖6 數(shù)值仿真所采用的加速通道初始?xì)怏wρ/(kg/m3)分布,數(shù)據(jù)取自文獻(xiàn)[26]Fig.6.Initial ρ/(kg/m3)distribution in the acceleration channel used in simulation,data from reference[26].

3.1 徑向磁場強(qiáng)度Br

通常情況下,Br將隨軸向距離z的增大而按指數(shù)衰減,但在如圖7所示的2μs時刻z=0—0.02 m區(qū)間卻形成了磁場強(qiáng)度達(dá)0.35 T的均勻磁場區(qū)域.這是由于等離子體電流片的存在壓縮了其與激勵線圈之間的磁場.伴隨電流片逐漸遠(yuǎn)離線圈,該均勻磁場區(qū)域的范圍逐漸擴(kuò)大,磁場強(qiáng)度逐漸減小.數(shù)值仿真結(jié)果重現(xiàn)了這一磁場演化趨勢,同時在磁場強(qiáng)度的大小及軸向分布上也與實(shí)驗(yàn)數(shù)據(jù)基本一致.

圖7 t=2,4,6μs時刻加速通道內(nèi)外徑中線處Br計算結(jié)果與文獻(xiàn)[28]實(shí)驗(yàn)數(shù)據(jù)Fig.7.Calculated and experiment measured[28]Bron the middle line of the acceleration channel at t=2,4,6μs.

3.2 軸向洛倫茲力密度(j×B)z

(j×B)z分布呈現(xiàn)出與Br不同的演化趨勢:其峰值區(qū)域在t=2μs首先出現(xiàn)在z=0附近;但在t=3μs,線圈附近的(j×B)z開始減小,其峰值區(qū)域向前推進(jìn),洛倫茲力的作用似乎穿透了均勻磁場區(qū)域而主要作用在電流片上;由于初始?xì)怏wρ分布的不均勻,靠近內(nèi)徑處的(j×B)z峰值區(qū)域運(yùn)動得比靠近外徑處更快.圖8(b)所示的計算結(jié)果基本呈現(xiàn)出了與實(shí)驗(yàn)數(shù)據(jù)一致的演化規(guī)律,其(j×B)z分布在t=2μs時刻與實(shí)驗(yàn)數(shù)據(jù)符合較好,在t=3μs存在一定差異.

圖8 t=2,3μs時刻(j×B)z/(N/m3)分布文獻(xiàn)[26]實(shí)驗(yàn)數(shù)據(jù)(a)與本文計算結(jié)果(b)Fig.8.Calculated(a)and experiment measured[26](b)data of(j×B)z/(N/m3)at t=2,3μs.

3.3 推力-時間曲線

將工質(zhì)所受洛倫茲力在加速通道內(nèi)做體積積分即可得到推力器的推力-時間曲線.由圖9可見,推力器的推力主要產(chǎn)生在0—9μs時刻,在9—10μs期間短暫為負(fù),10μs之后再次為正并逐漸衰減.數(shù)值計算結(jié)果在趨勢及大小上均與實(shí)驗(yàn)數(shù)據(jù)一致,表明數(shù)值仿真正確反映了推力器在不同時刻的整體工作狀態(tài).

圖9 推力器-時間曲線計算結(jié)果與文獻(xiàn)[28]實(shí)驗(yàn)數(shù)據(jù)對比Fig.9. Calculated and experiment measured[28]thrust vs.time.

3.4 不同放電電壓下的推力器性能

圖10給出了MK-1推力器在V0=18—24 kV下的推進(jìn)性能實(shí)驗(yàn)數(shù)據(jù)及計算結(jié)果對比,其中,單脈沖沖量I通過推力在0—20μs計算時間內(nèi)的積分得到,比沖Isp=I/mbit,效率η為工質(zhì)總動能與電容組儲能之比:

作為對比,圖中同時給出了Lovberg等[26]采用一維機(jī)電模型的計算結(jié)果.圖10表明,本文計算所得推力器性能略低于與實(shí)驗(yàn)測量數(shù)據(jù),但相較于“雪耙”模型的計算結(jié)果更加接近實(shí)驗(yàn)結(jié)果,推力器性能隨V0的變化趨勢則與實(shí)驗(yàn)數(shù)據(jù)基本一致.

圖10 不同放電電壓下的推力器比沖Isp、效率η性能計算結(jié)果、實(shí)驗(yàn)測試數(shù)據(jù)及一維“雪耙”機(jī)電模型計算結(jié)果對比Fig.10.Comparison of the thruster’s speci fi c impulse Isp and efficiency η:numerical results in this paper;experiment measured data in[26];numerical result from 1D circuit model in[8].

綜上,本文所采用的計算模型能較好地反映推力器工作過程中的各類物理過程,數(shù)值仿真成功復(fù)現(xiàn)了PIT的工作物理圖景,計算結(jié)果在等離子體瞬態(tài)參數(shù)分布、推力器工作狀態(tài)、推力器推進(jìn)性能三個方面均與實(shí)驗(yàn)取得了較好的一致.

4 計算結(jié)果分析與討論

4.1 推力器系統(tǒng)的工作物理圖景

圖11 計算所得等離子體ρ分布(a)與jθ分布(b)Fig.11.Calculated plasma ρ (a)and jθ (b)distribution.

為了對推力器工作物理圖景建立更深刻的認(rèn)識,圖11給出了不同時刻ρ和jθ的二維分布.如圖所示,t=2μs時激勵線圈附近產(chǎn)生了與線圈電流反向的(負(fù)號)等離子體電流jθ,jθ受洛倫茲力作用被軸向加速,帶動工質(zhì)壓縮成片狀,同時在其后方留下了近乎真空的低ρ;由于這一區(qū)域ρ極低,載流電子較少,其σ也較低,因此電流片開始逐漸脫離線圈表面,線圈磁場也得以穿過這一區(qū)域,持續(xù)地對高ρ,高jθ的電流片區(qū)域進(jìn)行加速;洛倫茲力在前期將等離子體壓縮成致密的片狀,伴隨電流片不斷遠(yuǎn)離,激勵線圈磁場逐漸衰減,6μs時刻洛倫茲力不再足以維持其致密的片狀結(jié)構(gòu),電流片開始破裂,等離子體則向真空中自由膨脹.

4.2 電路-等離子體的雙向耦合作用

借助本文所建立的模型,可以實(shí)現(xiàn)對推力器工作過程中電路狀態(tài)與等離子體狀態(tài)的同步分析.圖12在同一坐標(biāo)系內(nèi)給出了驅(qū)動電路電流I1與等離子體總電流的相反數(shù)?IP隨時間變化的曲線,其中IP通過jθ在D1域的積分得到.由圖可見,在放電前期I1與?IP幾乎完全重合,隨后逐漸分離,但在變化趨勢上?IP始終跟隨I1.對應(yīng)圖9所給出的推力-時間曲線可以得出,等離子體的感應(yīng)加速主要集中在I1的前1/2個周期內(nèi)I1與IP異號的階段;在9—10μs期間,I1與IP同號,推力為負(fù);在10μs之后,I1與IP再次異號推力為正并逐漸衰減.這些現(xiàn)象表明,電路電流實(shí)際上決定了等離子體電流及等離子體的加速進(jìn)程.另一方面,為了分析等離子體對電路電流的影響,本文通過在計算域中移除D1域計算了激勵線圈在真空下放電的空載電流曲線I?1.在圖12中對比I1與I?1可以見,耦合等離子體之后,I1相較于I?1振蕩周期明顯縮短,相當(dāng)于其等效電感減小;圖13則給出了對應(yīng)的空、負(fù)載放電狀態(tài)下的激勵線圈電勢降VP和V?P.由圖可見,耦合等離子體后VP顯著降低,且在t=0—2μs期間保持基本穩(wěn)定,這一特殊現(xiàn)象說明VP受到了等離子體瞬時流場形態(tài)的影響.以上分析表明,在PIT工作過程中其驅(qū)動電路與等離子體相互影響,存在強(qiáng)烈的雙向耦合作用.

圖12 計算所得等離子體總電流IP與電路總電流I1隨時間的變化Fig.12.Calculated total plasma current IPand circuit current I1vs.time.

圖13 計算所得空、負(fù)載放電狀態(tài)下的VPFig.13. Calculated drive-coil VPwhen plasma is loaded or unloaded.

為了對這一雙向耦合作用進(jìn)行定量化分析,首先將驅(qū)動電路和等離子體等效為如圖14的變壓器主次級.其中,LC為激勵線圈自感,LP為等離子體電流環(huán)自感,RP為等離子體等效總電阻;激勵線圈與等離子體之間的互感M體現(xiàn)了二者之間的感應(yīng)耦合強(qiáng)度;顯然,LP,RP和M均將隨等離子體狀態(tài)參數(shù)及流場形態(tài)的演化而改變[6].

圖14 PIT電路等離子體系統(tǒng)的變壓器等效電路圖Fig.14. Transformer equivalent circuit scheme for PIT’s drive circuit and plasma load.

采用前期研究提出的方法[11]計算所得RP和M結(jié)果如圖15和圖16,下面對照已給出的等離子體jθ瞬態(tài)分布對電路-等離子體之間的雙向耦合作用進(jìn)行定量分析:在放電初始時刻,激勵線圈表面附近建立起了環(huán)形的等離子體電流片,jθ因趨膚效應(yīng)作用聚集在激勵線圈附近,與線圈電流幾乎重合,因而M≈LC=0.75μH;此后伴隨電流片被加速遠(yuǎn)離線圈,M逐漸減小,其減小的速度也隨電流片運(yùn)動速度的增大而增大;t=10μs時刻I1過零點(diǎn)附近,因jθ再次聚集于激勵線圈表面,M又開始迅速增大,隨后再次減小.Rp則體現(xiàn)了等離子體的總歐姆耗散,放電前;放電后0—6μs期間,Rp因氣體被持續(xù)電離而不斷減小;6—8μs伴隨感應(yīng)耦合作用M的減弱,Rp開始增大;8—16μs再次下降.

圖15 激勵線圈與等離子體之間的互感MFig.15.Mutual inductance M between plasma load and drive-coil.

圖16 等離子體等效總電阻RPFig.16.Effective plasma total resistance RP.

綜上所述,耦合等離子體負(fù)載將驅(qū)動電路總體電阻增大,電感減小;等離子體與驅(qū)動電路的耦合作用強(qiáng)度隨等離子體瞬態(tài)參數(shù)分布及電路放電狀態(tài)的變化而變化;當(dāng)?shù)入x子體電流在激勵線圈附近集聚時,M將增大,反之減小.

5 結(jié) 論

建立了一種耦合外部電路的磁流體力學(xué)模型,實(shí)現(xiàn)了對加速通道內(nèi)等離子體二維流場結(jié)構(gòu)演化過程及驅(qū)動電路放電過程的同步耦合求解.數(shù)值仿真成功復(fù)現(xiàn)了PIT的工作物理圖景,計算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)基本一致,借助這一新模型,實(shí)現(xiàn)了對電路-等離子體雙向耦合作用的定量分析.模型可用于研究PIT的工作原理,預(yù)測其推進(jìn)性能,為推力器結(jié)構(gòu)設(shè)計、參數(shù)優(yōu)化提供技術(shù)支持和可行的計算方法.

[1]Polzin K A 2011J.Prop.Power27 3

[2]Martin A K,Dominguez A,Eskridge R H 201534th International Electric Propulsion ConferenceHyogo-Kobe,Japan,July 4–10,2015 p50

[3]Russell D,Poylio J H,Goldstein W 2004Space Conference and ExhibitSan Diego,America,September 28–30,2004 p6054

[4]Dailey C L,Loveberg R H 1987Pulsed Inductive Thruster Component TechnologyAFAL TR 07 012

[5]Dailey C L,Loveberg R H 1989AIAA/ASME/SAE/ASEE 25th Joint Propulsion ConferenceMonterey,America,July 10–12,1989 p2266

[6]Dailey C L,Lovberg R H 1993The PIT MkV Pulsed Inductive ThrusterNASA CR 19 1155

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

[8]Polzin K A 2006Ph.D.Dissertation(Princeton:Princeton University)

[9]Polzin K A,Sankaran K,Ritchie A G,Reneau J P 2013J.Phys.D:Appl.Phys.46 475201

[10]Martin A K 2016J.Phys.D:Appl.Phys.49 025201

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

[12]Mikellides P G,Neilly C 2007J.Prop.Power23 51

[13]Mikellides P G,Ratnayake N 2007J.Prop.Power23 854

[14]Mikellides P G,Villarreal J K 2007J.Appl.Phys.102 103301

[15]Cheng Y G,Xia G Q 2017Acta Phys.Sin.66 075204(in Chinese)[成玉國,夏廣慶 2017物理學(xué)報 66 075204]

[16]Cheng Y G 2015Ph.D.Dissertation(Changsha:National University of Defense Technology)(in Chinese)[成玉國2015博士學(xué)位論文(長沙:國防科技大學(xué))]

[17]ComsolMultiphysicsUsersGuide,Onlinevailable http://www comsol com/plasma-module/[2017-5-11]

[18]Li M,Liu H,Ning Z X 2015IEEE Trans.Plasma Sci.43 12

[19]John D A(translated by Yang Y)2011Hypersonic and High-Temperature Gas Dynamics(2nd Ed.)(Beijing:Aviation Industry Press)pp421–422(in Chinese)[小約翰D A著(楊永譯)2011高超聲速和高溫氣體動力學(xué)(第二版)(北京:航空工業(yè)出版社)第421—422頁]

[20]Cheng X 2009Thermal Plasma Heat Transfer and Flow(Bejiing:Science Press)pp50–55(in Chinese)[陳熙 2009熱等離子體傳熱與流動(北京:科學(xué)出版社)第50—55頁]

[21]Deb P,Agarwal R K 2001AIAA Aerospace Sciemces Meeting&ExhibitReno,America 2001,p794

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

[23]Ahangar M,Ebrahimi R,Shams M 2014Acta Astronaut.103 129

[24]Heiermann J 2002Ph.D.Dissertation(Stuttgart:Universitat Stuttgart)

[25]Sankaran K 2005Ph.D.Dissertation(Princeton:Princeton University)

[26]Glumb R J,Krier H 1986AIAA J.24 1331

[27]Lovberg R H,Dailey C L 1982AIAA/JSASS/DGLR 16th International Electric Propulsion ConferenceNew Orleans,America,November 17–19,1982 p1921

[28]Lovberg R H,Dailey C L 1982AIAA J.20 971