外磁場作用下的磁等離子體動力學過程仿真

羅楊 陳茂林 蘇冬冬 許諾 王忠晶 韓志聰 趙豪

1) (北京機械設備研究所,北京 100854)

2) (西北工業大學,燃燒、熱結構與內流場重點實驗室,西安 710072)

磁等離子體動力學推力器是空間高功率電推進裝置的典型代表,磁等離子體動力學過程是其核心工作機制.為深入理解外磁場對其工作特性的影響,本文采用粒子云(particle in cell,PIC)方法結合基于自相似準則的縮比模型,進行外加磁場作用下磁等離子體動力學推力器工作過程的建模仿真,通過與實驗結果對比驗證模型和方法的可靠性,并重點分析推力器點火啟動過程的等離子特性參數分布,以及外磁場和陰極電流對推力器工作性能的影響.研究結果表明:陰陽極放電電弧構建是推力器啟動和高效工作的關鍵步驟;外磁場強度較低工況不利于構建穩定放電電弧,等離子體束流集中于軸線附近,推力主要產生機制是自身場加速;外磁場強度較高時,陰陽極放電電弧穩定,推力產生主要機制是渦旋加速,推力、比沖隨外磁場強度線性增大;推力器效率隨陰極電流和外磁場強度增大而增大;放電電壓隨陰極電流增大而增大,但隨外磁場強度的增大表現出先減小后增大的趨勢.

1 引言

磁等離子體動力學推力器(magnetoplasmadynamic thruster,MPDT)是通過洛倫茲力產生推力的航天推進裝置,也稱洛倫茲力加速器,可實現兆瓦級功率和百牛級推力性能,是高功率、大推力電推力器的重要發展方向.MPDT 在美國、前蘇聯/俄羅斯、日本、德國、意大利、中國均得到了廣泛關注和重點研究[1,2],實驗室內實現了數千瓦至數兆瓦功率的MPDT 測試,最大測試推力達100 N 以上.

研究人員對MPDT 進行了大量的理論分析、仿真模擬和實驗測試等工作.在仿真模擬研究方面早期主要采用流體模型:Brushlinskii 和Morozov[3,4]最早于20 世紀60 年代采用流體模型對大功率MPDT 開展數值模擬,獲得了MPDT 軸向密度分布與電流的關系;Kimura 等[5]提出了單溫度的二維計算模型,并使用該模型對MPDT 的工作過程進行模擬;Niewood[6]和Sheppard 等[7]提出了雙溫度軸對稱數值模型來研究計算推力器內部的流場分布,考慮了非平衡電離過程和中性氣體的影響;Caldo 等[8]采用雙溫度模型來研究了反常運輸對MPDT 流動的影響.Rudolph[9]在MPDT 流體模擬過程中采用多重網格和多重時間步長法,解決了時間依賴性問題.Heiermann[10]、趙博強等[11]、Sary等[12]、Tahsini[13]采用流體模型對MPDT 內部的電離過程和多種輸運過程進行了詳細建模仿真,實現了對MPDT 推力器的整體性能計算.

流體模型適用于MPDT 的高密度等離子體特性,具有計算適用性強、計算約束條件少、計算速度快的特點.但是流體作為基于宏觀的平均化方法,流體模型無法捕捉到等離子體的一些瞬時現象.粒子云(particle in cell,PIC)[14]等基于動力學過程的等離子體粒子模擬方法是捕捉MPDT工作過程中等離子體演化過程和瞬時現象的優選方法.由于MPDT 高等離子體密度特性對PIC 仿真時間步長和空間步長的限制,使得其對計算資源和計算機時需求巨大,故早期粒子模擬方法在MPDT仿真中應用較少.Tang 等[15-17]開展了PIC方法在MPDT 模擬中的應用探索,采用改變重粒子質量的加速計算方法減少了計算機時,實現了低功率下MPDT 工作過程模擬;Li 等[18]采用PIC方法開展了無外加磁場作用下的MPDT 工作過程模擬.本文采用PIC 方法結合基于自相似準則的縮比模型進行外加磁場作用下的MPDT 工作過程模擬,解決MPDT 極高等離子體密度對模擬時間步長、空間步長的限制,在減小計算量的同時,進行真實工況下MPDT 工作過程模擬和性能影響因素分析.

2 模型與方法

2.1 MPDT 結構與工作原理

根據磁場產生方式的不同,MPDT 分為兩大類:自身場磁等離子體動力學推力器(self-field MPDT,SF-MPDT)和附加場磁等離子體動力學推力器(applied-field MPDT,AF-MPDT).本文主要開展AF-MPDT 工作過程的仿真研究,其主要結構包括陰極、陽極和電磁線圈,如圖1 所示.

圖1 AF-MPDT 結構和工作原理示意圖Fig.1.Schematic of AF-MPDT structure and operating principle.

AF-MPDT 主要有3 種加速機制產生推力:氣動加速、自身場加速和附加場加速,附加場加速又包括霍爾加速和渦旋加速.

1)氣動加速.氣動加速分量是通過電弧歐姆加熱將電能轉換為熱能來提高推進劑的總焓,進一步噴管將熱能轉換為動能而產生推力.當 MPDT在高功率工況工作時,電磁加速占主導,氣動加速分量可忽略.

2)自身場加速.圖2(a)為自身場加速原理.MPDT 工作時,電流軸向分量會產生一個周向自感磁場Bθ,該磁場又與電流相互作用產生洛倫茲力:其中徑向電流jr與周向自感磁場Bθ,產生一個軸向的力 (jr×Bθ),該作用力直接對推力產生貢獻;而軸向的電流jz與自感磁場Bθ相互作用產生徑向的洛倫茲力 (jz ×Bθ),該作用力會導致軸線附近的等離子體密度增大,并產生一個不平衡的壓力,間接對推力產生貢獻.由于感應磁場Bθ與電流線性相關,自身場加速產生的推力與電流表現為平方關系.自身場加速機制在SF-MPDT 中貢獻較大,而在AF-MPDT 中貢獻較小.

圖2 MPDT 加速機理 (a)自身場加速;(b)霍爾加速;(c)渦旋加速Fig.2.MPDT acceleration mechanism:(a) Self-field acceleration;(b) hall acceleration;(c) vortex acceleration.

3)霍爾加速.圖2(b)為霍爾加速原理.在AFMPDT 中,外部電磁線圈會在陰陽極間產生軸向磁場Bz和徑向磁場Br.根據歐姆定律,軸向磁場Bz和徑向電場作用下,等離子體會感生出周向電流jθ,與自身場機制相似,jθ會與附加磁場產生箍縮力 (jθ×Bz)和軸向力 (jθ×Br).

4)渦旋加速.圖2(c)為渦旋加速原理.在AFMPDT 中,陰陽極間的軸向和徑向電流分量,也會與外部電磁線圈產生的徑向磁場Br和軸向磁場Bz相互作用,產生周向電磁力 (jr×Bz) 和 (jz ×Br),使等離子體產生角向渦旋,宏觀上表現為自旋的等離子體徑向膨脹,在物理/磁噴管的作用下再轉化為軸向作用力,產生推力.渦旋加速產生的推力與電流和外加磁場的乘積IB成正比,Fradkin[19]的實驗也驗證了該推力關系.

2.2 物理模型

為了便于仿真結果與實驗結果對比分析,采用Myers 團隊設計的AF-MPDT 推力器[20]作為計算模型,其幾何結構如圖3 所示,圖中標注單位為mm,虛線框內為仿真模擬區域.

圖3 MPDT 幾何結構與計算區域Fig.3.MPDT geometry and simulation area.

2.3 仿真方法

2.3.1 PIC 方法

PIC 方法是等離子體粒子模擬的常用方法,其核心過程包括電磁場求解、帶電粒子在電磁場作用下運動、以及粒子碰撞三部分[21].

MPDT 工作過程中的電磁場分布及演化可以通過麥克斯韋方程描述,由于MPDT 工作過程中電磁場的時變量較小,可以忽略,故可將麥克斯韋方程簡化為

式中,B=?×A為磁感應強度,E=-?φ為電場強度,μ0為真空磁導率,J為電流,ρ為電荷密度,ε0為真空介電常數,A為矢量磁位,φ為電勢.進一步將(1)式改寫為泊松方程形式,便于數值求解:

帶電粒子的運動通過牛頓-洛倫茲方程進行描述,如下:

式中,m為粒子質量,v為粒子速度,e為單位電荷,s為粒子位置,t為時間.

PIC 方法處理粒子間的碰撞一般使用蒙特卡羅碰撞(Monte Carlo collision,MCC)方法[22].在MPDT 模擬過程中,考慮如下碰撞類型:電子與中性分子間的彈性碰撞、激發碰撞、電離碰撞,離子與中性分子間的彈性碰撞和電荷交換(charge exchange,CEX)碰撞,電子-電子、電子-離子、離子-離子間的庫侖碰撞,以及電子的Bohm 碰撞.

2.3.2 邊界條件設置

仿真中,主要包含陰極、陽極、絕緣壁面、對稱軸和自由邊界5 種邊界類型.陰極邊界設置為Neumann 邊界,固定法向場強E⊥=σ/ε0,切向場強為E//=0,其中σ為陰極表面電荷密度;固定陰極的軸向電流密度為,徑向電流密度和周向電流密度均為0,其中,I為陰極電流,Rcathode為陰極半徑;陰極作為粒子吸收邊界,粒子碰撞陰極后注銷.考慮離子碰撞陰極壁面的二次電子發射過程,對于鎢陰極,當氬離子的入射能量介于10—100 eV 之間時二次電子發射系數為0.095 左右.

陽極設置為零電勢的Drichlet 邊界和粒子吸收邊界,粒子碰撞陽極后注銷.陰極陽極間的絕緣壁面為Neumann 邊界,可認為壁面表面的場強是由累積電荷得到,且垂直于壁面,即法向場強E⊥=σ/ε0,切向場強為E//=0.對于離子,絕緣壁面為吸收邊界,離子碰撞絕緣壁面后注銷;對于電子,分別考慮在壁面沉積、反射、激發一個二次電子、以及激發兩個二次電子的情況,4 種情況的概率分別為

式中,ξ為電子能量,Pdep(ξ)為電子沉積壁面的概率,Pref(ξ)為電子發射反射的概率,P1se(ξ)為電子撞擊出一個二次電子的概率,P2se(ξ)為電子撞擊出兩個二次電子的概率.

自由邊界主要是粒子的出口,粒子運動越過自由邊界后注銷;計算模型是二維軸對稱模型,中心軸線為對稱邊界,粒子運動出對稱邊界時進行鏡面反射處理;電磁場方程均使用Neumann 邊界為自由邊界,可認為電場、磁場在自由邊界處的法相分量為零.

2.4 基于自相似準則的縮比模型

PIC 仿真過程中,空間步長需要小于德拜長度,時間步長需要小于等離子體頻率,并同時滿足CFL 條件.對于MPDT,電子數密度一般大于1020m—3,對應德拜長度一般小于10—6m,單個自由度的仿真網格數一般需要數萬至數十萬,等離子體頻率約1012Hz,對應時間步長約10—12s.要同時滿足空間步長和時間步長需求,也會給計算量帶來極大挑戰,因此需要一些加速方法在保證仿真準確性的同時,減少計算量.

目前對全粒子PIC 模擬有兩大類加速方法,一種是直接改變介電常數以及重粒子質量來加速計算的方法[23],另一種是Khayms[24]和Taccogna等[25]在霍爾推力器仿真中使用基于自相似準則的縮比模型方法,其基本特征是保證推力器性能參數成比例縮放.

鑒于MPDT 和霍爾推力器工作過程和物理特征相似,本文采用基于自相似準則的縮比模型方法,結合PIC 粒子模擬方法,進行MPDT 仿真模型構建和加速計算.為了獲得準確的推力器性能參數,縮比模型中,相似準則選取如下:1)推力器的比沖Isp不變;2)磁場對電子的約束效果不變;3)電子平均自由程與推力器特征尺寸比值不變.

縮比模型中,將推力器特征尺寸縮小至1/τ,則物理量電流I、質量流量、數密度ne,ni、電場E、磁場B、推力F、比沖Isp和效率η的縮比系數分別為1/τ,1/τ,1/τ,τ,τ,τ,1,1 和1/τ.其中自變量為計算輸入時需要調整的量,因變量表示根據調整后的輸入量而發生改變的物理量.縮比模型將推力器尺寸縮小為原來的1/τ后,電子數密度同樣也擴大為原來的τ倍,使得德拜長度縮小1/τ1/2,可以認為縮比系數τ對單方向的網格數量的縮小量為1/τ1/2,對于軸對稱模型,網格數量縮小1/τ.本文仿真中τ取為8 × 106.

3 仿真結果與分析

3.1 中性氣體流場分布

中性氣體流場分布一般采用直接蒙特卡羅模擬(direct simulation Monte Carlo,DSMC)方法求解,將等離子體放電和輸運過程模擬的PIC 方法與中性氣體流場模擬的DSMC 方法耦合,會大幅增大計算需求.本文將中性氣體流場與等離子體流場解耦處理,即先單獨采用DSMC 方法模擬MPDT 中的中性氣體流場,然后采用電推進仿真中常用的PIC/MCC 方法處理帶電粒子與中性氣體原子間的相互作用,即將具備一定分布的氣體原子作為背景,參與帶電粒子的碰撞處理,但其分布不隨模擬過程變化,這樣可在保證較好的計算精度的同時,大幅縮減計算機時.

針對Ar 工質,流量為=0.1 g/s 的情況下的氣體初始流場,Ar 原子數密度分布和流線如圖4所示,流場內中性粒子數密度量級約在1020—1021m—3量級.

圖4 Ar 原子數密度分布和流線Fig.4.The number density and streamlines of Ar.

3.2 等離子體密度分布演化

通過推力器內部等離子體特性參數的演化,可以分析推力器內部放電和帶電粒子的輸運過程.針對外磁場B=0.03 T,陰極電流I=1000 A 工作參數下的MPDT 的工作過程進行仿真,在縮比模型中,計算時間步長為1 × 10—18s,對應實際情況的時間步長8 × 10—12s.

3.2.1 電子數密度

MPDT 的放電由陰極發射的電子電流主導,故首先分析放電過程中的電子數密度分布隨時間的演化過程,不同時間步的電子數密度分布仿真結果如圖5 所示.

圖5 電子數密度分布Fig.5.Distribution of electron number density.

MPDT 點火啟動過程中,電子數密度分布隨時間的變化可初步分為3 個階段:放電初期,電子從陰極發射后,受外加磁場約束,沿磁力線螺旋運動,呈擴張型分布;放電中期,電子在陰極尖端下游堆積,逐漸形成一個高密度分布區域,并向陽極方向延伸,構成一條連通陰陽極的高密度電子分布帶,形成放電電弧;放電末期,放電電弧逐漸穩定,等離子體區域進一步擴張,并在陰陽極下游形成完整的高密度的穩態分布,推力器進入穩定工作狀態.推力器在達穩態工作狀態后,在推力器外部,電子束流分為兩個部分,一個集中在陰極尖端的軸線處向下游擴散,另一部分電子隨磁力線向下游擴散,這和實驗觀測到的現象一致[1].穩態時,MPDT陽極尖端的電子數密度達1021m—3,陰陽極下游的電子數密度達1020m—3.

3.2.2 離子數密度

陰極發射電子在電磁場共同作用下運動的同時,會碰撞中性氣體分子發生電離,產生新的電子和離子.由于離子完全由電離過程產生,離子數密度分布隨時間的演化過程更真實地反映MPDT 的電離過程,不同時間步的離子數密度分布仿真結果如圖6 所示.

由圖6 可知,在放電初期,離子數密度同樣由小增大,并隨著放電進行,其分布也逐漸穩定.在放電初期,離子數密度的分布與電子分布相似,其原因是離子由電子電離碰撞產生的,在高電子數密度處離子數密度高;隨著放電進行,離子數密度的分布逐漸偏離電子分布,呈現出高密度區域逐漸由陰極尖端附近向整個下游區域擴張并向軸線附近聚焦的趨勢,離子分布偏離電子分布的原因是離子質量大,受磁場的約束遠小于電子,故運動軌跡和最終分布也與電子不同.穩態工作時,離子數密度最高達到1021m—3,主要仿真區域離子數密度約1020m—3,與電子數密度相當.

圖6 離子數密度分布Fig.6.Distribution of ion number density.

3.3 外加磁場對MPD 性能影響及其與實驗數據的對比分析

外加磁場是影響MPDT 推力器性能的重要參數,在陰極電流I=1000 A,流量為=0.1 g/s工況下,開展外磁場B=0.008,0.017,0.025,0.034,0.052 和0.069 T 情況下的MPDT 工作過程仿真,分析不同磁場強度下的電子數密度分布和推力器性能參數,并進行仿真結果與實驗結果對比,驗證仿真模型和方法的有效性.

3.3.1 不同外磁場的電子數密度分布

不同外磁場作用下,MPDT 穩態工作時的電子數密度分布如圖7 所示.

圖7 不同磁場強度下的電子數密度分布Fig.7.Electron number density distribution with different applied magnetic field.

從圖7 可知,在B=0.008 T 情況下,電子僅在陰極下游的軸線附近形成高密度區域,未能形成高效的陰陽極電弧放電;隨著磁場強度的增大,中軸線高密度區域逐漸向徑向移動,在外磁場強度B=0.017—0.340 T 之間時,電子數密度出現兩個高密度區域,分別位于軸線附近和推力器出口下游區域,構成穩定的放電電弧;外磁場強度繼續增大,軸線附近的高密度區域逐漸縮小,當磁場B＞0.052 T 時,中心軸線的高密度區域基本不存在.分析認為:在外磁場強度較小時,推力器的自感磁場比重較大,自感磁場主要為周向磁場,約束電子在軸線附近聚集;隨著外磁場的增加,自感磁場的比重降低,外磁場對電子的約束成為主要部分,外磁場將使得電子沿著磁力線向下游運動.

3.3.2 外磁場對推力器性能的影響

根據不同外磁場作用下的仿真結果,計算了MPDT 穩態工作時的放電電壓U、推力F和效率η,如圖8 所示.

由圖8 可知,當B＜ 0.053 T 時,放電電壓隨外磁場增大而降低,B≥ 0.052 T 時放電電壓逐漸變大.其原因是陰極表面電荷密度決定了陰極電勢即放電電壓的大小,而陰極表面電子的運動受磁場影響嚴重:當外磁場較小時,自感磁場對電子的約束占據主導,電子向陰極聚集,增大了陰極表面電荷密度,放電電壓升高;隨著外加磁場增大,電子逐漸跳出自感磁場約束,遠離陰極,使得放電電壓降低;當外磁場進一步增大時,外加磁場對電子約束增強,電子難以跨越陰極磁力線沿徑向擴散,使得陰極表面電荷密度增大,放電電壓再次升高.

圖8 MPDT 性能隨外磁場的變化Fig.8.MPDT performance changes with applied magnetic field.

不同磁感應強度條件下的推力結果表現為隨磁感應強度B近似線性增大,與Myers[20]的實驗結論一致.在B=0.069 T 情況下的推力值為0.98 N,略小于Myers[20]的實驗結果1.16 N,誤差產生的原因可能是本文仿真忽略了氣動推力的貢獻;推力器的效率同樣表現為隨磁感應強度B增大而增大的趨勢.

3.3.3 與實驗測試數據對比分析

Myers[20]實驗測量了不同外加磁場情況下的MPDT 推力器的比沖、效率等性能參數,為了進一步驗證本文仿真程序的準確性,將仿真結果與Myers[20]的實驗結果進行對比,如圖9 所示.

圖9 仿真與實驗結果對比Fig.9.Comparison of simulation and experimental results.

由圖9 可知,仿真獲得的比沖、效率等性能參數隨外加磁場強度的變化趨勢與實驗結果一致.計算比沖值低于實驗結果約20%,其原因是計算比沖由Isp=F/m˙ 計算得來,在本文仿真中未考慮氣動力對總推力的貢獻,使得計算推力性能低于真實情況,并導致計算比沖低于實驗結果.計算效率值相比實驗結果有約50%的誤差,其原因可能是計算效率與推力和比沖的乘積成正比,計算推力和計算比沖均低于實驗結果,導致計算效率誤差進一步增大.

雖然推力、比沖和效率的仿真結果因忽略氣動加速貢獻略小于真實值,但誤差不是很大,并且能準確反映工況對性能的影響,故認為,本文建立的模型和方法在高磁場強度工況下有較好的準確性和可信度.

3.4 陰極電流對MPD 工作過程及性能的影響分析

陰極電流是影響MPDT 推力器性能的另一個重要參數,在外磁場B=0.025 T,流量=0.1 g/s工況下,開展陰極電流I=750,1000,1250,1500,1750 和2000 A 情況下的MPDT 工作過程仿真,分析陰極電流對電子數密度分布和推力器性能的影響.

3.4.1 電子數密度分布

不同陰極電流情況下,MPDT 穩態工作時的電子數密度分布如圖10 所示.

由圖10 可知,推力器內部的電子數密度隨著電流的增大而增大;對于B=0.025 T,流量=0.1 g/s 工況,750—2000 A 的陰極電流情況,都在陰極尖端下游位置都出現了電子高密度分布區域.陰極電流增大,高密度區域也逐漸擴大,具體表現為隨著電流增大,外磁場對放電等離子體束的約束效果逐漸減弱,低電流工況的等離子體束流發散角明顯小于高電流工況,其原因是電流增大,電子數密度增大,電子間的庫侖碰撞增多,會使得更多的電子沿跨越磁力線沿徑向擴散.

圖10 不同陰極電流下的電子數密度分布Fig.10.Distribution of electron number density under different cathode current.

3.4.2 陰極電流對推力器性能的影響

根據不同陰極電流條件下的仿真結果,計算了MPDT 穩態工作時的放電電壓U、推力F和效率η,如圖11 所示.

圖11 MPDT 性能隨陰極電流的變化Fig.11.MPDT performance changes with cathode current.

從圖11 可知,放電電壓隨陰極電流的增大而增大,原因是陰極電流越大,陰極表面的電子數密度越大,陰極表面電荷增加,放電電壓增大;另外,自感磁場隨陰極電流增大而增大,并進一步約束電子向陰極聚集,使得陰極表面電荷增加,放電電壓增大.

在750—2000 A 陰極電流范圍內,推力隨電流增大而近似線性增大.在I=750—2000 A 的電流工況下,計算區域內的自感磁場強度為0.001—0.040 T,量級上與外磁場相當,但自感磁場主要集中于陰極附近,在其余大部分區域外磁場明顯強于自感磁場,故此時MPDT 的推力由渦旋加速機制產生,正比于電流大小.

推力器的效率同樣表現為隨陰極電流I增加而增大的趨勢,這與MPDT 低功率工況效率低,高功率工況效率高的實驗結果一致.

4 結論

本文建立了基于縮比模型的MPDT 工作過程的粒子仿真模型,研究了MPDT 點火啟動時的放電電弧建立過程,以及外磁場強度和陰極電流對MPDT 性能參數的影響,并通過仿真結果與已有實驗數據對比,驗證了仿真模型的準確性.點火啟動過程研究表明,放電電弧建立后,在陰極下游和陽極附近形成兩處高電子數密度區并連通以維持穩定的放電電弧,離子主要集中在陰極下游并沿軸向噴射形成等離子體噴流.不同外磁場強度下的結果表明,提高外磁場強度有利于構建陰陽極放電電弧,隨著外磁場強度由0.008 T 提升至0.069 T,MPDT 的推力產生機制逐漸由自身場加速轉變為渦旋加速,推力和比沖性能隨外磁場強度增加近似線性增大.不同陰極電流的結果表明,在B=0.025 T,I=750—2000 A 條件下,在推力器內部絕大部分區域,外磁場強度大于自感磁場,推力產生機制為渦旋加速,推力隨電流增加而近似線性增加.MPDT 效率總體表現為隨外磁場強度和陰極電流增大而增大的趨勢,符合其適用于高功率工況的特性.