劉慶海，張紹華，鄭慧奇，王俊峰，任瓊英，趙 華，易 忠，黃建國(guó)

（北京衛(wèi)星環(huán)境工程研究所，北京 100094）

0 引言

高比沖、高效率、長(zhǎng)壽命的電推進(jìn)器在衛(wèi)星姿態(tài)控制、軌道保持及深空探測(cè)等領(lǐng)域具有性能優(yōu)勢(shì)和競(jìng)爭(zhēng)力[1-2]。其中螺旋波等離子體推進(jìn)器尤為引人關(guān)注。美國(guó)、法國(guó)、澳大利亞、荷蘭、意大利、日本等國(guó)家的相關(guān)研究機(jī)構(gòu)正在對(duì)螺旋波等離子體的產(chǎn)生機(jī)制、電雙層加速機(jī)制等開展理論研究，并通過實(shí)驗(yàn)手段診斷螺旋波等離子體的特性參數(shù)[2-6]。北京衛(wèi)星環(huán)境工程研究所正致力于螺旋波等離子體推進(jìn)器的研制，已取得了初步實(shí)驗(yàn)成果[7-10]。

本文主要針對(duì)螺旋波等離子體推進(jìn)器系統(tǒng)對(duì)磁場(chǎng)的要求，采取數(shù)值仿真的手段開展推進(jìn)器磁場(chǎng)位形的研究，建立三維磁場(chǎng)仿真模型。

1 螺旋波等離子體推進(jìn)器基本原理

圖1所示為螺旋波等離子體推進(jìn)器基本原理。圖中白色虛線是磁力線，彩色點(diǎn)線是由推進(jìn)工質(zhì)電離的離子之運(yùn)動(dòng)軌跡，右端噴口附近白色霧狀區(qū)域顯示的是電雙層加速區(qū)。推進(jìn)工質(zhì)從左端送入帶有螺旋波天線的放電管被電離，電離后的離子順著磁力線作螺旋運(yùn)動(dòng)，而后進(jìn)入出口處的磁場(chǎng)中，由電雙層加速為高速離子束流，從右端噴出而產(chǎn)生反推力。電雙層由正電荷層和負(fù)電荷層組成，盡管整體顯電中性，但在電雙層中存在電勢(shì)梯度，因而可以用于離子加速[11-13]。

圖1 螺旋波等離子體推進(jìn)器的基本原理示意圖 Fig.1 The working principle of the HPT

在螺旋波等離子體推進(jìn)器中，推進(jìn)工質(zhì)的電離、加速、噴出過程以及電羽流的污染控制等都與磁場(chǎng)密切相關(guān)：1）螺旋波是在徑向被約束的磁化等離子體中傳播的右旋極化電磁波，在螺旋波等離子體產(chǎn)生機(jī)制中，磁場(chǎng)束縛自由電子的徑向運(yùn)動(dòng)；2）通道內(nèi)強(qiáng)磁場(chǎng)的約束作用可減少高溫等離子體與推進(jìn)器壁面的接觸，削弱離子的壁面濺射侵蝕，延長(zhǎng)推進(jìn)器壽命；3）在噴口附近利用喇叭口形狀的磁場(chǎng)結(jié)構(gòu)產(chǎn)生穩(wěn)定的電雙層結(jié)構(gòu)，是加速離子和產(chǎn)生超高速噴流的關(guān)鍵部位，喇叭口狀的磁場(chǎng)結(jié)構(gòu)也稱為磁噴嘴。De Laval 化學(xué)噴嘴把燃燒室內(nèi)高溫、高壓、低速的化學(xué)氣體變成低溫、低壓、高速的氣流噴出，產(chǎn)生反推力[14-15]。因此，De Laval 噴嘴的尺寸設(shè)計(jì)必須非常精確，才能有效完成等離子體能量和動(dòng)量的轉(zhuǎn)換，提高電推進(jìn)器效率?；诨瘜W(xué)噴嘴尺寸設(shè)計(jì)要求，應(yīng)進(jìn)行磁噴嘴結(jié)構(gòu)精確設(shè)計(jì)。磁場(chǎng)位形的設(shè)計(jì)對(duì)于控制螺旋波等離子體的產(chǎn)生、通道內(nèi)螺旋波的傳播與吸收、磁噴嘴出口處離子的加速和逃逸，以及提高電推進(jìn)器的性能至關(guān)重要。

電推進(jìn)器所需的磁場(chǎng)可用銅線圈、永磁體或者超導(dǎo)磁體技術(shù)來實(shí)現(xiàn)[6]。銅線圈方案可通過電流調(diào)節(jié)來控制磁場(chǎng)，但發(fā)熱嚴(yán)重，需要配置冷卻系統(tǒng)；永磁體方案不需要電力和冷卻系統(tǒng)，但磁場(chǎng)位形設(shè)計(jì)難度大，磁場(chǎng)強(qiáng)度不易調(diào)節(jié)；超導(dǎo)線圈方案不發(fā)熱、體積小、重量輕，產(chǎn)生的磁場(chǎng)強(qiáng)度高，但需要低溫設(shè)施，技術(shù)實(shí)現(xiàn)難度較大。綜合考慮后采用永磁體產(chǎn)生螺旋波等離子體推進(jìn)器所需的磁場(chǎng)。

2 磁場(chǎng)仿真模型

2.1 等效磁荷原理

對(duì)于簡(jiǎn)單形狀的永磁體，可直接給出磁場(chǎng)的解析表達(dá)式；但對(duì)于復(fù)雜形狀的永磁體，則需要借助于數(shù)值計(jì)算來得到磁場(chǎng)位形分布。永磁體磁場(chǎng)仿真計(jì)算的主要方法有等效磁荷法、等效電流法和有限元法[16-17]。其中，等效磁荷法所用的等效磁荷模型認(rèn)為磁化體內(nèi)部存在密度為ρm的體磁荷，磁化邊界處存在密度為ρms的面磁荷；等效電流法所用的等效電流模型認(rèn)為存在體電流和面電流。

本研究中采用等效磁荷法。等效磁荷模型給出永磁體周圍空間磁場(chǎng)強(qiáng)度表達(dá)式為

式中：R是永磁體內(nèi)源點(diǎn)到磁場(chǎng)求解點(diǎn)的矢徑；R是源點(diǎn)到求解點(diǎn)的距離；V是永磁體積分區(qū)域；A是永磁體積分區(qū)域表面邊界。體磁荷密度ρm=-μ0?·M，其中μ0是真空磁導(dǎo)率，M是磁化強(qiáng)度。在均勻磁體內(nèi)部，由于?·M=0，所以體磁荷密度ρm=0。在磁體表面邊界，由于磁化強(qiáng)度M不連續(xù)，存在面磁荷密度ρms=μ0en·M，en是磁體邊界的外法向單位矢量。

2.2 模型驗(yàn)證

用Fortran 程序開發(fā)出三維等效磁荷模型軟件，并進(jìn)行驗(yàn)證。在三維仿真計(jì)算中，只有表面磁荷對(duì)磁場(chǎng)有貢獻(xiàn)，因此，利用式(1)第二項(xiàng)求解，磁場(chǎng)在空間某一點(diǎn)的值是所有表面離散網(wǎng)格點(diǎn)的積分值[17]。對(duì)三維空間區(qū)域進(jìn)行網(wǎng)格劃分，求解每個(gè)網(wǎng)格點(diǎn)的磁場(chǎng)值，就得到三維空間的磁場(chǎng)位形分布。

實(shí)驗(yàn)中使用的永磁體是銣鐵硼永磁體，它的質(zhì)量密度是8.0 g/cm3，面磁荷密度是1.0 Wb/cm2。為了驗(yàn)證模型的正確性，計(jì)算了直徑40 mm、高40 mm 的圓柱形永磁體周圍的磁力線及磁場(chǎng)強(qiáng)度分布。圓柱體充磁方向沿著z軸方向。圖2是圓柱形永磁體軸心截面的磁矢（黑色箭頭）、磁場(chǎng)分布及磁力線（白色帶箭頭線）?？梢钥闯觯艌?chǎng)強(qiáng)度在永磁體附近最強(qiáng)（約1 T 左右），在離磁體50 mm 之外衰減到0.1 T 以下，這是永磁體磁場(chǎng)強(qiáng)度分布的典型特征。

圖2 圓柱形永磁體軸心截面（x-z）的磁場(chǎng)分布 Fig.2 The distribution of magnetic field in x-z section

3 磁場(chǎng)仿真結(jié)果

首先，針對(duì)圓柱和圓環(huán)永磁體，構(gòu)造螺旋波等離子體推進(jìn)器所需的磁場(chǎng)位形，各磁體同軸分布，且磁體相互之間的軸向距離可調(diào)（見圖3）。磁體的具體尺寸和質(zhì)量參見表1。

圖3 磁體形狀及位置示意圖 Fig.3 The permanent magnet configuration of the HPT system

表1 磁體尺寸及質(zhì)量 Table1 The sizes and weights of the permanent magnets

圖4是圓柱磁體B1 與圓環(huán)磁體B2 之間，以及4 個(gè)圓環(huán)磁體（B2～B5）之間的軸向間距均為10 mm 時(shí)，軸截面的磁矢（白色箭頭）、磁場(chǎng)分布和磁力線（紅色實(shí)線）。圖中的黃色方框代表放電管所在位置?？梢钥闯觯诖朋w內(nèi)側(cè)和外側(cè)與磁體相鄰的地方以及放電管的右端有磁場(chǎng)較弱（小于0.01 T）的區(qū)域，這是不同方向磁場(chǎng)交匯的結(jié)果。從追蹤出的幾條典型磁力線可以看出，磁場(chǎng)在放電管內(nèi)部磁場(chǎng)基本是軸向均勻的，在放電管右端出口處有發(fā)散的趨勢(shì)。

圖4 磁場(chǎng)疊加磁力線的分布 Fig.4 The contour of magnetic field and magnetic force line

為更清楚地展示電推進(jìn)器工作區(qū)域的磁場(chǎng)強(qiáng)度分布，圖5給出了磁場(chǎng)強(qiáng)度沿著中心軸的分布及局部放大圖，其中黃色線段示意的是圓環(huán)磁體內(nèi)部區(qū)域，即放電管所在位置。從圖5(b)可以看出，圓柱形永磁體內(nèi)部的磁場(chǎng)最強(qiáng)，達(dá)到約1 T 量級(jí)。圓環(huán)磁體區(qū)域的磁場(chǎng)強(qiáng)度，在最左端是0.2 T，到中心區(qū)域下降至約0.03 T，在放電管右端出口處接近于0；遠(yuǎn)離放電管出口50 mm 處上升至約0.02 T，之后往更遠(yuǎn)處磁場(chǎng)強(qiáng)度再次緩慢衰減至0。

圖5 磁場(chǎng)沿著z 軸的分布及局部放大圖 Fig.5 The profiles of magnetic field along z-axis and the local zoom-in plot

對(duì)于螺旋波等離子體推進(jìn)器，若以13.56 MHz的射頻源產(chǎn)生1018m-3的高密度等離子體，則需要磁場(chǎng)強(qiáng)度約為0.02T[3-5]。圖5(b)顯示永磁體組合產(chǎn)生的磁場(chǎng)在放電管內(nèi)均勻度低于20%，強(qiáng)度也基本滿足產(chǎn)生螺旋波等離子體所需條件。但電推進(jìn)器噴口附近存在的弱磁場(chǎng)區(qū)導(dǎo)致磁場(chǎng)不連續(xù)，可能對(duì)等離子體羽流噴射造成一定的不利影響[18]，需要深入研究分析。

本研究中的三維磁場(chǎng)仿真模型是專門針對(duì)螺旋波等離子體推進(jìn)器開發(fā)設(shè)計(jì)的。圖6給出優(yōu)化改進(jìn)磁體組合之后，再仿真得到的磁場(chǎng)位形。由于永磁體組合是軸對(duì)稱分布的，產(chǎn)生的磁場(chǎng)也是三維空間軸對(duì)稱分布的，所以只給出軸截面的磁場(chǎng)強(qiáng)度和磁力線的分布。

圖6 優(yōu)化磁體組合之后的電推進(jìn)器磁場(chǎng)位形分布 Fig.6 The magnetic field distribution for the optimized assembly of permanent magnets

基于仿真得到的三維磁場(chǎng)位形，進(jìn)行了幾十萬個(gè)離子加速的測(cè)試模擬，給定離子初始位置和速度分布，求解離子運(yùn)動(dòng)控制方程，追蹤離子運(yùn)動(dòng)軌跡。離子加速的測(cè)試模擬中采用的是相對(duì)論離子運(yùn)動(dòng)方程：

式中：P是離子動(dòng)量；q是離子電荷量；r是離子在三維空間的位置坐標(biāo)；γ是相對(duì)論因子；mp是靜止離子質(zhì)量；c是光速；B是離子所在位置的磁感應(yīng)強(qiáng)度，而B隨著離子運(yùn)動(dòng)在空間維度中不斷更新，由計(jì)算網(wǎng)格點(diǎn)上的磁場(chǎng)通過雙線性插值得到。

4 模擬結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比

圖7是模擬結(jié)果與實(shí)驗(yàn)結(jié)果的對(duì)比，兩者形狀基本吻合。觀察實(shí)驗(yàn)點(diǎn)火結(jié)果發(fā)現(xiàn)，放電室右端出口附近的等離子體羽流周圍較亮，中心區(qū)域偏暗。而理論上講，如果右端出口處的磁場(chǎng)較大，離子會(huì)被磁場(chǎng)有效束縛[14]，等離子體羽流中心密度高，應(yīng)較周圍區(qū)域偏亮。因此，實(shí)驗(yàn)點(diǎn)火的結(jié)果與理論預(yù)期的結(jié)果不一致。其原因可能是放電室右端出口區(qū)域的磁場(chǎng)較弱，使得較重的離子不能被有效約束而是向周圍甩出，造成羽流中心的離子密度低，因而亮度偏暗。

圖7 測(cè)試粒子模擬結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比 Fig.7 The comparison of results between test particle simulation and discharge experiment

5 結(jié)束語

本文基于等效磁荷方法，建立了專門針對(duì)螺旋波等離子體推進(jìn)器的三維磁場(chǎng)仿真模型。該模型的計(jì)算方法簡(jiǎn)單，計(jì)算效率高，且使用靈活。在驗(yàn)證磁場(chǎng)仿真模型正確性基礎(chǔ)上，模擬出圓柱形與圓環(huán)形永磁體組合產(chǎn)生的磁場(chǎng)位形以及磁場(chǎng)強(qiáng)度分布特征。給定不同的磁鐵構(gòu)造，通過大量仿真計(jì)算，給出放電室內(nèi)的三維磁場(chǎng)位形和磁場(chǎng)強(qiáng)度分布，從中優(yōu)選出滿足螺旋波等離子體推進(jìn)器的磁鐵構(gòu)造。

三維磁場(chǎng)仿真模型和磁鐵構(gòu)造方案還有待在電推力測(cè)試和羽流特性試驗(yàn)中不斷改進(jìn)和優(yōu)化，其中，解決出口磁場(chǎng)強(qiáng)度偏弱導(dǎo)致羽流發(fā)散的問題，是提高電推進(jìn)器工作效率的關(guān)鍵。

（References）

[1] 于達(dá)仁,武志文,鄂鵬,等.穩(wěn)態(tài)等離子體發(fā)動(dòng)機(jī)磁場(chǎng)設(shè)計(jì)的發(fā)展及其展望[J].中國(guó)工程科學(xué),2005,7(4): 22-29 Yu Daren,Wu Zhiwen,E Peng,et al.Research status and development of the magnetic field design for stationary plasma thruster[J].Engineering Science,2005,7(4): 22-29

[2] 夏廣慶,王冬雪,薛偉華,等.螺旋波等離子體推進(jìn)研究進(jìn)展[J].推進(jìn)技術(shù),2011,32(6): 857-863 Xia Guangqing,Wang Dongxue,Xue Weihua,et al.Progress on the research of helicon plasma thruster[J].Journal of Propulsion Technology,2011,32(6): 857-863

[3] Chen F F,Boswell R W.Helicons: the past decade[J].IEEE Transactions on Plasma Science,1997,25(6): 1245-1257

[4] 房同珍.螺旋波激發(fā)等離子體源的原理和應(yīng)用[J].物理,1999,28(3): 162-167 Fang Tongzhen.Principle and applications of helicon wave excited plasma sources[J].Physics,1999,28(3): 162-167

[5] Chen F F.Physics of helicon discharges[J].Physics Plasmas,1996,3(5): 1783-1793

[6] 李波,王一白,張普卓,等.VASIMR 中螺旋波等離子體源設(shè)計(jì)[J].北京航空航天大學(xué)學(xué)報(bào),2012,38(6): 720-725 Li Bo,Wang Yibai,Zhang Puzhuo,et al.Design of helicon plasma source for VASIMR[J].Journal of Beijing University of Aeronautics and Astronautics,2012,38(6): 720-725

[7] 趙華,任瓊英,張紹華,等.螺旋波等離子體推力器機(jī)制的實(shí)驗(yàn)研究[C]//2013年（第九屆）中國(guó)電推力器學(xué)術(shù)研討會(huì).蘭州,2013: 260-265

[8] 任瓊英,丁永杰,彭輝,等.螺旋波電推進(jìn)器氙氣工質(zhì)電離率實(shí)驗(yàn)測(cè)量分析[C]//2013年（第九屆）中國(guó)電推力器學(xué)術(shù)研討會(huì).蘭州,2013: 282-293

[9] 劉慶海,任瓊英,王俊峰,等.螺旋波電推力器氙氣、氮?dú)庥鹆鞅容^實(shí)驗(yàn)研究[C]//2013年（第九屆）中國(guó)電推力器學(xué)術(shù)研討會(huì).蘭州,2013: 274-281

[10] 趙華,嚴(yán)世林,任瓊英,等.螺旋波電推力器氙氣工質(zhì)電離狀態(tài)轉(zhuǎn)換光譜實(shí)驗(yàn)研究[C]//2013年（第九屆）中國(guó)電推力器學(xué)術(shù)研討會(huì).蘭州,2013: 294-300

[11] Charles C,Boswell R W.Laboratory evidence of supersonic ion beam generated by a current-free “helicon” double layer[J].Physics of Plasma,2004,11(4): 1706-1714

[12] Charles C,Boswell R W.Time development of a current-free double-layer[J].Physics of Plasma,2004,11(8): 3808-3812

[13] Charles C,Boswell R W.Effect of exhaust magnetic field in a helicon double-layer thruster operating in xenon[J].IEEE Transactions on Plasma Science,2008,36(5): 2141-2146

[14] Ebersohn F H,Girimaji S S.Magnetic nozzle plasma plume: review of crucial physical phenomena,AIAA 2012-4274[R]

[15] Winglee R,Ziemba T,Giersch L.Simulation and laboratory validation of magnetic nozzle effects for the high power helicon thruster[J].Physics of Plasmas,2007,14(6): 3501

[16] 湯雙清,陳習(xí)坤,唐波.永磁體空間磁場(chǎng)的分析計(jì)算及其在永磁磁力軸承中的應(yīng)用[J].大學(xué)物理,2005,24(3): 32-36 Tang Shuangqing,Chen Xikun,Tang Bo.Analysis and calculation on spatial magnetic field emanated from a permanent magnet and its application to permanent magnetic bearing[J].College Physics,2005,24(3): 32-36

[17] 李金,鄭小林,侯文生,等.用于磁定位的永磁體磁場(chǎng)仿真與實(shí)驗(yàn)研究[J].系統(tǒng)仿真學(xué)報(bào),2009,21(18): 5919-5922 Li Jin,Zheng Xiaolin,Hou Wensheng,et al.Simulation and experimental research of magnetic field produced by permanent magnet used in magnetic localization[J].Journal of System Simulation,2009,21(18): 5919-5922

[18] Chen F F.Permanent magnet helicon source for ion propulsion[J].IEEE Transactions on Plasma Science,2008,36(5): 2095-2110