軌道平行柵推力器工作機(jī)理的數(shù)值研究

汪 忠，李艷武，車清論

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

1 計算模型

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

圖1 平行雙柵電場加速電離層等離子體圖

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

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

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

圖2 計算域及邊界條件圖

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

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

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

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

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

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

時間步長的選取必須滿足：

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

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

2 結(jié)果分析

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

表1 柵極參數(shù)

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

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

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

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

圖3 電勢分布圖

圖4 粒子軌跡圖

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

圖5 軸線處電勢隨z軸變化曲線圖

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

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

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

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

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

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

3 結(jié)論

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

[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ù)的評述[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ā)動機(jī)光學(xué)系統(tǒng)粒子模擬[J].北京航空航天大學(xué)學(xué)報，2006，32（4）：382-386.

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

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