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

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

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

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

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

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

1 仿真方法

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

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

(1)

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

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

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

2 實驗系統(tǒng)

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

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

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

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

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

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

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

圖3 計算域示意圖Fig.3 Schematic diagram of computational domain

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

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

為了出口分布的準確性，針對離子推力器在實驗中的工作情況進行了仿真，并與實驗結(jié)果進行對比。仿真結(jié)果如圖5所示，實驗數(shù)據(jù)由Faraday探針測得，實驗在真空羽流實驗設(shè)備(PES)[7，20]中開展?？梢钥闯觯抡娼Y(jié)果與實驗結(jié)果符合較好。因為近場羽流中的電流密度大部分來自于推力器直接噴出的高速離子，因此近場電流密度仿真結(jié)果與實驗結(jié)果相符能夠說明仿真所用的推力器出口參數(shù)是合理的。

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

圖5 羽流場電流密度仿真與實驗對比Fig.5 Comparison of current density in plume flow field between simulation and experiment

4 計算結(jié)果分析

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

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

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

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

圖6 軸向位置熱效應實驗和仿真對比Fig.6 Comparison of thermal effect in axial direction between simulation and experiment

圖7 距離推力器出口0.5 m和0.9 m處徑向熱流對比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)，實驗測量值全面高于平均滯止熱流仿真值和平均熱流仿真值。平均滯止熱流相對誤差絕對值在18.4%以內(nèi)，平均熱流相對誤差絕對值在26.5%以內(nèi)，誤差在可接受范圍。但同時能夠看出，在0.9 m處徑向小角度內(nèi)誤差普遍較高，說明仿真中采用的束流模型在0.5 m后下降速度較快，與實際分布相差較大。

圖8(a)給出了熱流計處于距離推力器出口0.9 m處流場中的一價Xe離子分布情況?？梢钥闯觯瑢τ谝粌rXe離子來說，熱流計的影響范圍僅在熱流計后約0.1 m的范圍內(nèi)，熱流計前方流場基本不受影響，且對流場整體影響較小。

圖8 距離推力器出口0.9 m處放置熱流計時的一價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)給出了有熱流計情況下流場中的Xe原子分布情況。可以看出，模擬熱流計對Xe原子的影響同時存在于熱流計的前方及后方部分區(qū)域。在熱流計前方，滯止于熱流計表面的Xe離子被中和成為Xe原子，使小范圍內(nèi)的粒子數(shù)密度升高，數(shù)密度級別與推力器出口相近。在熱流計后方同時出現(xiàn)數(shù)密度較低的區(qū)域，但由于Xe原子運動速度較低，且粒子數(shù)密度相比Xe離子較高，易于實現(xiàn)空間積聚，因此并未出現(xiàn)粒子數(shù)密度為0的區(qū)域。

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

圖9 距離推力器出口0.9 m處放置熱流計時的一價Xe離子、Xe原子及無熱流計時數(shù)密度對比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é) 論

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

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

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

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

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

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