基于環(huán)境工質(zhì)的超低軌吸氣式螺旋波電推進(jìn)器仿真分析

丁 亮，彭毓川，鄭慧奇，唐振宇，任瓊英，趙 華，秦 瑋

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

0 引言

大氣阻力是制約航天器在超低軌道工作壽命的主要原因之一。而如果航天器能夠充分利用軌道環(huán)境中的殘余大氣作為工質(zhì)，經(jīng)過收集、電離、加速后為航天器提供空間推進(jìn)動(dòng)力，維持航天器的軌道速度，將大大減少航天器所需攜帶的推進(jìn)劑量，延長(zhǎng)航天器的在軌運(yùn)行壽命。Nishiyama[1]提出利用電子回旋共振（ECR）技術(shù)的吸氣式離子發(fā)動(dòng)機(jī)概念，以期維持航天器的超低軌長(zhǎng)時(shí)間運(yùn)行。Pekker[2]研究霍爾電推進(jìn)技術(shù)，以超低軌道上殘余氣體為工質(zhì)，理論上可以將航天器維持在80～90 km高度的軌道上。這2種方案都是先將軌道殘余大氣經(jīng)過壓縮收集到高壓容器中，然后向推進(jìn)器供氣。其主要難點(diǎn)一是如何高效收集軌道殘余大氣并實(shí)現(xiàn)上萬倍的密度壓縮；二是軌道殘余大氣的主要成分為氧原子，而在離子電推進(jìn)和霍爾電推進(jìn)這些有電極加速機(jī)制的推力器上，氧原子對(duì)電極的剝蝕作用將導(dǎo)致電推進(jìn)器的壽命大幅縮短，難以實(shí)現(xiàn)長(zhǎng)時(shí)間的軌道維持[3-5]。

螺旋波電推進(jìn)采用射頻功率耦合，無電極腐蝕過程，從原理上能夠避免電推進(jìn)器的上述壽命缺陷。國(guó)內(nèi)外針對(duì)螺旋波電推進(jìn)技術(shù)已開展了大量研究[6-12]：Chen最先開展了螺旋波等離子體的研究，分析了螺旋波電離、放電模式、功率吸收與模式激發(fā)過程；Charles等人將螺旋波等離子體應(yīng)用到了電推進(jìn)領(lǐng)域，并分析了推進(jìn)模式和加速機(jī)制。但目前還未見關(guān)于螺旋波等離子體在吸氣式電推進(jìn)技術(shù)中應(yīng)用的相關(guān)報(bào)道。

為了驗(yàn)證超低軌道高度環(huán)境下吸氣式螺旋波電推進(jìn)的可行性，本文進(jìn)行以原子氧為工質(zhì)，綜合考慮碰撞、激發(fā)、電離等過程，通過13.56 MHz射頻加熱，由磁噴口完成推力輸出全流程的仿真模擬；建立特定形態(tài)吸氣式螺旋波電推進(jìn)結(jié)構(gòu)模型，施加不同功率的輸入，分析功率沉積、等離子體參數(shù)分布和推力輸出。

1 模擬分析的原理

在超低軌道環(huán)境中，以180 km高度為例，氣體主要成分為原子氧（O）。收集的軌道殘余氣體通過射頻波加熱后呈等離子體態(tài)，故可采用等離子體物理模型模擬該氣體流。在等離子體模型中考慮電子（e）、氧離子（O+）、激發(fā)態(tài)氧原子（Os）以及本底氧原子（O）4種粒子的電化學(xué)反應(yīng)：

對(duì)于彈性碰撞，有

對(duì)于激發(fā)過程，有

對(duì)于退激過程，有

對(duì)于電離過程，有

同時(shí)，取吸附系數(shù)為1時(shí)壁面電化學(xué)反應(yīng)為：

在所采用的等離子體模型中，電子密度ne、平均電子能量nε的漂移?擴(kuò)散方程為：

式(7)～(9)中：μe為電子遷移率；E為等離子體電場(chǎng)強(qiáng)度；De為電子擴(kuò)散系數(shù)；Re為電子源項(xiàng)；με為電子能量遷移率；Dε為電子能量擴(kuò)散系數(shù)；Γe為電子通量；Rε為非彈性碰撞損失項(xiàng)；Qrh為波加熱項(xiàng)貢獻(xiàn)；J為等離子電流；E*為射頻波場(chǎng)。其中，De、Dε、με的計(jì)算式為

假設(shè)有M種反應(yīng)凈生成電子，則電子源項(xiàng)為

假設(shè)有P種非彈性碰撞反應(yīng)導(dǎo)致能量變化，則非彈性碰撞損失項(xiàng)為

式(10)～(12)中：Te為電子溫度；xj為碰撞目標(biāo)粒子的份額；kj為類反應(yīng)的速度系數(shù)；Nn為總的中性粒子密度；Δεj為j類型反應(yīng)轉(zhuǎn)移的能量。

對(duì)于氧離子，可通過等離子體電化學(xué)反應(yīng)計(jì)算獲得其密度分布。對(duì)于非帶電的激發(fā)態(tài)原子氧，采用重粒子的對(duì)流擴(kuò)散方程進(jìn)行模擬計(jì)算，

其中：ρ為電荷密度；wk為總粒子數(shù)的份額；u為等離子體宏觀流體速度矢量；jk為激發(fā)態(tài)原子氧擴(kuò)散通量；Rk為激發(fā)態(tài)原子氧的生成速率源項(xiàng)，由電化學(xué)反應(yīng)決定。假設(shè)有n種帶電重粒子，則等離子體電場(chǎng)的計(jì)算方程為

其中：ε0為真空介電常數(shù)；εr為相對(duì)介電常數(shù)；V為電勢(shì)；電荷密度ρ由等離子體化學(xué)反應(yīng)過程耦合計(jì)算獲得；q為單位電荷；Zk為帶電粒子電荷數(shù)；nk為帶電重粒子密度。

采用Maxwell方程組模擬加熱等離子體的電磁場(chǎng)：

式(15)～(18)中：B為磁感應(yīng)強(qiáng)度；J為等離子體電流；A為磁矢勢(shì)；D為電位移矢量；假設(shè)等離子體各向同性，電導(dǎo)率σ為

其中：me為電子靜止質(zhì)量；γe為電子?中性粒子碰撞頻率。

采用可膨脹層流模型方程組模擬原子氧流體變化：

式(20)～(22)中：U為原子氧流體的速度矢量；P為流體壓力；I為單位張量；μ為等離子體黏度系數(shù)，通過各粒子的份額綜合計(jì)算獲得；Cp為原子氧流體熱容；Q為熱源項(xiàng)，由等離子體電化學(xué)反應(yīng)方程耦合獲得；ρ是由各粒子的份額綜合計(jì)算獲得的密度值。

整體上將洛倫茲、Maxwell方程組耦合的熱流磁流體方程對(duì)裝置吸入氣體的熱場(chǎng)、流速場(chǎng)進(jìn)行模擬，并通過流體的物性與前述的等離子體模型中的粒子密度、溫度參數(shù)進(jìn)行耦合，可形成一個(gè)評(píng)估波加熱等離子體效應(yīng)對(duì)宏觀流體力學(xué)性能影響的物理模型，進(jìn)而得出裝置的推力效果。該模型通過有限元COMSOL-MULTIPHYSIC軟件平臺(tái)進(jìn)行求解。模型網(wǎng)格劃分為：等離子體壁面邊界采用四邊形網(wǎng)格，線圈區(qū)域采用四邊形網(wǎng)格，其他區(qū)域采用三角形網(wǎng)格；總網(wǎng)格數(shù) 750 000，時(shí)間步長(zhǎng) 1×10-5s。

當(dāng)達(dá)到熱平衡時(shí)，推進(jìn)器產(chǎn)生的推力Fthrust由以下方程計(jì)算：

質(zhì)量流量為

沉積功率為

式 (23)～(25)中：Vo、ρo、So分別為出口的流體速度、密度、截面積；ΔW為沉積在原子氧流體中的功率。

2 幾何模型及計(jì)算條件

吸氣式螺旋波電推進(jìn)器的幾何模型如圖1所示，其裝置尺寸見表1。

圖1 吸氣式螺旋波電推進(jìn)器仿真結(jié)構(gòu)Fig. 1 Simulated structure of the air-breathing helicon plasma propulsion

表1 吸氣式螺旋波電推進(jìn)器仿真結(jié)構(gòu)裝置尺寸Table 1 Sizes of simulated structure of the air-breathing helicon plasma propulsion

推進(jìn)器運(yùn)行軌道高度設(shè)計(jì)為180 km。在該軌道環(huán)境下的殘余氣體粒子密度為15.4×1015m-3，飛行器的飛行速度為7800 m/s，則氣體收集口的粒子通量為 1.2×1020m-2·s-1。設(shè)計(jì)裝置的入口截面積為0.020 1 m2，相應(yīng)的，單位時(shí)間吸入的粒子數(shù)為 2.4×1018s-1，約為 4.0×10-6mol/s，對(duì)應(yīng)為 178 sccm。根據(jù)粒子密度推算大氣壓力約為1.9×10-5Pa，環(huán)境溫度約為650 K。用以約束粒子的偏轉(zhuǎn)磁場(chǎng)，最強(qiáng)處設(shè)計(jì)為1500 Gs，磁場(chǎng)位型如圖2所示。軸向0點(diǎn)設(shè)定為收集口喉部位置，為在直觀反映計(jì)算結(jié)果的同時(shí)減小計(jì)算量，計(jì)算過程示以半剖面。

圖2 約束磁場(chǎng)位型Fig. 2 The constrained magnetic field

3 不同驅(qū)動(dòng)功率下的仿真計(jì)算結(jié)果

200 W的輸入驅(qū)動(dòng)功率水平下，沉積在原子氧流體中的功率為2.478 W，出口處質(zhì)量流量為7.511 9×10-6kg/s，計(jì)算得到推力值約為 6 mN。相應(yīng)的，達(dá)到穩(wěn)態(tài)時(shí)功率沉積、電子密度分布、電子溫度分布、流體速度分布如圖3所示，電子密度峰值1.24×1018m-3，電子溫度峰值3.25 eV。

800 W的輸入驅(qū)動(dòng)功率水平下，沉積在原子氧流體中的功率為6.452 7 W，出口處質(zhì)量流量為7.511 9×10-6kg/s，計(jì)算得到推力值為 9.85 mN。相應(yīng)的，達(dá)到穩(wěn)態(tài)時(shí)功率沉積、電子密度分布、電子溫度分布、流體速度分布如圖4所示，電子密度峰值3.52×1018m-3，電子溫度峰值 3.45 eV。

圖3 200 W 輸入功率仿真結(jié)果Fig. 3 Simulation result with input power of 200 W

圖4 800 W 輸入功率仿真結(jié)果Fig. 4 Simulation result with input power of 800 W

1200 W的輸入驅(qū)動(dòng)功率水平下，沉積在原子氧流體中的功率為8.285 4 W，出口處質(zhì)量流量為7.511 9×10-6kg/s，計(jì)算得到推力值為 11.16 mN。相應(yīng)的，達(dá)到穩(wěn)態(tài)時(shí)功率沉積、電子密度分布、電子溫度分布、流體速度分布如圖5所示，電子密度峰值5.15×1018m-3，電子溫度峰值 3.61 eV。

圖5 1200 W 輸入功率仿真結(jié)果Fig. 5 Simulation result with input power of 1200 W

2000 W的輸入驅(qū)動(dòng)功率水平下，沉積在原子氧流體中的功率為11.71 W，出口處質(zhì)量流量為7.511 9×10-6kg/s，計(jì)算得到推力值為 13.23 mN。相應(yīng)的，達(dá)到穩(wěn)態(tài)時(shí)功率沉積、電子密度分布、電子溫度分布、流體速度分布如圖6所示，電子密度峰值7.53×1018m-3，電子溫度峰值 3.72 eV。

流體沉積功率、推力隨驅(qū)動(dòng)功率的變化如圖7所示。

圖6 2000 W 輸入功率仿真結(jié)果Fig. 6 Simulation result with input power of 2000 W

圖7 流體沉積功率、推力隨驅(qū)動(dòng)功率的變化Fig. 7 Deposited power and thrust versus input power

4 結(jié)束語

通過以上仿真分析發(fā)現(xiàn)，隨著驅(qū)動(dòng)功率的增加，推進(jìn)器的推力和流體電子密度均顯著增加——輸入功率從200 W增加至2000 W的過程中，推力從6.00 mN增加至13.23 mN。而吸氣式螺旋波電推進(jìn)器的設(shè)計(jì)應(yīng)用軌道環(huán)境截面阻力約0.226 mN，因此該推進(jìn)器能夠滿足阻力補(bǔ)償?shù)囊蟆ｋ娮用芏确逯祻?1.24×1018m-3增加至 7.53×1018m-3，增幅較大；電子溫度峰值緩慢增加，從3.25 eV增加至3.72 eV。電磁輻照的功率沉積量較低，但沉積功率的密度較高，在 2000 W 時(shí)到達(dá)了最大值 105W/m3。以上數(shù)據(jù)表明：隨著輸入功率的增加，較多功率用于電離，而溫度升高得較少，說明能量獲得能力較低，從而導(dǎo)致推力漲幅不大，僅增加了1倍多。其中可能的原因是通道直徑較小，粒子復(fù)合率較高，輸入的功率大部分消耗在將復(fù)合后的粒子重新電離上。因此，需要考慮加大加熱段的容積來增加功率的沉積值，以增大裝置的推力。

總之，根據(jù)上述推力計(jì)算結(jié)果，在下一步工作中，將考慮加大入口截面積、增加氣體的質(zhì)量流量、增大沉積功率，繼而增大推進(jìn)器的推力。