，，，,*

1.北京航空航天大學 宇航學院，北京 100191 2.北京控制工程研究所，北京 100094 3.空間物理重點實驗室，北京 100076

同軸微陰極電弧推力器的粒子網(wǎng)格法數(shù)值模擬

熊子昌1，耿金越2，陳軒3，王海興1,*

1.北京航空航天大學 宇航學院，北京 100191 2.北京控制工程研究所，北京 100094 3.空間物理重點實驗室，北京 100076

采用粒子網(wǎng)格(Particle-In-Cell,PIC)方法對同軸微陰極電弧推力器(μCAT)工作過程進行了模擬研究，并應用自相似方法對模型進行簡化，獲得了推力器羽流區(qū)的電子數(shù)密度分布、離子數(shù)密度分布、電勢分布及離子軸向平均速度，通過改變磁感應強度和位形分析磁場對推力器內等離子體運動特性及推力器性能的影響。計算結果表明，電子被外加磁場捕獲約束在磁力線附近，低速離子與高速電子形成的雙極擴散電場加速離子噴出；在相同流量情況下，磁感應強度0.02 T時，離子返流嚴重，磁感應強度0.05～0.30 T時，磁感應強度變化對速度影響較小；磁場位形對離子運動和推力器性能有較大影響，磁力線與軸線夾角較小時離子速度下降明顯，夾角較大時離子返流嚴重。

同軸微陰極電弧推力器；粒子網(wǎng)格法；自相似法；磁感應強度；磁場位形

微陰極電弧推力器是美國喬治華盛頓大學近年來研制的一種新型推力器，主要有同軸和環(huán)形兩種結構[1]，同軸微陰極電弧推力器結構示意如圖1所示。其主要工作原理是電弧燒蝕陰極材料產生較高電離度的等離子體，在電場和磁場的共同作用下高速噴出而產生推力，外加磁場能夠提高等離子體軸向速度,減小羽流擴散角,提高推力器比沖。該推力器結構簡單，不需要復雜的推進劑供給系統(tǒng)，比沖較高，因而受到廣泛關注，在微小衛(wèi)星推進系統(tǒng)應用方面有廣泛的應用前景[2-5]。

圖1 同軸微陰極電弧推力器結構示意Fig.1 Schematic drawing of co-axial micro- cathode arc thruster

從目前已有的研究來看，推力器結構、磁感應強度及分布、脈沖放電能量的大小和持續(xù)時間、陰極材料的選取等都會對推力器的推力、比沖等性能參數(shù)產生較大影響[6-9]。但對該影響的物理機制和物理過程的深入研究尚無報道[10]。為了深入了解推力器結構與工作參數(shù)的改變對推力器性能影響，為推力器的優(yōu)化設計提供參考，縮短推力器研制周期，有必要對微陰極電弧推力器工作過程進行系統(tǒng)的數(shù)值模擬研究。

本文采用粒子網(wǎng)格(Particle-In-Cell,PIC)方法對同軸微陰極電弧推力器進行模擬，應用自相似方法對模型進行簡化，獲得電子數(shù)密度分布、離子數(shù)密度分布、電勢分布、離子軸向平均速度分布等結果，對微陰極電弧推力器的工作過程進行分析，估計推力器比沖；并且對微陰極電弧推力器影響最大的因素——磁場進行進一步分析，研究了磁感應強度和磁場位形對推力器性能的影響。

1 模擬方法

1.1 PIC法

PIC法[11]：通過向計算域內放置一些宏粒子來模擬真實的等離子體粒子的運動，每個宏粒子可以看作一起運動的等離子體團，其荷質比與真實粒子相等。

PIC方法的主要流程為：1)將網(wǎng)格內粒子的電荷分配到網(wǎng)格節(jié)點上獲得電荷分布；2)在網(wǎng)格節(jié)點上求解泊松方程計算電場分布；3)在電磁場分布的基礎上，按照離子所受電場和磁場作用力來求解宏粒子的運動狀況。

因為微陰極電弧推力器中等離子體電離度很高，可以看作完全電離等離子體，所以，本文不考慮激發(fā)、電離等碰撞反應。微陰極電弧推力器中的外加磁感應強度遠遠大于放電所產生的自感磁場磁感應強度，所以認為推力器中磁場不變。在每個計算周期求解運動方程(式(1)(2))與泊松方程(式(3))以推動粒子運動和求解電場。

式中：m為粒子質量；v為粒子速度；E,B分別為電場強度和磁感應強度;x為粒子位置；φ為空間電勢；ρ為空間電荷密度；ε為介電常數(shù)。

為了使得數(shù)值算法精確可靠，時間步長的選取依賴于等離子體德拜長度，空間步長的選取依賴于等離子體頻率[7]。

式中：Δx為空間步長；Δt為時間步長；λD為等離子體德拜長度；ωpe為電子振蕩頻率；ωce為電子回旋頻率。

1.2 自相似方法

在同軸微陰極電弧推力器二維模擬中，德拜長度為10 μm量級，推力器羽流特征長度為厘米量級，對于全尺寸的二維PIC模擬，模型計算域需要劃分的網(wǎng)格節(jié)點數(shù)為百萬量級，為保證模擬精度，所需要模擬的粒子數(shù)將為上億個，這樣龐大的計算量對于計算機的存儲和運算速度是一個巨大的挑戰(zhàn)。因此，模型的簡化是很有必要的。

傳統(tǒng)的模型簡化方法有減少重粒子質量法和增大介電常數(shù)法，這兩種方法可以有效地減少收斂所需要的循環(huán)步數(shù)和求解所需要的網(wǎng)格數(shù)，但是會破壞等離子體的物理過程，使結果有較大失真。對此，Taccogna等提出了自相似法，并驗證了其可行性[12-14]。自相似法可以通過保證一些重要的物理、性能參數(shù)不變，以保證縮比后的模型物理過程相似，很大程度上保留了縮放后物理特性，故本文選擇自相似法對同軸微陰極電弧推力器模型進行簡化。但由于霍爾推力器的PIC模擬需要考慮電離過程，所以努森數(shù)為Taccogna自相似方法中的重要縮放參數(shù)。而本文中，因為等離子體電離度很高，假設為完全電離等離子體，不考慮電離過程，所以并不需要維持努森數(shù)不變。

(1)離子速度不變

(2)電流密度不變

故離子數(shù)密度也不變，推力器的粒子電流I=jA，A為面積。所以I～L2。

(3)帶電粒子回旋半徑和長度的比值不變

所以有B～1/L。

表1 各物理量對應的相似準則數(shù)Table 1 Physical quantities and their self-similarity factors

1.3 模型建立

微陰極電弧推力器放電過程中，外加磁場磁感應強度大小遠大于放電過程所產生的自感磁場磁感應強度大小，因此認為放電過程中磁感應強度不變。

(1)磁場計算

磁場計算在商業(yè)軟件COMSOL中的“磁場”模塊下進行，根據(jù)推力器線圈的位置、結構和材料對推力器通道內及羽流區(qū)磁場進行有限元計算，得到磁場分布如圖2所示。

圖2 同軸微陰極電弧推力器磁場分布Fig.2 Distribution of magnetic field of the co-axial micro-cathode arc thruster

(2)計算域選取

為了盡可能保留推力器的羽流變化特性，本文選取0.06 m×0.03 m的矩形計算域，下游出口邊界為零電勢邊界，軸線為軸對稱邊界，其余邊界為開放邊界。并在此基礎上對模型采用自相似方法進行簡化，縮放比為0.05。計算域網(wǎng)格數(shù)為300×150，計算時間步長為1×10-12s。

平均自由程為離子在兩次碰撞之間平均運動的距離，平均自由程越大，則粒子間發(fā)生碰撞的概率越小[15]。假設種類1的粒子與作為背景的種類2粒子發(fā)生碰撞，則兩種粒子碰撞的平均自由程為：

式中：n2為背景粒子的數(shù)密度；Q為粒子間的碰撞截面。

推力器電子溫度為幾電子伏特，電子-離子碰撞截面為10-18m2量級，推力器等離子體數(shù)密度為1018m-3量級，故電子-離子碰撞自由程為1 m量級，大于推力器特征長度毫米量級，因而本文中忽略電子-離子碰撞導致的復合。

由于推力器羽流電離度很高，羽流等離子體99%被電離[16]，近似完全電離，所以模型中，假設等離子體只由Ti二價離子與電子組成。額定工況下，磁感應強度為0.3 T，Ti離子的質量流量為3.8×10-6mg/120 μs，陰陽極之間的電勢差為50 V，文獻[17]表明，即使磁場會影響推力器的最終速度，但是它不會改變Ti離子在電極間的初始速度，這個速度為2×104m/s，本文將此速度作為模擬中離子的入口速度。

2 數(shù)值模擬結果

2.1模型驗證

由于同軸微陰極電弧推力器尚無詳細的試驗數(shù)據(jù)報道，本文首先對不同磁場下的環(huán)形微陰極電弧推力器進行模擬，得到不同磁場及對應流量下離子平均速度，并與文獻報道的試驗測量結果[18]進行比較，如圖3所示。模擬結果與文獻報道試驗結果符合較好，驗證了模型的正確性。

圖3 離子平均速度對比Fig.3 Comparison of average ion velocity

2.2 磁感應強度為0.3 T額定工況模擬結果

圖4和圖5分別為電子數(shù)密度分布和離子數(shù)密度分布。圖中左下角黃色區(qū)域為推力器位置。電子質量很小，在磁場中做回旋運動的Larmor半徑很小，因此被外加磁場捕獲約束在磁力線附近，形成電子聚集區(qū)。而質量大、速度慢的離子，則由于雙極擴散電場的作用，聚集于電子聚集區(qū)附近并加速噴出。圖5可以觀察到明顯的沿電子運動區(qū)域的離子主流區(qū)。

圖4 電子數(shù)密度分布Fig.4 Distribution of electron number density

圖5 離子數(shù)密度分布Fig.5 Distribution of ion number density

圖6為電勢分布。由圖6可知，磁力線捕獲約束電子，并在此區(qū)域形成了狹長的低電勢通道，電子和大部分離子均沿此通道運動。在圖6電勢分布中，沿軸線方向存在高電勢區(qū)，這是因為陰極表面高電勢與電子低電勢區(qū)存在徑向電場，使得陰極表面發(fā)出的離子在電子低電勢區(qū)附近震蕩，使得離子數(shù)密度在軸線較高，而此處并沒有電子與之中和，從而產生局部高電勢區(qū)。

圖6 電勢分布Fig.6 Distribution of electric potential

圖7給出了離子軸向平均速度沿軸線變化。離子軸向平均速度通過統(tǒng)計各軸向位置所在徑向平面內離子的軸向平均速度而得到。由圖7中0.3 T條件下離子軸向平均速度可知，離子從陰極表面噴出，經(jīng)過陰極表面附近的高電勢區(qū)域速度急劇下降，之后進入低電勢通道后速度逐漸增加至25.3 km/s左右,對應的比沖約為25 800 m/s。

圖7 不同磁感應強度離子軸向平均速度Fig.7 Axial average ion velocity in different magnetic flux density

由以上仿真結果可知，磁場在微陰極電弧推力器中有至關重要的作用，將直接影響等離子體的運動過程及推力器性能。下面將分別改變外加磁場磁感應強度大小及位形，研究其影響規(guī)律。

2.3 磁感應強度大小影響

為了研究磁感應強度大小對推力器內等離子體運動過程及推力器性能的影響，本文還采用圖2相同的磁場位形，僅同比減小磁感應強度，使陰極斑點附近磁感應強度分別為0.15 T、0.05 T和0.02 T，并分別對其進行模擬。結果如圖7～圖9所示。

圖8和圖9分別為電子數(shù)密度分布和離子數(shù)密度分布。電子的分布情況將通過雙極擴散電場影響離子分布，由圖8和圖9可知，隨著磁感應強度的減弱，磁場對電子的束縛能力逐漸減弱，導致電子聚集區(qū)的“寬度”逐漸增加；在磁感應強度為0.02 T時電子聚集區(qū)寬度進一步增加，甚至部分電子開始脫離磁力線的束縛，同時離子也出現(xiàn)了明顯的返流現(xiàn)象。

圖7為不同磁感應強度下離子軸向平均速度沿軸向變化。由圖7可知，當磁感應強度為0.05～0.30 T時，隨著磁感應強度變化，離子軸向平均速度變化趨勢基本一致，且數(shù)值相當。這與圖4中平均速度隨磁感應強度增大而增大不同，主要是因為圖4中，不同磁感應強度采用離子流量不同，造成了各磁場條件下離子速度的差異。而本文中為了分析磁感應強度單一條件對等離子體流動的影響，采用了相同流量。磁感應強度為0.02 T時離子軸向平均速度與0.05～0.30 T有較大差別，主要是部分離子出現(xiàn)返流引起的。

圖8 不同磁感應強度下的電子數(shù)密度分布Fig.8 Distribution of electron number density in different magnetic flux density

總的來說，當磁感應強度為0.05～0.30 T時，隨著磁感應強度變化，離子軸向平均速度變化趨勢基本一致，隨著磁感應強度的增大而小幅增大；當磁感應強度為0.02 T時，離子數(shù)密度降低嚴重，甚至出現(xiàn)返流的現(xiàn)象。

圖9 不同磁感應強度下的離子數(shù)密度分布Fig.9 Distribution of ion number density in different magnetic flux density

2.4 磁場位形影響

除了磁感應強度大小，磁場位形也對等離子體運動及推力器性能有較大影響。

通過改變電磁線圈與推力器陰陽極之間的相對位置以改變磁場位形，即分別沿z軸正負方向移動3 mm，使得穿過陰極端面與絕緣介質端面交界點的磁力線發(fā)生變化，磁力線與z軸的夾角由額定工況下的11°分別變化為22°、2°，磁場與推力器相對位置示意如圖10所示，圖中藍色、灰色、紅色部分為推力器移動后的位置，帶箭頭的曲線表示磁力線。分別模擬這兩種磁場條件下額定工況的結果，模擬中改變計算域中磁場位置，而推力器在計算域中的位置不變。

圖10 穿過陰極端面與絕緣介質端面交界點的磁力線與軸線夾角示意Fig.10 Schematic drawing of 2 degrees between axis and the magnetic field line passing through the junction point of cathode end face and insulator end face

圖11為磁力線與軸線夾角為22°和2°情況下的電子數(shù)密度分布。由圖11可知，由于電子被磁化，電子運動軌跡與穿過陰極端面與絕緣介質端面交界點的磁力線一致。磁力線與軸線夾角度數(shù)越大表明磁力線越靠近z軸，導致電子運動也越靠近z軸。

圖12為磁力線與軸線夾角為22°和2°情況下的離子數(shù)密度分布。由圖12可知，離子主流區(qū)隨電子運動軌跡改變而改變，即隨著磁力線的改變而改變。當夾角很大時(22°)，磁力線與陽極很靠近，有一部分離子會與陽極發(fā)生碰撞，導致碰撞后的離子很有可能形成返流。當夾角很小時(2°)，磁力線與陽極距離很遠，不但不會形成返流，還會減輕陽極附近的離子震蕩幅度，形成穩(wěn)定束流。但是因為磁力線位置遠離z軸，造成離子主流區(qū)出口速度方向與z軸夾角變大。

圖11 不同磁場位形電子數(shù)密度分布Fig.11 Distribution of electron number density in different magnetic field topology

圖12 不同磁場位形離子數(shù)密度分布Fig.12 Distribution of ion number density in different magnetic field topology

圖13為不同磁場位形情況下離子軸向平均速度沿軸向變化。由圖13可知，當夾角很大時，雖然出口速度有明顯提升，但是會造成離子返流現(xiàn)象。當夾角為2°時，離子主流區(qū)出口速度與z軸夾角變大，導致離子軸向平均速度減小明顯，對推力貢獻不利。而這3種工況之間相對位置僅僅相差3 mm。

圖13 不同磁場位形離子軸向平均速度Fig.13 Axial average ion velocity in different magnetic field topology

3 結束語

本文采用2D軸對稱粒子網(wǎng)格法對同軸微陰極電弧推力器等離子體運動過程進行數(shù)值模擬，為了加速算法的收斂時間，將一種可以縮小發(fā)動機幾何尺寸，同時保證模擬時重要的參數(shù)不變的自相似法應用其中，以保證縮放前后結果一致。

計算結果發(fā)現(xiàn)：

1)推力器產生的電子被外加磁場捕獲約束在磁力線附近，低速的離子由于雙極擴散電場的作用，聚集于電子聚集區(qū)附近并加速噴出。

2)當磁感應強度為0.05～0.30 T時，隨著磁感應強度的減弱，離子軸向出口速度逐漸減小，但總體趨勢基本一致，當磁感應強度為0.02 T時，磁場對電子的束縛減弱明顯，離子返流嚴重，出口流量降低明顯。

3)磁場位形對推力器性能影響明顯，穿過陰極端面與絕緣介質端面交界點的磁力線與軸線夾角較小導致軸向速度分量下降明顯，而夾角較大則會導致返流嚴重。

為了進一步了解同軸微陰極電弧推力器系統(tǒng)的工作過程，今后的研究需要建立更加精細的推力器模型例如考慮壁面效應、陽極效應等。

References)

[1] KEIDAR M,ZHUANG T,SHASHURIN A,et al.Electric propulsion for small satellites[J].Plasma Physics & Controlled Fusion,2015,57(1):014005.

[2] ZHUANG T,SHASHURIN A,DENZ T,et al.Micro-vacuum arc thruster with extended lifetime[C].45th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit,Denver,August 2-5,2009.

[3] ZHUANG T,SHASHURIN A,HAQUE S,et al.Performance characterization of the micro-cathode arc thruster and propulsion system for space applications[C].46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit,Nashville,July 25-28,2010.

[4] HAQUE S,KEIDAR M,LEE T.Low-thrust orbital maneuver analysis for cubesat spacecraft with a micro-cathode arc thruster subsystem[C].33rd International Electric Propulsion Conference,The George Washington University,Washington D.C.,October 6-10,2013.

[5] SHASHURIN A,KEIDAR M.Comparative analysis of micro-cathode arc thruster performance[C].33rd International Electric Propulsion Conference,The George Washington University,Washington,D.C.,October 6-10,2013.

[6] ZHUANG T,SHASHURIN A,KEIDAR M.Micro cathode thruster (μCT) plume characterization[J].IEEE Transactions on Plasma Science,2011,39(11): 2936-2937.

[7] VAIL P,PANCOTTI A.Performance characterization of micro-cathode arc thruster (μCAT)[C].47th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit,San Diego,July 31-August 3,2011.

[8] HURLEY S,KEIDAR M.Linear actuated micro-cathode arc thruster system [C].52nd AIAA/SAE/ASEE Joint Propulsion Conference,Salt Lake City,July 25-27,2016.

[9] BRIEDA L,ZHUANG T,KEIDAR M.Towards near plume modeling of a micro cathode arc thruster[C].49th AIAA/ASME/SAE/ASEE Joint Propulsion Conference,San Jose,July 14-17,2013.

[10] ZHUANG T,SHASHURIN A.Circular periodic motion of plasma produced by a small-scale vacuum arc[J].Plasma Sources Science and Technology,2011,20:015009.

[11] VERBONCOEUR J P.Particle simulation of plasmas: review and advances[J].Plasma Physics & Controlled Fusion,2005,47(47):A231.

[12] TACCOGNA F,LONGO S,CAPITELLI M,et al.Plasma flow in a Hall thruster[J].Physics of Plasmas,2005,12(4):78.

[13] TACCOGNA F,LONGO S,CAPITELLI M,et al.Self-similarity in Hall plasma discharges: applications to particle models[J].Physics of Plasmas,2005,12(5):78.

[14] TACCOGNA F,SCHNEIDER R,LONGO S,et al.Kinetic simulations of a plasma thruster[J].Plasma Sources Science and Technology,2008,17(2):1-16.

[15] SZABO J.Fully kinetic numerical modeling of a plasma thruster [D].Cambridge: Massachusetts Institute of Technology,2001:92-106.

[16] KEIDAR M,HAQUE S,ZHUANG T,et al.Micro-cathode arc thruster for PhoneSat propulsion [C].27th Annual AIAA/USU Conference on Small Satellites,Logan,2013: 1-9.

[17] ZHUANG T,SHASHURIN A,BEILIS I,et al.Ion velocities in a micro-cathode arc thruster[J].Physics of Plasmas,2012,19(6):661.

[18] ZHUANG T.Micro-cathode arc thruster system for cube satellite[D].Washington,D.C.: George Washington University,2012:83-85.

(編輯：高珍)

Particle-in-cellsimulationofaco-axialmicro-cathodearcthruster

XIONG Zichang1，GENG Jinyue2，CHEN Xuan3，WANG Haixing1,*

1.SchoolofAstronautics,BeihangUniversity，Beijing100191,China2.BeijingInstituteofControlEngineering，Beijing100094，China3.ScienceandTechnologyonSpacePhysicsLaboratory，Beijing100076，China

A modeling study was performed to investigate the working process of a co-axial micro-cathode arc thruster with particle-in-cell (PIC) method,and the self-similar method was applied to simplify the simulation model.The distribution of electron number density,ion number density,electric potential and axial average ion velocity in the thruster plume region was obtained.The influences of magnetic field on plasma motion characteristics and thruster performance were studied by changing magnetic flux density and magnetic field topology.The results show that electrons are constrained by the applied magnetic field,cycling around magnetic field lines,and that ions are accelerated by the ambipolar diffusion electric field produced by slower ions and faster electrons.With the same mass flow rate,the ion reflux is severe when the magnetic flux density is 0.02 T,while the ion velocity is less affected by the magnetic flux density when the magnetic flux density is 0.05～0.30 T.Magnetic field topology has great influence on the ion motion and thruster performance.The ion velocity decreases obviously when the angle between the axis and magnetic field lines is small,while the ion reflux is severe when the angle is large.

co-axial micro-cathode arc thruster; particle-in-cell; self-similarity; magnetic flux density; magnetic field topology

http://zgkj.cast.cn

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

V439

A

2017-05-04；

2017-07-24；錄用日期2017-09-12；< class="emphasis_bold">網(wǎng)絡出版時間

時間：2017-09-25 14:07:15

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

國家自然科學基金(11275021，11575019,11702021)；民用航天項目(混合模式推進系統(tǒng)優(yōu)化技術)

熊子昌(1995-)，男，碩士研究生，xiongzichang@buaa.edu.cn，研究方向為空間電推進

*通訊作者：王海興(1969-)，男，教授，whx@buaa.edu.cn，研究方向為空間電推進

熊子昌，耿金越，陳軒，等.同軸微陰極電弧推力器的粒子網(wǎng)格法數(shù)值模擬[J].中國空間科學技術,2017,37(5)：81-88.XIONGZC,GENGJY,CHENX,etal.Particle-in-cellsimulationofaco-axialmicro-cathodearcthruster[J].ChineseSpaceScienceandTechnology,2017,37(5)：81-88(inChinese).