張建華， 李晶華， 尤鳳儀， 鄭鴻儒

(北京航空航天大學(xué) 宇航學(xué)院, 北京 100083)

電推力器以其比沖高、壽命長(zhǎng)和可控性強(qiáng)等優(yōu)點(diǎn)在深空探測(cè)中發(fā)揮著重要作用，已經(jīng)被應(yīng)用于姿態(tài)控制、軌道轉(zhuǎn)移、深空探測(cè)等空間任務(wù)中[1-2]。電推進(jìn)技術(shù)的在軌應(yīng)用為航天器空間環(huán)境研究帶來(lái)了許多無(wú)法忽視的新難題，電推進(jìn)羽流與航天器的相互作用是其中重要的內(nèi)容之一。

離子推力器羽流主要由離子推力器高速噴出的高能工質(zhì)離子與中和器發(fā)射的電子構(gòu)成的束流等離子體及與電荷交換產(chǎn)生的交換電荷等離子體兩部分組成。羽流會(huì)撞擊航天器表面，對(duì)航天器產(chǎn)生力效應(yīng)、熱效應(yīng)和污染效應(yīng)等，影響航天器部件的正常工作和壽命，嚴(yán)重時(shí)會(huì)導(dǎo)致航天任務(wù)失敗[3]。其中，熱效應(yīng)主要是由于離子推力器羽流會(huì)直接或者返流與航天器表面發(fā)生碰撞，這種碰撞會(huì)引起航天器敏感材料熱變形，而熱流的累積也會(huì)造成航天器重要元器件的損壞，進(jìn)而導(dǎo)致航天器重要部件功能的喪失。因此，開(kāi)展離子推力器的羽流熱效應(yīng)研究十分必要。

目前，國(guó)內(nèi)外電推進(jìn)的診斷測(cè)量主要利用靜電探針、發(fā)射探針、阻滯勢(shì)分析儀(RPA)、法拉第探針等接觸式診斷設(shè)備直接放置在等離子羽流中，對(duì)電子密度、空間電位、離子能量分布、等離子體電勢(shì)、電流密度等加以測(cè)量，然而關(guān)于等離子體的熱效應(yīng)研究正式發(fā)表的文獻(xiàn)很少。2007年，上海交通大學(xué)周志雄等[4]對(duì)真空條件下氣體與固體表面的相互作用過(guò)程的熱適應(yīng)系數(shù)進(jìn)行了數(shù)值模擬仿真研究。2008年，美國(guó)奧斯汀大學(xué)Berisford等[5]以螺旋波等離子體為研究對(duì)象，對(duì)氬放電過(guò)程進(jìn)行了熱效應(yīng)診斷實(shí)驗(yàn)研究，用熱電偶和熱輻射探針?lè)謩e安裝在氣流的下游和真空艙尾部測(cè)量熱流值。

對(duì)于推力器熱效應(yīng)的仿真研究主要集中于化學(xué)推力器。2014年，王黎珍等[6]對(duì)化學(xué)推力器真空羽流熱效應(yīng)的計(jì)算模型進(jìn)行了修正和誤差分析。目前，尚未見(jiàn)針對(duì)電推進(jìn)羽流熱效應(yīng)的仿真研究，而仿真研究對(duì)于電推力器的實(shí)際在軌應(yīng)用具有重要意義。如果分析過(guò)于保守，則導(dǎo)致熱防護(hù)過(guò)度；估計(jì)不足，則會(huì)使熱量積聚，可能導(dǎo)致設(shè)備失效。在電推進(jìn)羽流熱效應(yīng)分析中，如何準(zhǔn)確獲得羽流流場(chǎng)及高能粒子與表面作用的精確模型是2個(gè)主要難點(diǎn)。本文針對(duì)LIPS-200型離子推力器地面真空艙內(nèi)羽流熱效應(yīng)實(shí)驗(yàn)中的部分工況，采用三維粒子網(wǎng)格-直接模擬蒙特卡羅(PIC-DSMC)方法進(jìn)行數(shù)值模擬研究，并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比分析。

1 仿真方法

本文涉及到的仿真算例均采用EX-PWS(Extension of Plume Work Station)[7]軟件實(shí)現(xiàn)。EX-PWS是羽流工作站(PWS)[8]軟件的電推進(jìn)擴(kuò)展程序包，使用PIC方法[9]處理電場(chǎng)解算、粒子運(yùn)動(dòng)等過(guò)程，使用DSMC方法[10]處理粒子間碰撞。EX-PWS軟件使用三維非結(jié)構(gòu)網(wǎng)格適應(yīng)復(fù)雜空間環(huán)境。此外，為提高計(jì)算效率，軟件使用了MPI并行計(jì)算方法。

為了減少計(jì)算量，EX-PWS將等離子體中的電子視為流體，將離子及推力器出口噴出的Xe原子視為計(jì)算粒子。電子運(yùn)動(dòng)采用準(zhǔn)中性假設(shè)[11]和Boltzmann關(guān)系式描述。電子溫度使用等溫模型，忽略磁場(chǎng)影響。Boltzmann關(guān)系式[12]可表述為

(1)

式中：ne為電子密度；nref為電勢(shì)φ為零處的參考電子密度，這里取1×1016m-3；Te為電子溫度，這里取為常數(shù)3.5 eV；k為Boltzmann常數(shù)。

粒子間碰撞包括彈性碰撞及電荷交換碰撞，其是影響流場(chǎng)的一個(gè)主要因素。在軟件中，粒子間碰撞使用DSMC方法進(jìn)行處理。忽略Xe原子、Xe離子與碳原子之間的彈性碰撞和電荷交換碰撞。各種粒子間的碰撞截面選取與文獻(xiàn)[13]相同。

本文在模擬環(huán)境壓強(qiáng)時(shí)采用虛擬粒子的方法，將背景氣體視為虛擬粒子，僅參與碰撞計(jì)算。當(dāng)計(jì)算粒子運(yùn)動(dòng)到某一網(wǎng)格中將要與背景Xe原子發(fā)生碰撞時(shí)，按照設(shè)置的環(huán)境壓強(qiáng)給出數(shù)密度，按照環(huán)境溫度根據(jù)Maxwellian 分布給出熱運(yùn)動(dòng)速度，碰撞后背景粒子狀態(tài)不變，計(jì)算粒子運(yùn)動(dòng)速度及方向，按照碰撞模型發(fā)生變化。計(jì)算粒子與背景粒子間的碰撞類(lèi)型包括彈性碰撞及電荷交換碰撞。

2 實(shí)驗(yàn)系統(tǒng)

實(shí)驗(yàn)系統(tǒng)由真空羽流實(shí)驗(yàn)設(shè)備、離子推力器系統(tǒng)、熱流傳感器系統(tǒng)和Faraday探針系統(tǒng)組成，本文主要介紹測(cè)量熱流的傳感器系統(tǒng)。熱流傳感器系統(tǒng)主要由總熱流計(jì)、輻射熱流計(jì)、信號(hào)轉(zhuǎn)換裝置和電流信號(hào)采集系統(tǒng)組成。Schmidt-Boelter[14-16]熱流計(jì)是基于空間溫度梯度原理，利用熱電堆測(cè)量溫度差的方法，進(jìn)而獲得熱流的大小。實(shí)驗(yàn)測(cè)量系統(tǒng)的示意圖如圖 1所示，本次實(shí)驗(yàn)過(guò)程中使用的熱流傳感器實(shí)物圖如圖 2所示。

本文進(jìn)行離子推力器羽流熱效應(yīng)診斷采用的熱流計(jì)為Schmidt-Boelter熱流計(jì)。熱流計(jì)的量程為10 kW/m2，精度滿量程±3%，因此熱流計(jì)的診斷誤差在±0.3 kW/m2左右。同時(shí)需要注意的是，熱流計(jì)最高容許傳感器溫度為204℃，在具體的實(shí)驗(yàn)過(guò)程中采用兩路熱電偶對(duì)熱流傳感器的溫度進(jìn)行監(jiān)控，并采用降溫處理使熱流傳感器的溫度保持在合理的溫度范圍內(nèi)。

圖1 離子推力器羽流熱效應(yīng)測(cè)量系統(tǒng)示意圖Fig.1 Schematic diagram of ion thruster plume thermal effect measurement system

圖2 兩種熱流傳感器實(shí)物圖Fig.2 Photo of two kinds of heat flow sensor

3 網(wǎng)格及邊界條件

仿真計(jì)算區(qū)域設(shè)定為包括推力器前方1.5 m和推力器后方0.5 m、直徑2.0 m的臥式圓柱體，網(wǎng)格數(shù)為200萬(wàn)左右，沿x方向劃分為24個(gè)并行分區(qū)。如圖3所示，計(jì)算域內(nèi)部有一個(gè)小圓柱體用來(lái)模擬推力器，模擬熱流計(jì)為直徑30 mm、長(zhǎng)40 mm的小圓柱體，模擬熱流計(jì)與推力器間的相對(duì)位置如圖4所示。

在處理粒子與邊界面之間的相互作用時(shí)，采用適應(yīng)系數(shù)衡量來(lái)流粒子速度對(duì)壁面溫度的適應(yīng)程度。在粒子與表面發(fā)生碰撞時(shí)，程序會(huì)生成一個(gè)隨機(jī)數(shù)，其大于適應(yīng)系數(shù)，則粒子采用鏡面反射，否則，采用漫反射。文獻(xiàn)[17]中提出當(dāng)離子能量高于400 eV時(shí)，適應(yīng)系數(shù)將接近于1.0，為保守起見(jiàn)，本文中熱流計(jì)表面適應(yīng)系數(shù)取為0.9，即認(rèn)為90%的粒子進(jìn)行漫反射，10%的粒子進(jìn)行鏡面反射。計(jì)算域邊界設(shè)置為自由邊界，當(dāng)粒子運(yùn)動(dòng)到邊界后即從流場(chǎng)中移除該粒子。

圖3 計(jì)算域示意圖Fig.3 Schematic diagram of computational domain

圖4 推力器與模擬熱流計(jì)相對(duì)位置示意圖Fig.4 Schematic diagram of relative position of thruster and simulated heat flow meter

實(shí)驗(yàn)中采用的推力器為蘭州空間技術(shù)物理研究所研制的1.3 kW的LIPS-200型離子推力器[18]。推力器工質(zhì)為氙氣，工作時(shí)的流量為1.37 mg/s,柵極電壓為1 000 V。仿真中設(shè)置推力器電離率為90%，羽流發(fā)散半角為16°，二價(jià)離子占離子總數(shù)的10%。推力器仿真參數(shù)設(shè)置如表 1所示。推力器出口射出的粒子包括Xe原子、一價(jià)離子和二價(jià)離子3種，出口流量分布設(shè)置為高斯分布[19]。計(jì)算中忽略了中和器帶來(lái)的影響。

為了出口分布的準(zhǔn)確性，針對(duì)離子推力器在實(shí)驗(yàn)中的工作情況進(jìn)行了仿真，并與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。仿真結(jié)果如圖5所示，實(shí)驗(yàn)數(shù)據(jù)由Faraday探針測(cè)得，實(shí)驗(yàn)在真空羽流實(shí)驗(yàn)設(shè)備(PES)[7，20]中開(kāi)展?？梢钥闯?，仿真結(jié)果與實(shí)驗(yàn)結(jié)果符合較好。因?yàn)榻鼒?chǎng)羽流中的電流密度大部分來(lái)自于推力器直接噴出的高速離子，因此近場(chǎng)電流密度仿真結(jié)果與實(shí)驗(yàn)結(jié)果相符能夠說(shuō)明仿真所用的推力器出口參數(shù)是合理的。

表1 LIPS-200型離子推力器仿真基本參數(shù)

圖5 羽流場(chǎng)電流密度仿真與實(shí)驗(yàn)對(duì)比Fig.5 Comparison of current density in plume flow field between simulation and experiment

4 計(jì)算結(jié)果分析

本節(jié)針對(duì)LIPS-200型離子推力器，開(kāi)展電推進(jìn)羽流熱效應(yīng)仿真分析，并與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。本文對(duì)羽流熱效應(yīng)的計(jì)算綜合考慮了熱流計(jì)置于流場(chǎng)中對(duì)羽流的影響，尤其是高能粒子撞擊熱流計(jì)表面后被中和并損失能量成為慢速中性原子后，仍然參與Xe原子流場(chǎng)的計(jì)算，使得仿真條件與實(shí)驗(yàn)工況盡量吻合。計(jì)算中，時(shí)間步長(zhǎng)設(shè)置為1×10-7s，粒子權(quán)重設(shè)置為5×108，環(huán)境及各壁面溫度為300 K，環(huán)境壓強(qiáng)為2 mPa， 并認(rèn)為推力器及熱流計(jì)接地，所有表面電勢(shì)設(shè)置為0 V。

圖6給出了軸向位置上(距離推力器出口0.5～1.0 m)實(shí)驗(yàn)和仿真的對(duì)比。其中，平均滯止熱流仿真值曲線顯示粒子與熱流計(jì)表面進(jìn)行了完全的熱交換(熱適應(yīng)系數(shù)為1.0)；平均熱流仿真值曲線即帶有熱適應(yīng)系數(shù)的仿真熱流值(本文中的適應(yīng)系數(shù)取值為0.9)。需要注意的是，仿真值中的平均指的是在30 mm的模擬熱流計(jì)表面進(jìn)行能量統(tǒng)計(jì)后在表面上進(jìn)行平均，避免在單個(gè)網(wǎng)格中由于粒子數(shù)較少造成熱流值畸高。

可以看出，在軸線上仿真值與實(shí)驗(yàn)測(cè)量值均沿著與推力器的距離升高而降低，實(shí)驗(yàn)測(cè)量值與平均滯止熱流仿真值相近，最大誤差出現(xiàn)在0.9 m處，與平均滯止熱流仿真值的誤差值為17.0%，與平均熱流仿真值的誤差為25.3%，在等離子體測(cè)量領(lǐng)域認(rèn)為是可以接受的。但同時(shí)可以看出，本文算例中采用的適應(yīng)系數(shù)偏低，當(dāng)粒子能量全部傳遞(適應(yīng)系數(shù)為1.0)，即粒子滯止時(shí)熱流仿真值與實(shí)驗(yàn)測(cè)量值符合較好。仿真與實(shí)驗(yàn)誤差產(chǎn)生的可能原因包括束流模型偏差、適應(yīng)系數(shù)選取及實(shí)驗(yàn)誤差等方面。

圖7(a)給出了距離推力器出口0.5 m處徑向熱流隨角度的分布情況。可以看出，在角度為0°～15°的羽流束流區(qū)域內(nèi)，實(shí)驗(yàn)測(cè)量值基本小于平均滯止熱流仿真值，而高于平均熱流仿真值。平均滯止熱流相對(duì)誤差絕對(duì)值在11.1%以內(nèi)，平均熱流相對(duì)誤差絕對(duì)值在21.1%以內(nèi)，誤差較小。平均滯止熱流相對(duì)誤差較大的位置在15°左右，此處為束流半角，即束流區(qū)和非束流區(qū)的交界處，此處的等離子體數(shù)密度較為稀薄，仿真和實(shí)驗(yàn)的難度都很大，因此相對(duì)誤差較高。但同時(shí)也能看出，無(wú)論是平均滯止熱流仿真值還是平均熱流仿真值，其沿徑向下降的速度與實(shí)驗(yàn)測(cè)量值相比較為平緩，說(shuō)明在采用Maxwell模型處理粒子與表面能量交換時(shí)使用的適應(yīng)系數(shù)不應(yīng)是固定值，而是應(yīng)該隨著角度進(jìn)行變化。從粒子撞擊的角度來(lái)說(shuō)，正入射時(shí)粒子垂直撞擊熱流計(jì)表面，絕大部分能量轉(zhuǎn)化為熱能，而斜入射粒子存在一定的出射速度，從而降低了能量交換。因此，當(dāng)角度升高時(shí)，適應(yīng)系數(shù)應(yīng)隨之降低。

圖6 軸向位置熱效應(yīng)實(shí)驗(yàn)和仿真對(duì)比Fig.6 Comparison of thermal effect in axial direction between simulation and experiment

圖7 距離推力器出口0.5 m和0.9 m處徑向熱流對(duì)比Fig.7 Comparison of radial heat flow at 0.5 m and 0.9 m from thruster exit

圖7(b)給出了距離推力器出口0.9 m處徑向熱流隨角度的分布。可以看出，在羽流束流區(qū)域內(nèi)，實(shí)驗(yàn)測(cè)量值全面高于平均滯止熱流仿真值和平均熱流仿真值。平均滯止熱流相對(duì)誤差絕對(duì)值在18.4%以內(nèi)，平均熱流相對(duì)誤差絕對(duì)值在26.5%以內(nèi)，誤差在可接受范圍。但同時(shí)能夠看出，在0.9 m處徑向小角度內(nèi)誤差普遍較高，說(shuō)明仿真中采用的束流模型在0.5 m后下降速度較快，與實(shí)際分布相差較大。

圖8(a)給出了熱流計(jì)處于距離推力器出口0.9 m處流場(chǎng)中的一價(jià)Xe離子分布情況?？梢钥闯?，對(duì)于一價(jià)Xe離子來(lái)說(shuō)，熱流計(jì)的影響范圍僅在熱流計(jì)后約0.1 m的范圍內(nèi)，熱流計(jì)前方流場(chǎng)基本不受影響，且對(duì)流場(chǎng)整體影響較小。

圖8 距離推力器出口0.9 m處放置熱流計(jì)時(shí)的一價(jià)Xe離子和Xe原子分布Fig.8 Number density distribution of Xe atom and monovalent xenon ion when heat flow meter is located at 0.9 m away from thruster exit

圖 8(b)給出了有熱流計(jì)情況下流場(chǎng)中的Xe原子分布情況?？梢钥闯?，模擬熱流計(jì)對(duì)Xe原子的影響同時(shí)存在于熱流計(jì)的前方及后方部分區(qū)域。在熱流計(jì)前方，滯止于熱流計(jì)表面的Xe離子被中和成為Xe原子，使小范圍內(nèi)的粒子數(shù)密度升高，數(shù)密度級(jí)別與推力器出口相近。在熱流計(jì)后方同時(shí)出現(xiàn)數(shù)密度較低的區(qū)域，但由于Xe原子運(yùn)動(dòng)速度較低，且粒子數(shù)密度相比Xe離子較高，易于實(shí)現(xiàn)空間積聚，因此并未出現(xiàn)粒子數(shù)密度為0的區(qū)域。

圖9給出了距離推力器出口0.9 m處放置熱流計(jì)后軸線上的一價(jià)Xe離子數(shù)密度、Xe原子數(shù)密度分布與無(wú)熱流計(jì)時(shí)數(shù)密度的對(duì)比?？梢钥闯觯瑢?duì)于一價(jià)Xe離子來(lái)說(shuō)，影響較大的區(qū)域在熱流計(jì)測(cè)量平面后的0.1 m范圍內(nèi)，在熱流計(jì)后方距離較近區(qū)域接近為0，然后逐漸恢復(fù)至正常量級(jí)，但在計(jì)算域內(nèi)仍低于無(wú)熱流計(jì)的情況。這是由于一價(jià)Xe離子的速度較高，約為3.8×104m/s，在無(wú)干擾的情況下基本沿直線運(yùn)動(dòng)，很難受到影響，當(dāng)高速運(yùn)動(dòng)的Xe離子撞擊熱流計(jì)表面滯止后，被中和成為Xe原子，同時(shí)能量傳遞給熱流計(jì)，因此造成了熱流計(jì)后方小范圍內(nèi)粒子數(shù)為0。而對(duì)于Xe原子，由于一部分Xe離子被中和后返回Xe原子流場(chǎng)中，造成了熱流計(jì)前方0.1 m范圍內(nèi)Xe原子數(shù)密度升高約2個(gè)量級(jí)，但由于在稀薄氣體中粒子間的碰撞概率較低，小范圍內(nèi)Xe原子數(shù)密度的升高對(duì)流場(chǎng)影響不大。

圖9 距離推力器出口0.9 m處放置熱流計(jì)時(shí)的一價(jià)Xe離子、Xe原子及無(wú)熱流計(jì)時(shí)數(shù)密度對(duì)比Fig.9 Comparison of xenon atom and monovalent xenon ion number density distribution on axis with or without heat folw meter located at 0.9 m away from thruster exit

5 結(jié) 論

本文針對(duì)LIPS-200型離子推力器羽流熱效應(yīng)進(jìn)行了三維數(shù)值模擬，采用準(zhǔn)中性假設(shè)和Boltzmann關(guān)系式描述電子分布，采用Maxwell模型對(duì)粒子與表面能量交換過(guò)程，得到如下結(jié)論：

1) 采用混合PIC-DSMC方法可以對(duì)電推進(jìn)羽流熱效應(yīng)進(jìn)行有效的數(shù)值模擬分析，推力器出口軸線上滯止熱流仿真與實(shí)驗(yàn)對(duì)比，誤差小于17.0%。

2) 在進(jìn)行徑向熱效應(yīng)仿真時(shí)，由于粒子具有一定的出射速度，熱適應(yīng)系數(shù)不是恒定不變的，而是應(yīng)隨著角度的升高而減小。

3) 在進(jìn)行三維粒子模擬時(shí)，模擬熱流計(jì)的介入對(duì)流場(chǎng)整體影響較小。對(duì)Xe離子的影響主要集中在熱流計(jì)后方0.1 m范圍內(nèi)，并出現(xiàn)了小范圍內(nèi)數(shù)密度為0的區(qū)域。熱流計(jì)后方Xe離子數(shù)密度低于無(wú)熱流計(jì)時(shí)的流場(chǎng)。

4) 模擬熱流計(jì)對(duì)于Xe原子的影響同時(shí)存在于熱流計(jì)前方和后方，在熱流計(jì)前方中性粒子數(shù)密度升高約2個(gè)量級(jí)，與推力器出口附近的數(shù)密度相近。熱流計(jì)后方Xe原子數(shù)密度低于無(wú)熱流計(jì)時(shí)的流場(chǎng)。

綜上所述，本文在驗(yàn)證羽流流場(chǎng)的基礎(chǔ)上，對(duì)比了推力器出口1.0 m范圍內(nèi)的羽流熱效應(yīng)實(shí)驗(yàn)測(cè)量結(jié)果，驗(yàn)證了仿真模型精確性，可為進(jìn)一步分析衛(wèi)星及敏感儀器受熱流的影響提供可信的輸入條件。