汪 忠，李艷武，車清論

（蘭州空間技術(shù)物理研究所 真空技術(shù)與物理重點(diǎn)實(shí)驗(yàn)室，蘭州 730000）

0 引言

電推進(jìn)技術(shù)是目前空間推進(jìn)的先進(jìn)技術(shù)，已成為航天器性能先進(jìn)性的標(biāo)志[1-2]。電推進(jìn)的原理為空間推進(jìn)技術(shù)的發(fā)展提供了進(jìn)一步的空間。對(duì)近地軌道有一種可能性，就是利用地球高層大氣（電離層中F層）中自然存在的稀薄等離子體作為空間發(fā)動(dòng)機(jī)的推進(jìn)劑。太陽(yáng)輻射的能量致使高層大氣部分電離，形成稀薄的等離子體。原理上如果航天器上的適當(dāng)裝置可以直接加速這些軌道環(huán)境中的等離子體，則航天器將受到一個(gè)反向推力。有關(guān)這個(gè)設(shè)想的基本方案是在航天器上添加一個(gè)加有電壓的平行雙柵結(jié)構(gòu)，當(dāng)航天器在電離層飛行時(shí)，電離層中的等離子體以6 000～8 000 m/s的相對(duì)速度進(jìn)入航天器上的平行柵，其中正離子被雙柵之間的電場(chǎng)加速，電子則被柵間電場(chǎng)排斥，由于電子質(zhì)量遠(yuǎn)小于正離子，所以產(chǎn)生的阻力可以忽略。2006年美國(guó)Northrop Grumman公司提出了利用平行柵加速軌道環(huán)境離子的電推進(jìn)技術(shù)方案[3]（The Ambient At?mosphere Ion Thruster，AAIT）。此后，美國(guó)喬治亞州理工學(xué)院研究小組報(bào)告了研究結(jié)果，建立了利用軌道環(huán)境離子的平行柵推力器的數(shù)學(xué)模型，并計(jì)算了推力、阻力等[4-5]。平行柵電場(chǎng)加速軌道離子能否產(chǎn)生推力是研究平行柵推力器要解決的一個(gè)核心問(wèn)題，用二維PIC方法對(duì)軌道離子被雙柵引出的過(guò)程進(jìn)行數(shù)值模擬[6]。

1 計(jì)算模型

加有電壓的平行雙柵的結(jié)構(gòu)類似于離子推力器的屏柵和加速柵，雙柵上打有規(guī)則分布的同軸小孔，因柵極上的每個(gè)孔都具有軸對(duì)稱性，采用二維軸對(duì)稱模型[7-8]對(duì)單孔進(jìn)行模擬，如圖1所示。

圖1 平行雙柵電場(chǎng)加速電離層等離子體圖

計(jì)算域r、z分別為徑向和軸向位置，R0、Z0分別為計(jì)算域徑向和軸向長(zhǎng)度，rf為前柵孔半徑，rb為后柵孔半徑，Φ為電勢(shì)，下邊界為柵孔軸線，坐標(biāo)原點(diǎn)選在計(jì)算域左下角，計(jì)算區(qū)域的大小為7.5 mm×1.1 mm，劃分為300×44的正交等距網(wǎng)格，網(wǎng)格大小為0.025 mm，滿足△z≤λD，△r≤λD的要求，其中λD為等離子體的德拜長(zhǎng)度。計(jì)算域的邊界條件分為電場(chǎng)邊界條件和粒子運(yùn)動(dòng)邊界條件兩個(gè)部分，如圖2所示。計(jì)算域左右邊界的電場(chǎng)邊界條件取為Φ=0，前柵電勢(shì)取為0，后柵電勢(shì)取為-1 000 V，上下邊界的電場(chǎng)邊界條件取為Neumann型邊界條件。粒子運(yùn)動(dòng)邊界條件分為兩種。上下邊界為反射邊界，左右邊界和柵極為吸收邊界。模擬過(guò)程中，如果粒子超出左右邊界或者撞到柵極，則將粒子刪除，如果達(dá)到上下邊界，則將粒子鏡面反射回計(jì)算域，即粒子z向速度不變，r向速度相反。

雙柵產(chǎn)生的靜電勢(shì)和運(yùn)動(dòng)等離子體產(chǎn)生的自洽電勢(shì)可通過(guò)求解Possion方程來(lái)得到：

式中：e為電子電量；ε0為真空介電常數(shù)；ni為離子數(shù)密度；ne為電子數(shù)密度。ni通過(guò)以網(wǎng)格點(diǎn)為中心的一個(gè)網(wǎng)格大小的方形區(qū)域的離子個(gè)數(shù)除以網(wǎng)格面積求得。將電子視為熱平衡流體，ne通過(guò)Boltzmann方程求得。采用有限差分法求解出電勢(shì)后，可求出電場(chǎng)強(qiáng)度：

圖2 計(jì)算域及邊界條件圖

PIC方法的基本思想為：通過(guò)求解電場(chǎng)Poisson方程得到電勢(shì)分布，用插值法求得粒子所在位置處的電場(chǎng)力，然后加速粒子，再用新的粒子位置求解新的電場(chǎng)分布，上述過(guò)程循環(huán)進(jìn)行，直至收斂。

粒子以初始位置和初始速度從計(jì)算域左邊界進(jìn)入，經(jīng)電場(chǎng)作用后再?gòu)倪吔珉x開(kāi)。在粒子初始化中，須確定粒子的初始位置與初始速度，粒子初始位置都位于計(jì)算域左邊界發(fā)射面上，其在發(fā)射面上位置隨機(jī)給定，并保證粒子在整個(gè)圓形發(fā)射面上均勻分布，粒子初始發(fā)射速度z向分量由航天器在軌速度和Maxwell速度確定：

發(fā)射速度徑向分量滿足Maxwell速度分布：

式中：vM為Maxwell速度；θ是在0到π間均勻分布的隨機(jī)數(shù)。

粒子加速方法采用蛙跳格式，粒子運(yùn)動(dòng)差分方程為：

式中：m為粒子質(zhì)量；v為粒子速度；F為電場(chǎng)力；q為電量；x為粒子位置，電場(chǎng)強(qiáng)度由與粒子鄰近的4個(gè)網(wǎng)格點(diǎn)上的電場(chǎng)強(qiáng)度插值得到。

時(shí)間步長(zhǎng)的選取必須滿足：

式中：ωp為等離子體頻率；選取tp=1.0×10-9s。

為簡(jiǎn)化問(wèn)題，在粒子模擬過(guò)程中，只考慮軌道等離子體環(huán)境中的一價(jià)單原子氧離子、電子，不考慮中性原子。等離子體中存在多種粒子之間的碰撞，如帶電粒子之間的庫(kù)侖碰撞，離子與中性粒子間的電荷交換碰撞等，但在低軌道環(huán)境中，粒子密度較低，分子自由程較大，因此，不考慮粒子間的碰撞。

2 結(jié)果分析

模擬程序需要輸入的參數(shù)包括：柵極的幾何參數(shù)、柵極電壓、等離子體參數(shù)。所取參數(shù)如表1、2所列，其中等離子體參數(shù)以300 km高度軌道的電離層等離子體為例[9]。

表1 柵極參數(shù)

表2 等離子體參數(shù)

在300 km高度，主要的離子成分為單原子氧離子，在模擬過(guò)程中，只考慮O+而忽略其他離子成分（如He+）。平行柵上游離子密度niu是正常電離層300 km高度的數(shù)值，下游離子密度根據(jù)上游粒子密度計(jì)算得到，具體如式（7）：

式中：vu為上游離子速度，即離子與航天器的相對(duì)速度為7 700 m/s；vd為下游離子速度；ΔV為雙柵電勢(shì)差；m為氧離子質(zhì)量。

從圖3、4中可以看出，柵極間的電勢(shì)線幾乎垂直于z軸，因此在z軸方向上形成極強(qiáng)的電場(chǎng)將進(jìn)入的離子加速引出，而當(dāng)離子被引出后柵，電勢(shì)則開(kāi)始上升，離子受到與z軸相反方向的電場(chǎng)開(kāi)始減速，這是因?yàn)榧佑胸?fù)高壓的后柵在等離子體中將形成較大尺寸的鞘層（幾分米至幾米），鞘層處的電位與空間電位相當(dāng)，可認(rèn)為是0電位，后柵和鞘層之間形成與z相反方向的電場(chǎng)。值得說(shuō)明的是在模擬中，計(jì)算域z軸方向的尺度應(yīng)與鞘層尺度相當(dāng)，但這將耗費(fèi)相當(dāng)長(zhǎng)的機(jī)時(shí)。因此，模擬過(guò)程未將計(jì)算域z軸取到鞘層尺度，但z軸方向尺度已經(jīng)取得較大，可認(rèn)為右邊界處即為鞘層，這樣的做法對(duì)模擬結(jié)果的影響可以忽略。

圖3 電勢(shì)分布圖

圖4 粒子軌跡圖

如圖5所示，左邊界電位為0，在平行柵中，電位下降△U，到平行柵下游，電位開(kāi)始上升，到右邊界，電位恢復(fù)到0，離子在這一過(guò)程中電場(chǎng)勢(shì)能沒(méi)有改變，因此并沒(méi)有從平行柵電場(chǎng)中獲得能量，即離子在這一過(guò)程中沒(méi)有被加速。

圖5 軸線處電勢(shì)隨z軸變化曲線圖

King等[5]提出在平行柵前增加一個(gè)離子勢(shì)能提升區(qū)，用來(lái)提升離子勢(shì)能，這樣，當(dāng)離子到達(dá)平行柵下游的鞘層處，多余的勢(shì)能轉(zhuǎn)化為動(dòng)能，即離子最終被加速噴出，航天器獲得反向推力。

圖6和圖7是在離子勢(shì)能提升區(qū)將電勢(shì)提升到50 V后的模擬結(jié)果。比較圖4和圖7可以看出，加了離子勢(shì)能提升區(qū)后，離子在加速柵中的聚焦變的明顯。區(qū)的辦法，平行柵則可以加速軌道中的離子以獲得反向推力。

圖6 加離子勢(shì)能提升區(qū)的電勢(shì)分布圖

圖7 加離子勢(shì)能提升區(qū)的粒子軌跡圖

如圖8所示，由于離子勢(shì)能提升，到下游鞘層處，電勢(shì)恢復(fù)到0時(shí)離子還剩余部分動(dòng)能，即離子被加速噴出，航天器將獲得反向推力。

圖8 加離子勢(shì)能提升區(qū)后軸線處電勢(shì)隨z軸變化曲線圖

3 結(jié)論

模擬了不加離子勢(shì)能提升區(qū)和加上離子勢(shì)能提升區(qū)兩種情況下離子被平行柵電場(chǎng)引出的情況。模擬結(jié)果表明：加電壓的平行柵無(wú)法直接加速軌道離子，利用在平行柵前增加一個(gè)離子勢(shì)能提升

[1]張?zhí)炱?，張雪?空間電推進(jìn)技術(shù)及應(yīng)用新進(jìn)展[J].真空與低溫，2013，19（4）：187-193.

[2]張?zhí)炱剑荜怀?，孫小菁，等.小衛(wèi)星領(lǐng)域應(yīng)用電推進(jìn)技術(shù)的評(píng)述[J].真空與低溫，2014，20（4）：187-192.

[3]Dressler G A.Spacecraft propulsive device using ambient up?per atmospheric constituents for reaction mass[C]//42nd Joint Propulsion Conference&Exhibit，2006.

[4]King S T，Walker M L R，Chianese S G.Ambient Atmosphere Ion Thruster（AAIT）Proof-of-Concept Modeling[C]//47thJoint Propulsion Conference&Exhibit，2011.

[5]King S T，Walker M L R，Silvio G.Chianese.Atmospheric Elec?tric Propulsion Mission Performance Tool[J].Journal of Space?craft and Rockets，2014，51（3）：931-937.

[6]邵福球.等離子體離子模擬[M].北京：科學(xué)出版社，2002：9-11.

[7]劉暢，湯海濱，顧佐，等.基于PIC方法的離子發(fā)動(dòng)機(jī)光學(xué)系統(tǒng)粒子模擬[J].北京航空航天大學(xué)學(xué)報(bào)，2006，32（4）：382-386.

[8]劉洋，李娟，楚豫川，等.考慮電荷交換的柵極區(qū)離子流數(shù)值模擬[J].真空與低溫，2011，17（4）：198-208.

[9]莊洪春.宇航空間環(huán)境手冊(cè)[M].北京：中國(guó)科學(xué)技術(shù)出版社，2000：116-117.