多級會切場推力器振蕩特性模擬研究

，，，

哈爾濱工業大學 先進動力研究所，哈爾濱 150001

多級會切場推力器振蕩特性模擬研究

牛翔，伍環，劉輝*，于達仁

哈爾濱工業大學 先進動力研究所，哈爾濱 150001

目前多級會切場推力器(High Efficiency Multi-Plasma Thruster，HEMPT)數值模擬研究采用的主要手段是全粒子(Paticle In Cell，PIC)模擬方法，但是這種全粒子模擬對計算機的要求較高，因此文章通過建立HEMPT的一維流體模型，再現HEMPT的振蕩現象，對放電電流和等離子體參數之間的關聯性及等離子體參數分布規律進行研究。沿著兩條典型的電子傳導路徑,利用MATLAB中Simulink模塊建立一維流體模型。結果分析表明，推力器內部等離子體參數的振蕩周期與放電電流振蕩周期一致。推力器的主電離區在兩個磁尖端之間，電勢降主要集中在出口附近。

多級會切場推力器；一維流體模型；振蕩；放電電流；電子傳導；電勢分布

與化學推進相比，電推進具有高比沖、小推力等特點。霍爾推力器是一種常見的電推力器，離子與壁面劇烈的碰撞是導致其壽命下降的重要原因，而會切場的磁場位型可以避免這一點，因而其具有長壽命的特點[1]。此外試驗表明，會切場推力器可實現推力大范圍連續可調[2]。它的工作原理如圖1所示。

圖1 多級會切場推力器結構原理Fig.1 Illustrative diagram of HEMPT

多級會切場推力器(High Efficiency Multi-Plasma Thruster，HEMPT)內部存在一系列的復雜物理過程，在研究這些物理過程中，數值模擬因其成本低，能夠直接獲得難以測量的參數等優點而被廣泛應用。基于微觀機理的全粒子(Paticle In Cell，PIC)模擬是研究推力器內部復雜物理過程的一種重要手段，文獻[3]用PIC模擬了霍爾推力器與HEMPT推力器，并得到了HEMPT推力器內部電子與壁面碰撞少的結論。

PIC作為一種以粒子為對象的模擬方法，能夠反映等離子體參數在空間的二維分布。文獻[4]就利用PIC模型關于發射電子的徑向分布對等離子體參數空間分布的影響進行了研究。但是目前PIC模擬存在對計算機性能要求較高、計算量較大等難點。與PIC相比，一維流體模型雖然不能反映等離子體參數的二維分布，但是算例的運行時間明顯減少。另外，文獻[5]在穩態等離子體推進器(Stationary Plasma Thruster，SPT)的研究中也驗證了一維流體模型研究放電電流振蕩特性的準確性，并且利用該模型得到的結論與胡俊鋒等的試驗結論相符[6]。因此在放電電流振蕩模擬研究方面，使用一維流體模型可以保證在一定準確程度的情況下，顯著降低算例運行時間。流體模型中電子采用穩態假設，可以很快得到推力器內部參數值，極大地減小了運算量。文獻[7-10]利用流體模型對霍爾推力器內部呼吸模式振蕩進行了模擬，文獻[11-12]不僅利用流體模型對霍爾推力器的模式轉換進行研究，還利用混合直接動力學模擬以獲得更為準確的離子速度分布。

劉輝等利用PIC研究表明在有兩個磁尖端的多級會切場推力器內部，電子有兩條典型的電子路徑[13]。當電子經過出口磁尖端進入通道內部后，在中部尖端區域電子路徑開始分叉，有一部分電子會沿著磁感線運動至磁尖端，然后跨越磁感線繼續運動，這類電子稱之為捕獲電子；還存在一部分電子并沒有沿磁感線運動，而是直接通過通道中央從而到達陽極，這部分電子稱之為漏失電子。漏失電子的產生是由于一部分電子在跨越第一個磁尖端后變成高能電子，在做拉莫爾回旋運動時，所受的離心力較大，因而能夠擺脫后面磁尖端的束縛，從中央中軸線直接逃逸。本文沿著兩條典型電子路徑建立一維流體模型，在抓住電子主要的傳導路徑基礎上簡化了模型，減小了運算量。本模型不僅可以研究會切場推力器的低頻振蕩，還可以通過計算通道內部等離子體參數的軸向分布，為更精確的模型提供初始分布。但是，模型計算得到的等離子體參數分布是通道截面的平均值，因而無法獲得在某一截面處的更為精確的等離子體參數分布，這是目前一維流體模型的不足之處。另外，由于一維流體模型主要用來研究會切場推力器的低頻振蕩，為了縮短收斂時間，這里沒有用高維度的模型進行研究。

1 流體模型的建立過程

本文通過建立流體模型來對等離子體內部各個參數進行模擬。電子路徑如圖2所示，通過沿著兩條電子路徑建模可以將問題簡化為一維。根據電子、原子、離子的運動過程作出如下假設：

1)為了避免雙流體模擬帶來的時間尺度問題，流體模擬采用對電子進行穩態假設，即忽略電子的時間偏導項。

2)電子在運動過程中，始終沿著磁感線運動。

3)模型僅考慮氙原子的一次電離碰撞，不考慮高價氙離子的產生、原子激發、庫倫碰撞。

圖2 電子路徑Fig.2 Electron path diagram

對玻爾茲曼方程取速度的零階矩，可以得到原子、離子的連續性方程；對玻爾茲曼方程取速度的一階矩可以得到離子動量方程：

(3)

式中：n為等離子體體密度；na為原子密度；Va為原子速度，取常數Va=200 m/s；Vi為離子速度；mi為氙離子的質量，取值2.2×10-25kg；e為元電荷，取值為1.6×10-19C；E為電場強度；βi為電離率。

回路方程及電流方程如下：

式中：I為放電電流；l為通道長度，取值為0.1 m；A為通道面積，取值為5.3×10-4m2；Lc為通道電感，取值為10-6H；R為內阻，取值為10Ω；Vex為電子速度。

由于假設電子一直沿著磁感線運動，因此電子運動方程如下：

式中：me為電子質量，取值為9.1×10-31kg；Te為電子溫度；B為磁感應強度；μ為比例系數，取值為1/3；νeff為有效電子碰撞頻率：

式中：βm為電子與氙原子的動量碰撞交換系數；αB為波母系數，取值為1/32；ωce為電子自旋頻率，取值(eB)/me。

由圖2可知，通道內存在兩條電子路徑，為了將二維問題簡化為一維問題，提出下面合并電子路徑的公式：

式中：α為給定的參數，用來表示在中部電子通過中央路徑引起的電流密度與總電流密度之比。jex1表示逃逸電子引起的電子電流密度，jex2表示捕獲電子引起的電子電流密度。合并式(6)、(8)可得：

式中：νeff1為逃逸電子有效碰撞頻率；νeff2為捕獲電子有效碰撞頻率。

由之前的電子處于穩態可忽略偏導項的能量方程：

(10)

式中：Δε為單次電離碰撞所引起的能量損失，由于只考慮一次碰撞，這里取值為12.1 eV。右邊第一項為歐姆加熱項，第二項為電離碰撞項，由于忽略電子碰壁及熱傳導造成的損失，因此在第二項前面加一個系數γi，使模擬結果更接近真實情況，這里取系數γi為3，文獻[8]也做了相同處理。

邊界條件處理如下：陽極位置的原子密度，離子密度和離子速度分別為2×1019m-3，6×1017m-3，2×103m/s。陰極處的電子溫度為1 eV。模型的初始條件即初始原子和離子密度分布，離子速度分布分別如圖3、圖4所示。

圖5為計算所用磁場。利用捕獲電子與漏失電子數量之比α反映兩種不同類型電子的相對大小。非尖端區不考慮波母傳導即計算電子有效碰撞頻率時不包括電子自旋頻率的平方項。而在尖端區則考慮該項，尖端區與非尖端區的分界是以磁場強度導數為0定義的。

圖3 初始原子、離子密度分布Fig.3 Initial atom and ion distribution

圖4 初始離子速度分布Fig.4 Initial ion velocity distribution

圖5 磁場分布Fig.5 Magnetic field distribution

2 計算結果及分析

在電壓Us=300 V，α=0.1的條件下模擬推力器的振蕩問題。如圖6所示，放電電流產生劇烈振蕩，時均放電電流約為7 A，電流峰值約為10.8 A，振幅約為5.8 A，振蕩頻率約為40 kHz。

圖6 電流的振蕩曲線Fig.6 Current oscillation

圖7、圖8分別為電子溫度及歸一化的方程源項的分布情況，其中點畫線為磁分界面，虛線為電流振蕩為峰值的時刻。方程源項是表征推力器內部電離程度的參數。源項越大，電離程度就越大。從圖8中可以看出方程源項在時間上具有周期性，且周期與放電電流的振蕩周期一致；在空間上，方程源項在歸一化的軸向位置為0.4處達到最大，所以這個位置應為推力器的主電離區。

從時間上可以看出，當源項恰好處于極小值時，如在0.17 ms左右，振蕩電流達到峰值，這是之前主電離區大量電離的結果。之后從圖8中可以看出，電離程度變得很弱，因此電流也隨之下降。在一個振蕩周期里，主電離區內存在兩個源項極大值，且其數值一小一大，也就是說電離會有兩次不同程度增強。在第一次電離增強區域，電流隨之增加，但是由于電離程度較小，所以電流所達到的極大值小于峰值;而第二次由于電離程度更大，導致第二次電流達到峰值，隨之進入下一個周期。

從空間上看，主電離區的位置位于兩個磁尖端之間，與高電子溫度區域一致。在主電離區上游，由于電子溫度較低，電離程度下降，所以源項較小。而在主電離區下游，一方面電子溫度比較小，另一方面大部分原子已經被電離，所以在這一區域電離程度也很弱。

值得注意的是，電子溫度分布也存在類似于源項分布的周期性，并且在一個周期內存在一大一小兩個極值。

圖9為瞬時電場分布。可以看出，在出口磁尖端附近電場強度總是處于極大值，這主要是由于出口磁尖端存在磁鏡效應，電子跨越磁感線的運動很難實現，因而在出口處存在較大的電勢降。除此之外，在通道中部至出口處存在隨時間周期性變化的場強，而通過與電子溫度的分布圖7對比可以發現，在場強取得極大值處電子溫度均較小，這可能是由于局部電子溫度較低導致電離不充分而形成的電勢降。

圖7 電子溫度分布Fig.7 Electron temperature distribution

圖8 歸一化方程源項Fig.8 Normalized source term

值得注意的是，在中部磁尖端之前存在著較小的電勢降，而非磁尖端。而通過圖2可以看出，中部磁尖端電離比較強，因而電勢降提前可能是由于中部磁尖端附近電離增強造成的。

原子、離子密度分布如圖10、11所示。可以看出，原子密度在歸一化的軸向位置為0.4左右時會有振蕩現象發生，而粒子密度分布從歸一化軸向位置為0.4后也有明顯振蕩現象。從放電電流振蕩圖像中容易得到，多級會切場推力器的振蕩頻率在5 kHz左右，這與霍爾推力器內部低頻振蕩頻率類似，研究表明霍爾推力器內部的振蕩可能是由于原子供給與電離過程的不匹配造成的[14-17]。當局部電離增強時離子密度上升，原子密度下降，從而放電電流上升。此后離子經電場加速噴出，使離子密度下降，待原子供給充足后，電離繼續增強，進入到下一輪循環。所以這也可能是多級會切場推力器內部低頻振蕩發生原因。

從圖10和圖11通道中部就可以看出原子密度下降的同時,離子密度上升，同時對應放電電流上升并且達到最大值。

圖10 原子密度分布Fig.10 Atom density distribution

圖11 離子密度分布Fig.11 Ion density distribution

3 結束語

本文沿著兩條典型電子路徑建立了HEMPT的一維流體模型，復現了放電電流振蕩現象。模擬結果分析表明，HEMPT放電電流與等離子體參數振蕩頻率一致，推力器的主電離區在兩個磁尖端之間，主要電勢將集中在通道出口附近。利用建立的一維流體模型可以在保證一定精度的前提下，較快地得到推力器內部振蕩參數及等離子體參數分布，為推力器振蕩特性的研究提供便利。但是一維流體模型各個位置處均采用了截面平均的假設，這在一定程度上限制了模擬精確程度的提高，因此后續可以通過建立HEMPT的二維流體模型，在更高精度上對HEMPT內部振蕩現象進行研究。

References)

[1] GILDEA S R，MATLOCK T S，MARTNEZSNCHEZ M，et al.Erosion measurements in a low-power cusped-field plasma thruster[J].Journal of Propulsion & Power，2013，29(4):906-918.

[2] GENOVESE A，LAZURENKO A，KOCH N，et al.Endurance testing of HEMPT-based ion propulsion modules for small GEO[C].32nd International Electric Propulsion Conference,Wiesbaden，Germany：IEPC，2011: 15.

[3] MATYASH K,SCHNEIDER R,MUTZKE A,et al.Kinetic simulations of SPT and HEMP thrusters including the near-field plume region[J].IEEE Transactions on Plasma Science,2009,38(9):2274-2280.

[4] BRANDT T，TROTTENBERG T，GROLL R，et al.Particle-in-cell simulation of the plasma properties and ion acceleration of a down-scaled HEMP-Thruster[C].30th International Spacecraft Propulsion Conference,Cologne,Germany：ISPC，2014.

[5] BOEUF J P，GARRIGUES L.Low frequency oscilla-tions in a stationary plasma thruster[J].Journal of Applied Physics，1998，84(7):3541-3554.

[6] 胡俊鋒，劉輝，李建置，等.會切磁場推力器低頻振蕩特性[J].中國空間科學技術,2016，36(1)：26-34.

HU J F，LIU H，LI J Z，et al.Low frequency oscillation characteristics in a cusped field thruster[J].Chinese Space Science and Technology，2016，36(1):26-34(in Chinese).

[7] MOROZOV A I，SAVEL′EV V V.One-dimensional hybrid model of a stationary plasma thruster[J].Plasma Physics Reports，2000，26(10):875-880.

[8] MOROZOV A I，SAVEL′EV V V.One-dimensional hydrody-namic model of the atom and ion dynamics in a stationary plasma thruster[J].Plasma Physics Reports，2000，26(3):219-224.

[9] BARRAL S，AHEDO E.Theoretical study of the breathing mode in Hall thrusters[C]//42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference，Scacramento,California,July 9-12,2006.

[10] BARRAL S，AHEDO E.Low-frequency model of breathing oscillations in Hall discharges[J].Physical Review E,2009，79(4):046401.

[11] HARA K，BOYD I D，KOLOBOV V I.One-dimens-ional hybrid-direct kinetic simulation of the discharge plasma in a Hall thruster[J].Physics of Plasmas，2012，19(11):530-541.

[12] HARA K，SEKERAK M J，BOYD I D，et al.Mode transition of a Hall thruster discharge plasma[J].Journal of Applied Physics，2014，115(20):203304-2033010.

[13] LIU H，WU H，ZHAO Y，et al.Study of the electric field formation in a multi-cusped magnetic field[J].Physics of Plasmas，2014，21(9):214.

[14] MOROZOV A I，SAVEL′EV V V.One-dimensional hydrodynamic model of the atom and ion dynamics in a stationary plasma thruster[J].Plasma Physics Reports，2000，26(3): 219-224.

[15] MOROZOV A I，SAVEL′EV V V.One-dimensional hybrid model of a stationary plasma thruster[J].Plasma Physics Reports，2000，26(10): 875-880.

[16] FIFE J M.Hybrid-PIC modeling and electrostatic probe survey of Hall thrusters[D].Massachusetts: Massachusetts Institute of Technology，1998:189-195.

[17] BARRAL S，AHEDO E.Theoretical study of the breathing mode in Hall thrusters[C]//42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit,Scacramento,California,July 9-12,2006: 5172.

(編輯：車曉玲)

Simulationresearchonoscillationcharacteristicinhighefficiencymulti-plasmathruster

NIU Xiang，WU Huan，LIU Hui*，YU Daren

InstituteofAdvancedPower，HarbinInstituteofTechnology，Harbin150001，China

Now the main numerical simulation method of high efficiency multi plasma thruster (HEMPT) is particle in cell.But this method puts higher demands on computer.Therefore，this method aims at reappearing the oscillation phenomenon in high efficiency multi-plasma thruster by building up one dimension fluid model and studying the plasma parameters distribution and the relationship between the current and plasma parameters.A one dimension fluid model was built up along two distinguishing typical electron paths in thrusters by Simulink.The simulation results show that the plasma parameters oscillation period is consistent with that of the current.The main ionization zone is between the two cusped fields.The main area where potential drops is near the channel exit.

high efficiency multi-plasma thruster；one dimension fluid model；oscillation；current；electron conductivity；potential distribution

http://zgkj.cast.cn

10.16708/j.cnki.1000-758X.2017.0072

V430

A

2017-05-04；

2017-08-03；錄用日期2017-09-12；< class="emphasis_bold">網絡出版時間

時間：2017-09-24 16:01:06

http://kns.cnki.net/kcms/detail/11.1859.V.20170924.1601.008.html

國家自然科學基金(51421063，11505041，11275055)

牛翔(1995-)，男，碩士研究生，15776867067@163.com，研究方向為振蕩與噪聲分析

*通訊作者：劉輝(1981-)，男，副教授，huiliu@hit.edu.cn，研究方向為空間電推進

牛翔，伍環，劉輝，等.多級會切場推力器振蕩特性模擬研究[J].中國空間科學技術，2017，37(5)：75-80.NIUX，WUH，LIUH，etal.Simulationresearchonoscillationcharacteristicinhighefficiencymulti-plasmathruster[J].ChineseSpaceScienceandTechnology，2017，37(5)：75-80(inChinese).