劉 祺,楊 磊,黃玉平,鄭再平,趙 絮

0 引 言

脈沖等離子體推力器(pulsed plasma thruster, PPT)的概念最早在20世紀30年代被提出，自1960年開始進入研究階段，至今已有幾十年的發(fā)展歷史.在這幾十年的研究中，PPT已發(fā)展成為固、液、氣3種工質(zhì)，平行軌道、同軸、T管等電極結構，可外加磁場的多種分支類型[1].作為世界上第一臺執(zhí)行空間任務的電推力器，PPT在姿態(tài)精確控制與軌道修正任務中具有較大的優(yōu)勢和潛力.美國空軍研究實驗室通過樣機XPPT-1的試驗發(fā)現(xiàn)，他們利用電子顯微鏡觀測到發(fā)射的40±3%的推進劑粒子僅貢獻1%的推力[2]，因此PPT性能的提升潛力很大，值得進一步研究.

PPT的研究主要有實驗研究與仿真研究，研究主要集中在推力測量、推力器等離子體探針及光譜的診斷、推力器放電特性的仿真與實驗研究等.目前對PPT的仿真研究主要使用MHD方法以及機電模型，難以直接對電離過程進行仿真，也難以直觀反映PPT放電室內(nèi)粒子密度等參數(shù)分布.隨著電推進技術的發(fā)展與研究的深入，以粒子為對象的方法——粒子網(wǎng)格-蒙特卡洛碰撞(PIC-MCC)方法已應用于霍爾推力器和離子推力器的仿真工作中.而對于脈沖等離子體推力器的仿真研究，多采用前述兩種方法，PIC-MCC方法在PPT工作過程仿真中的研究仍非常有限.雖然有限元軟件COMSOL等應用廣泛，已應用于某些非平衡放電的仿真中，但對于真空放電問題，至今尚未開發(fā)基于PIC方法的計算模塊.

PPT內(nèi)等離子體稀薄，粒子并非連續(xù)，文獻[3]中的實驗結果也表明粒子團不連續(xù).因此基于連續(xù)流體假設的MHD方法并不是計算此類問題的最佳選擇.PIC方法仿真非常接近實際過程，這是PIC方法的優(yōu)點，但這也帶來了一些問題.PPT中的粒子密度高，因此計算量也隨之增大.PPT中電磁場及等離子體密度都是非穩(wěn)態(tài)的，變化迅速，這也給仿真帶來困難.應對這個問題，NEUDORFER等[4]在對PPT放電室的PIC仿真計算中進行了很多簡化，并利用反饋控制的思想加強計算的穩(wěn)定性和收斂性.

NEUDORFER等[4]利用PIC方法對PPT進行了三維模擬，不考慮電離碰撞，將離子與電子在相同位置加入到計算域中.然而，電離碰撞是推力器工作時非常重要的一個過程，也是PIC方法與MHD方法的不同點之一.不考慮電離過程則失去了粒子仿真的部分意義.HUANG等[5]利用電路模型、燒蝕模型，基于等離子體片質(zhì)量不變的假設，進行了仿真研究，他在計算中加入電離模型，依據(jù)電流片能量是否達到一定值來判斷整個電流片是否電離，但是假設電流片形狀質(zhì)量不變，這與實驗結果不符，因為實際過程中基本不可能出現(xiàn)整個電流片中的粒子全部同時電離的情況.ANTIPOV等[6]認為PPT工作初期等離子體密度很低，等離子體密度分布對電磁場影響小，因此他們僅對電磁場中的單粒子運動軌跡進行研究.事實上，雖然等離子體是由單個電子和離子組成的，但僅研究單個粒子的運動不足以體現(xiàn)等離子體團的特性.

本文利用PIC-MCC方法，結合電路模型，建立LES-6 PPT的放電室內(nèi)粒子網(wǎng)格模型，得到燒蝕的推進劑粒子在放電室內(nèi)的分布，電流電阻的變化情況，以及PPT放電過程中等離子體密度分布情況.對加入與不加入粒子預分布這兩種條件下的電離過程進行仿真計算，并將二者結果進行對比.

1 模型建立

1.1 電路模型

PPT的電源為脈沖電源，工作時，首先使儲能電容器充電到它的工作電壓，然后放電使點火塞點火[1].將這一過程簡化為RLC無源脈沖放電回路，方程寫為[5]：

(1)

其中，I、Q、U分別為放電電流、電容帶電量以及電容器的電壓，L、R分別為電路的總電感和總電阻，C、U0分別為電容和電容的初始電壓.電流可由經(jīng)驗公式計算，經(jīng)驗公式中需提供等離子體電阻[7].

1.2 碰撞模型

PPT的推進劑是聚四氟乙烯，包含兩種元素——碳和氟.由于其碰撞過程較為復雜，所以推力器在工作過程中可能產(chǎn)生多種價態(tài)的離子[8].實驗測試已證明，在推力器放電室中電離產(chǎn)生的離子主要以一價離子為主[9]，在燒蝕型PPT的仿真研究中，也是通常只考慮一價離子，忽略多價離子[10].因而我們將推力器放電時推進劑的分解簡化為[11]：

aC2F4→bC+cC++dF+eF+

(2)

電離碰撞發(fā)生在電子與原子、電子與一價離子之間，這里主要討論電子與原子之間的電離碰撞.當電子的能量達到電離能時，此電子與原子可能發(fā)生電離碰撞.一個電子與一個原子發(fā)生電離碰撞后，產(chǎn)生一個低能的二次電子，原子失去一個電子變成一價離子，原初電子發(fā)生能量損失，速度重新分配.碰撞后原初電子的能量Escat,e為初始能量Einc,e減去二次電子的能量Ese和電離能Eionize：

Escat,e=Einc,e-Ese-Eionize

(3)

電離碰撞后產(chǎn)生離子的速度按照碰撞前原子的速度給定，產(chǎn)生二次電子的能量為[12]

(4)

其中B(Einc,e)為已知函數(shù).

陰極為金屬材質(zhì)，當電子到達陰極，視為被漫反射；粒子到達右邊界時將被視為離開計算域.

1.3 電磁場計算

電勢與電場的關系寫為[13]：

E=-φ

(5)

E為電場強度，φ為電勢.電場與等離子體電荷密度的關系表示為：

E=ρ/ε0

(6)

其中，ε0是真空介電常數(shù)，ρ為凈電荷密度：

ρ=e(ni-ne)

(7)

其中，e為單位電荷量，ni和ne分別表示離子和電子的粒子數(shù)密度.將式(5)代入式(6)，可得泊松方程：

2φ=-e(ni-ne)/ε0

(8)

將泊松方程寫為直角坐標系下偏微分方程的形式(二維)：

(9)

為了使用計算機對方程進行數(shù)值求解，將式(9)寫成有限差分格式：

(10)

本文采用超松弛迭代方法(SOR)進行電場的求解，迭代公式為

(11)

磁場的計算方法與電場相似，也采用SOR迭代方法，所不同的是，磁場計算考慮感生磁場，加入電流對磁場的影響.在電勢的計算中，陰極邊界的電勢由電源電壓的變化而定；放電室右邊界為外推邊界.

2 幾何模型及參數(shù)選取

2.1 推力器模型

本文選取美國林肯實驗通信衛(wèi)星6號平臺LES-6的東西位保PPT作為研究對象進行仿真計算，其相關幾何參數(shù)及電路參數(shù)如表1所示，脈沖燒蝕質(zhì)量來自美國宇航局(NASA)的實驗數(shù)據(jù)[14-15].

表1 LES-6 PPT的幾何參數(shù)及電路參數(shù)Tab.1 Geometric and circuit parameters of LES-6 PPT

圖1為LES-6 PPT的基本結構，它是平行軌道電極型PPT，電極外接儲能電容器.本文研究的計算域為圖中的放電室區(qū)域.

圖1 LES-6 PPT基本幾何結構示意圖Fig.1 The structure of LES-6 PPT

2.2 仿真參數(shù)的選取

仿真中選取的主要參數(shù)列于表2.帶電粒子的權重為5×108，即一個仿真粒子代表5×108個真實粒子；由于原子的數(shù)量遠高于帶電粒子的數(shù)量，故將原子的權重取為1011.由于PIC方法的計算量非常大，本文使用了PIC方法常用的兩個簡化方法，目的分別是減小內(nèi)存占用及加快計算收斂速度：減小重粒子(離子和原子)質(zhì)量和增大真空介電常數(shù)ε0的方法.將重粒子質(zhì)量由原來的M減小為M′，可以使計算域中速度最慢的原子運動加快，從而縮短計算時間.

時間步長Δt的選取依據(jù)德拜長度，空間步長Δx的選取依據(jù)等離子體振蕩頻率.

Δx<λd

(12)

(13)

其中λd為德拜長度，ωp為電子振蕩頻率.

表2 數(shù)值計算的參數(shù)Tab.2 Numerical parameters

3 結果分析

3.1 燒蝕中性粒子分布

推力器開始點火工作時，電極與推進劑表面之間產(chǎn)生微量放電，放電產(chǎn)生的電子在電場中向陽極方向加速運動.這些加速后的電子與推進劑表面碰撞，使工質(zhì)分解和電離，產(chǎn)生離子和更多的電子.這是PPT放電時的真實過程.而在計算電離碰撞時，考慮的是碳、氟原子分別與電子的碰撞，而非電子直接與特氟龍進行碰撞.因此，將此處進行簡化，在中性粒子進入計算域時，直接成為碳原子和氟原子.推進劑燒蝕0.5 μs后兩種原子的分布如圖2所示.

由于特氟龍中含氟元素的量高于碳元素含量，因此從圖2可以看出，氟原子的密度比較大，兩種原子數(shù)密度分布輪廓基本一致.

圖2 燒蝕產(chǎn)生碳、氟原子的數(shù)密度分布Fig.2 The amount densities of ablated carbon and fluorine atoms

3.2 等離子體電阻及電流變化規(guī)律

根據(jù)文獻[16]的放電模型進行計算，得到3 μs內(nèi)電路總電阻及電流隨時間的變化.從圖3中可以看出，等離子體電阻在初始時非常大，之后隨著推進劑的電離，等離子體在電極之間運動，使總電阻減到極小.由此可知，從PPT點火到開始放電的過程是非常迅速的.

圖3 放電過程中等離子體電阻及放電電流隨時間的變化Fig.3 The plasma resistance and discharge current during the discharge process

從圖3電流的變化可以看出，在放電過程中，回路為欠阻尼狀態(tài)，電流曲線呈振蕩衰減趨勢.電流曲線的第一個峰較大，電流產(chǎn)生的感生磁場也較大，因此等離子體加速主要發(fā)生在這一階段.

3.3 帶電粒子分布

本文的PIC仿真中包含3種帶電粒子：電子、C+和F+.分別在兩種條件下進行了仿真.Case 1：在計算前計算域內(nèi)沒有布置粒子，直接進行仿真計算；Case 2：輸入粒子在空間中的預分布，在預分布的基礎上進行計算.預分布粒子的布置方法為：將粒子以基于壁面溫度的麥克斯韋速度分布從左邊界入射，計算域中的粒子數(shù)收斂后，將得到的粒子分布作為預分布輸入到Case 2的計算程序中.

在t=0.5 μs時，在計算之前沒有進行粒子預分布的情況下，3種粒子的數(shù)密度分布云圖如圖4所示.從圖4(a)可以看出，電子從陰極發(fā)射，部分電子在電場作用下向陽極運動，大部分電子仍集中在陰極附近靠近推進劑表面的位置.圖4(b)和圖4(c)分別為C+和F+的粒子數(shù)密度分布情況，比較兩圖可以看出，雖然放電室內(nèi)氟元素所占比例應高于碳元素，但由于F原子的第一電離能高于C原子，仿真中電子沒有充分加速，能夠達到F原子電離能的電子數(shù)量較少，所以此時刻空間中的F+很少.

對于Case 2，從圖5可以看出，在計算之前將部分粒子預分布在計算域中，t=0.5 μs時帶電粒子密度比Case 1中的密度高出兩個量級.t=0.5 μs時刻在PPT內(nèi)放電是很劇烈的，離子密度應達到較高的值.根據(jù)文獻[17-18]，PPT羽流中離子數(shù)密度峰值的量級為1021/m3，因此Case 2中的計算結果與目前實驗及仿真的粒子密度分布更為接近.考慮PPT的工作過程，在推進劑大規(guī)模電離之前，少量電子首先產(chǎn)生微量電離，此點火過程約持續(xù)幾微秒，因此在推力器開始放電之前，空間中已有等離子體的分布.由此可知，在進行PPT放電的仿真時，空間中應預分布等離子體.

對比圖5(a)與圖4(a)，由于Case 2中空間分布的粒子與電子發(fā)生碰撞促進了電子的遷移運動，所以圖5(a)中的電子擴散范圍較大.相應的，也就有更多的原子被電離，電離產(chǎn)生的碳離子和氟離子的數(shù)密度分布云圖分別如圖5(b)、(c)所示.

圖4 Case 1 t=0.5 μs時刻帶電粒子數(shù)密度分布云圖：Fig.4 Case 1: The amount density of charged particles at 0.5 μs

圖5 Case 2 t=0.5 μs時刻帶電粒子數(shù)密度分布云圖：Fig.5 Case 2: The amount density of charged particles at 0.5 μs

由圖3電流隨時間的變化曲線可以看出，在t=0.5 μs和t=1.5 μs時，電流的方向是相反的.圖6和圖7分別為t=1.5 μs時帶電粒子的分布云圖.從圖6(a)與圖7(a)可以看出，由于此時陽極電勢低于陰極，所以電子在陽極產(chǎn)生并入射，在電場的作用下向陰極運動.電子在靠近陽極處較集中，在y=0.028 5 m 附近也相對集中，Case 1與Case 2電子的分布情況基本一致.在圖6(b)、(c)與圖7(b)、(c)中，t=1.5 μs時離子與電子的分布情況相似，Case 1與Case 2離子的分布情況也相似.由此可以看出，雖然Case 1與Case 2的初始條件不同，但在t=1.5 μs 時都具有相似的分布特征.相對而言，在計算前未進行預分布的算例(Case 1)中的粒子數(shù)密度普遍低于1020/m3，而在Case 2有預分布粒子的條件下粒子的數(shù)密度更接近粒子密度的實驗值.

圖6 Case 1 t=1.5 μs時刻帶電粒子數(shù)密度分布云圖：Fig.6 Case 1: The amount density of charged particles at 1.5 μs

圖7 Case 2 t=1.5 μs時刻帶電粒子數(shù)密度分布云圖:Fig.7 Case 2: The amount density of charged particles at 1. 5μs

4 結 論

本文針對PPT的放電過程，利用PIC-MCC方法結合電路模型，進行電離仿真研究，得到了一些初步結果.通過電路模型的計算結果得到電流、電阻變化規(guī)律；建立LES-6的放電室內(nèi)粒子網(wǎng)格模型，得到燒蝕的推進劑粒子在放電室內(nèi)的分布以及PPT放電過程中等離子體密度分布情況，對采用PIC方法研究PPT的工作過程具有一定的意義.本文通過在不加入粒子預分布與加入粒子預分布的兩種條件下進行仿真計算，對比粒子分布情況可知，加入粒子預分布的算例的密度計算結果更接近實驗結果.參考實際的PPT工作過程，在放電之前推力器內(nèi)是有等離子體存在的，所以在仿真中應該進行粒子的預分布.由于PPT放電過程中電磁場變化迅速、各參量的瞬態(tài)變化，計算易產(chǎn)生不收斂的情況，故本文做了較多的簡化處理，僅得到了初步的計算結果.下一步將進行更深入的研究，使仿真更接近實際的放電過程.

參 考 文 獻

[1] 毛根旺, 唐金蘭. 航天器推進系統(tǒng)及其應用[M]. 西安：西北工業(yè)大學出版社, 2009.

[2] SPANJERSG G, LOTSPEICH J S, MCFALL K A, et al. Propellant losses because of particulate emission in a pulsed plasma thruster[J]. Journal of Propulsion & Power, 1998,14(4):554-559.

[3] RYOSUKE N, KIMIYA K, YOSHIHIRO A. Study on plasma acceleration in an ablative pulsed plasma thruster[C]//Proceedings of the 43rdAIAA/ASME/ SAE/ASEE, Joint Propulsion Conference and Exhibit. Washington D.C.: AIAA 2007-5226, 2007.

[4] NEUDORFER J，STINDLB T，SCNNEIDER R. Three-dimensional particle-in-cell simulation of a pulsed plasma thruster: modeling and challenges[C]//The 32ndInternational Electric Propulsion Conference. Wiesbaden Germany: IEPC-2011-116,2011.

[5] HUANGT, WU Z, LIU X, et al. Modeling of gas ionization and plasma flow in ablative pulsed plasma thrusters[J]. Acta Astronautica, 2016, 129:309-315.

[6] ANTIPOVA A, BOGATYY A B. Simulation of the motion of charged particles in an ablative pulsed plasma thruster at the initial stage of the discharge[J]. Procedia Engineering, 2017, 185:61-68.

[7] 鄧健. 脈沖等離子體推進器脈沖電源及推進劑燒蝕特性的實驗研究和數(shù)值分析[D]. 南京理工大學, 2013.

[8] YANG L, ZENG G, TANG H, et al. Numerical studies of wall-plasma interactions and ionization phenomena in an ablative pulsed plasma thruster[J]. Physics of Plasmas, 2016, 23(7):716-15.

[9] SCH?NHERRT, NEES F, ARAKAWA Y, et al. Characteristics of plasma properties in an ablative pulsed plasma thruster[J]. Physics of Plasmas, 2013, 20(3):1310-563.

[10] BOYDI D, KEIDAR M, MCKEON W. Modeling of a pulsed plasma thruster from plasma generation to plume far field[J]. Journal of Spacecraft & Rockets, 1971, 37(3):399-407.

[11] GATSONISN A, YIN X. Hybrid (particle-fluid) modeling of pulsed plasma thruster plumes[J]. Journal of Propulsion & Power, 2015, 17(5):945-958.

[12] VAHEDI V, SURENDRA M. A Monte Carlo collision model for the particle-in-cell method: applications to argon and oxygen discharges[J]. Computer Physics Communications, 1995, 87(1-2): 179-198.

[13] BIRDSALLC K. Particle-in-cell charged-particle simulations, plus Monte Carlo collisions with neutral atoms, PIC-MCC[J]. IEEE Transactions on Plasma Science 1991, 19(2): 65-85.

[14] HAMIDIAN J P , DAHLGREN J B . Testing and evaluation of the LES-6 pulsed plasma thruster by means of a torsion pendulum system[R]. NASA-CR-135940, 1973.

[15] 楊樂. 脈沖等離子體推力器工作過程理論和實驗研究[D]. 長沙：國防科學技術大學, 2007.

[16] YANGL, LIU X Y, WANG S Y, et al. Theoretical and numerical analysis of discharge characteristics in pulsed electromagnetic accelerators[J]. Advanced Materials Research, 2013, 765-767:805-808.

[17] SOLBESA, THOMASSEN K, VONDRA R. Analysis of solid teflon pulsed plasma thruster[J]. Journal of Spacecraft & Rockets, 1970, 7(7):1402-1406.

[18] BOYDI D, KEIDAR M, MCKEON W. Modeling of a pulsed plasma thruster from plasma generation to plume far field[J]. Journal of Spacecraft & Rockets, 2000, 37(3):399-407.

[19] YANGL, LiU X Y, et al. Different ablation models for the wall-plasma interaction process in pulsed plasma thruster[J]. High Voltage Engineering, 2013, 39(9): 2301- 2308.

[20] KEITHI T. Exhaust velocity studies of a solid teflon pulsed plasma thruster[J]. Journal of Spacecraft & Rockets，1972，9 (1) ：61-64.