，，，，*

1. 上?？臻g推進(jìn)研究所，上海 201112 2. 電子科技大學(xué) 電子學(xué)院，成都 611731 3. 上?？臻g發(fā)動機(jī)工程技術(shù)研究中心，上海 201112

霍爾推力器羽流是一種由高速離子和電子組成的等離子體射流，是形成推力、電荷中和的主要區(qū)域，然而若羽流中介入任何懸浮物體時(shí)，必然會發(fā)生相應(yīng)的質(zhì)量、動量以及能量的交換，特別地，有些交換過程具有一定的破壞性。尤其在空間應(yīng)用工況中，羽流對航天器的影響是不可忽略的。根據(jù)國外已有報(bào)道，羽流對航天器的影響主要表現(xiàn)在高速離子對各元器件表面的濺射作用[1-3]，尤其是光學(xué)敏感元器件太陽電池翼的表面，由于濺射對太陽翼表面材料二氧化硅具有一定的削蝕改性作用，可能導(dǎo)致透光率的降低，從而影響太陽翼對光能的吸收，降低衛(wèi)星系統(tǒng)的能量供給量。

關(guān)于羽流對太陽翼的濺射影響，國內(nèi)外已有一些報(bào)道：Randolph等[4]對SPT-100霍爾推力器羽流半徑為1m的位置進(jìn)行探針診斷，認(rèn)為在該位置以外的區(qū)域，羽流對太陽翼的功率影響不會超過1%，而對翼板的濺射影響可以忽略；Fife等[5]研發(fā)了一個(gè)對羽流預(yù)估、優(yōu)化設(shè)計(jì)的仿真平臺COLISEUM，通過試驗(yàn)驗(yàn)證其計(jì)算精度后，對200W的霍爾推力器羽流進(jìn)行了模擬，結(jié)果顯示濺射對厚度的削蝕速率不超過6×10-12m/s；田立成等[6]對SPT-100霍爾推力器建立了束流計(jì)算模型以及衛(wèi)星敏感區(qū)的理論模型，并利用該模型對SPT-100進(jìn)行翼板布置角在30°～50°工況的數(shù)值計(jì)算，認(rèn)為在45°布置后的濺射與沉積作用可以忽略。實(shí)際上，上述學(xué)者或是通過試驗(yàn)，或是通過仿真，均以少量工況的預(yù)估為主，給出了一些布置工況下的預(yù)估，雖然為濺射影響的分析提供了一定參考，但是上述工況覆蓋并不全面，沒有形成明顯的變化規(guī)律，缺乏通用性，這對衛(wèi)星快速調(diào)整推力器布置帶來一定難度。為此，開展霍爾推力器羽流對太陽翼多種布置工況的濺射作用規(guī)律研究，以補(bǔ)充這部分研究的空白。

由于涉及大量在軌工況的預(yù)估，本文采用數(shù)值模擬的方法。以較為成熟的單元粒子/蒙特卡洛混合算法(PIC/DSMC)[7-9]對霍爾推力器羽流進(jìn)行參數(shù)分布計(jì)算，需要說明的是，在該模型中引入離子擴(kuò)散機(jī)制來進(jìn)一步提升羽流輸運(yùn)過程的仿真精度，接著，結(jié)合鞘層電勢模型和Yamamura模型[10-12]來對太陽翼表面的濺射作用進(jìn)行模擬。

1 模型

首先針對本文所采用的主要物理模型和計(jì)算流程進(jìn)行介紹，共涉及兩個(gè)物理過程：等離子體輸運(yùn)過程和羽流-翼板相互作用過程。在等離子體輸運(yùn)過程中：

1)以PIC算法求解計(jì)算粒子在每個(gè)時(shí)間步長內(nèi)的運(yùn)動狀態(tài)(力、速度、位移等物理參數(shù))；

2)完成各個(gè)計(jì)算粒子的物理參數(shù)對周圍節(jié)點(diǎn)的權(quán)重分配；

3)在完成運(yùn)動狀態(tài)的更新后，以DSMC碰撞模型對粒子間的碰撞過程發(fā)生概率及碰撞類型進(jìn)行判斷；

4)更新本時(shí)間步長內(nèi)，部分發(fā)生碰撞的粒子的速度和加速度，以此來實(shí)現(xiàn)羽流中等離子體粒子的運(yùn)動及碰撞過程的模擬。

在羽流-翼板相互作用過程中，以濺射模型來判斷粒子在碰撞太陽翼后的運(yùn)動狀態(tài)以及計(jì)算濺射產(chǎn)額，以此統(tǒng)計(jì)羽流粒子對太陽翼表面作用的分布數(shù)據(jù)。具體模型細(xì)節(jié)將在第1.1和1.2節(jié)中介紹。

1.1 等離子體輸運(yùn)過程

等離子體羽流的輸運(yùn)過程，又稱為廣義雙極擴(kuò)散運(yùn)動，主要包括帶電粒子在電磁場中的受力運(yùn)動和粒子間的擴(kuò)散運(yùn)動。對于霍爾推力器羽流而言，可以忽略自感磁場和外加磁場，只考慮空間電荷的電場作用。每一個(gè)網(wǎng)格內(nèi)，每一個(gè)帶電粒子在運(yùn)動中所受總合外力都將遵循如下描述：

FTotal=FE+FE,AF+Fd

(1)

式中：FE為外加電場力；FE,AF為空間電荷所形成的電場力；Fd為濃度梯度力。對于空間電勢，采用求解泊松公式的方法。

除帶電粒子在電磁場中的遷移運(yùn)動外，擴(kuò)散運(yùn)動也屬于等離子體的雙極擴(kuò)散。擴(kuò)散過程是一種對濃度梯度的響應(yīng)過程，尤其在真空環(huán)境下，這種機(jī)制就必須要考慮。盡管PIC/DSMC在計(jì)算過程中，通過對粒子溫度的熱速度計(jì)算可以描述粒子，但根據(jù)已報(bào)道的計(jì)算結(jié)果來看[13]，這種描述存在一定的誤差，考慮可能是適用于中性原子的熱速度模型對于帶電的離子來說，在求解輸運(yùn)過程時(shí)會產(chǎn)生兼容性的問題(帶電粒子主要基于對力的求解，這與熱速度的耦合在算法上會有兼容性問題)，為提高計(jì)算精度，這里采取一種經(jīng)驗(yàn)性的修正。

需要說明的是，在后文的計(jì)算中，中性氣體的速度求解依然應(yīng)用熱速度模型，而正離子的擴(kuò)散運(yùn)動采用基于濃度梯度力的計(jì)算思想(菲克定律)[14]：Xe+離子在等離子體環(huán)境中的擴(kuò)散過程，微元空間某點(diǎn)A對另一點(diǎn)B的濃度擴(kuò)散力Fd為：

(2)

式中：kd為擴(kuò)散系數(shù)，根據(jù)本文計(jì)算結(jié)果，在10-36N·m4的量級與試驗(yàn)可以有較好的吻合趨勢；Ni為離子濃度(即數(shù)密度)；l為兩點(diǎn)間的距離。

綜上，需要對2個(gè)參量進(jìn)行節(jié)點(diǎn)分配計(jì)算：電荷和數(shù)密度。為此，引入單元粒子法。

1.1.1 單元粒子法(PIC)

在PIC模擬中，通常要用節(jié)點(diǎn)的統(tǒng)計(jì)參數(shù)來代表網(wǎng)格內(nèi)的粒子狀態(tài)，因此需要以某種原則，將網(wǎng)格內(nèi)所有粒子的關(guān)注參數(shù)分配到節(jié)點(diǎn)上。一般采用權(quán)重分配，常用的處理方法是：對于二維問題，采用面積權(quán)重法；對于三維問題，采用體積權(quán)重法。以二維問題為例，電量為q的離子分配給網(wǎng)格節(jié)點(diǎn)i的電量為：

(3)

式中：Ai為各個(gè)部分面積，如圖1所示。

圖1 物理參數(shù)在節(jié)點(diǎn)分配的方法Fig.1 Assignment strategy of the physical parameters to mesh node

可以利用節(jié)點(diǎn)上各項(xiàng)物理參數(shù)(如電荷、數(shù)密度)等求解相應(yīng)的運(yùn)動狀態(tài)參數(shù)(電場力、濃度梯度力)。

1.1.2 直接蒙特卡洛法(DSMC)

本文只介紹所考慮到的碰撞類型：

1) Xe-Xe彈性碰撞，是一種自由程0.005～0.1 m內(nèi)的碰撞，發(fā)生頻率較高，不可忽略；

2) Xe-Xe+彈性碰撞，該碰撞的自由程在0.001 m左右，不可忽略；

3) Xe-Xe+CEX碰撞，該類型碰撞是影響羽流發(fā)散最為重要的碰撞，自由程在0.005 m左右，是羽流中最為重要碰撞類型。

3種類型碰撞的碰撞截面公式見表1。

表1 3種碰撞類型的碰撞截面

Table 1 Collision section equations of three types

碰撞類型碰撞截面Xe-Xe+彈性碰撞σelasticXe,Xe+ =6.42×10-16crXe-Xe+CEX碰撞σcexXe,Xe+ =(-23.10lncr+138.5)×0.8423×10-20Xe-Xe彈性碰撞σelasticXe,Xe =2.12×10-18c0.24r

需要說明的是，CEX碰撞截面公式中含有經(jīng)驗(yàn)參數(shù)，本文根據(jù)已有推力器試驗(yàn)件進(jìn)行了修正，所以該公式與先前文獻(xiàn)不完全一致。

1.2 粒子與壁面的相互作用過程

等離子體羽流與太陽翼的相互作用包括電荷與濺射兩方面，電荷方面主要是指翼板壁面所形成的鞘層電勢；而濺射方面主要是指高能離子對太陽翼壁面所產(chǎn)生的濺射產(chǎn)額。前者影響入射離子的能量，后者影響透光率。

1.2.1 懸浮電極的鞘層電勢模型

太陽翼壁面相對于霍爾推力器屬于懸浮電極，因此，該鞘層模型屬于懸浮電極下的電子鞘層模型，根據(jù)Langmuir-Child鞘層理論[15]，鞘層電勢的計(jì)算過程如下。這里存在一個(gè)假設(shè)，電極所在位置的等離子體處于均勻分布假設(shè)：np=ni=ne(電子數(shù)密度接近離子數(shù)密度)。

離子通量密度為：

Γi=npci

(4)

式中：ci為離子入射速度；np為等離子體數(shù)密度。而電子通量密度為：

(5)

對于懸浮電極來說，電子通量與離子通量相等，有：

Γi=Γe

(6)

于是有：

(7)

因此，Xe+離子入射能量為進(jìn)入鞘層時(shí)的動能加上鞘層電勢能，即

(8)

式中：vi為離子入射速度；Ei為離子撞擊壁面的動能(該物理量可轉(zhuǎn)為eV單位制)。

1.2.2 濺射產(chǎn)額計(jì)算

根據(jù)相關(guān)的濺射理論，只有當(dāng)撞擊壁面的粒子動能超過壁面材料的濺射閾值時(shí)，才會發(fā)生濺射。對于本文而言，太陽翼的壁面材料為二氧化硅，其濺射閾值與硅相近，取91eV。因此，對于入射動能超過91eV的Xe+離子：

1) 若入射角度為垂直于壁面時(shí)，濺射產(chǎn)額Y(0)為:

(9)

式中：P為與入射粒子種類相關(guān)的因子；sn(ε)為核阻止截面；se(ε)為電子阻止截面；ε為入射粒子當(dāng)前的簡并能量(無量綱)；Us為壁面材料的升華能；Eth為二氧化硅的濺射閾值,取91eV。該公式各項(xiàng)參數(shù)的具體取值可以參考文獻(xiàn)[10]和文獻(xiàn)[11]。

2)對于大多數(shù)的Xe+離子，入射角度都不是0°，但此時(shí)的濺射產(chǎn)額與Y(0)有關(guān)，那么對于入射角度為θ的Xe+離子而言，濺射產(chǎn)額Y(θ)為：

(10)

該公式中的各項(xiàng)參數(shù)取值依然可以參考文獻(xiàn)[10-11]。

2 驗(yàn)證試驗(yàn)

本文所使用的計(jì)算模型分為兩部分：1)等離子體流場輸運(yùn)模型(包含離子擴(kuò)散機(jī)制)；2)等離子體對太陽翼的濺射模型。其中，模型1是加入了新物理機(jī)制，而模型2已在以往的研究中得到驗(yàn)證。因此，本文僅針對模型1進(jìn)行試驗(yàn)驗(yàn)證以及相關(guān)參數(shù)修正。

在真空艙內(nèi)開展羽流診斷試驗(yàn)，如圖2所示。采用1.35 kW霍爾推力器作為試驗(yàn)和仿真中產(chǎn)生羽流的推力器，為對模型中經(jīng)驗(yàn)參數(shù)kd的修正具有一定覆蓋性，試驗(yàn)共選取了3組有代表性的工況：

1)工況1(標(biāo)準(zhǔn)工況)，放電電壓300 V，放電電流4.5 A，總氣體流率5.488 mg/s；

2)工況2(大推力工況)，放電電壓280 V，放電電流4.82 A，總氣體流率6.174 mg/s；

3)工況3(高比沖工況)，放電電壓380 V，放電電流3.56 A，總氣體流率3.724 mg/s。

以法拉第探針對距離推力器1 m半徑內(nèi)的10個(gè)角度位置進(jìn)行電流密度的測量，具體見圖2(b)。

各算法下的計(jì)算結(jié)果與試驗(yàn)結(jié)果的對比見圖3。在沒有修正離子擴(kuò)散模型前，離子在軸向輸運(yùn)過程中只有受到碰撞才會向兩側(cè)運(yùn)動，但在羽流遠(yuǎn)場中，濃度擴(kuò)散作用與電場遷移作用相當(dāng)，如果采用原有擴(kuò)散模型會有很大的誤差，所以圖3的黑色曲線數(shù)據(jù)會偏離試驗(yàn)結(jié)果較多。

在引入離子擴(kuò)散系數(shù)kd后，根據(jù)工況不同，該系數(shù)取值又會出現(xiàn)不同。這里給出了3種工況下最為吻合的3個(gè)kd取值，可以看出，從大推力到高比沖2種極限工況范圍內(nèi)，kd的浮動范圍在1.24～1.29×10-36N·m4。由此，在下文的數(shù)值計(jì)算中，取kd=1.26×10-36N·m4。需要說明的是，由于菲克定律是基于流體模型，而粒子模型采用這種思想會有一定誤差，所以擴(kuò)散系數(shù)kd是一種經(jīng)驗(yàn)參數(shù)，并不具備對所有工況的通用性。本文認(rèn)為在10-36N·m4量級時(shí)會有較好的吻合結(jié)果，但其他學(xué)者在使用該模型前，依然要根據(jù)實(shí)際推力器及環(huán)境條件來修正該參數(shù)。在該修正下，3種工況所對應(yīng)的平均誤差可以在8.7%左右，可滿足下文數(shù)值計(jì)算需要。

圖2 驗(yàn)證試驗(yàn)的系統(tǒng)布置Fig.2 Systematic diagram of the rectification test

圖3 各工況下試驗(yàn)與計(jì)算結(jié)果對比Fig.3 Comparisons of test and calculation results in each case

3 計(jì)算結(jié)果與討論

針對衛(wèi)星在軌工作情況，進(jìn)行在多種推力器布置下的太陽翼濺射產(chǎn)額分布的數(shù)值計(jì)算，具體推力器的工作參數(shù)和布置方位見圖4。

針對上述20組工況，進(jìn)行霍爾推力器羽流的相關(guān)參數(shù)計(jì)算，羽流在穩(wěn)態(tài)下的各物理參數(shù)分布結(jié)果見圖5，該結(jié)果作為濺射計(jì)算的輸入條件。濺射質(zhì)量密度的分布計(jì)算結(jié)果見圖6。

圖4 推力器各布置工況和計(jì)算輸入Fig.4 Layout cases and input conditions of the thruster

圖5 羽流各物理參數(shù)分布云圖Fig.5 Distribution contours of plume physical parameters

圖6 各工況下的羽流對太陽翼的濺射產(chǎn)額密度分布(累計(jì)5 000 h)Fig.6 Sputtering yield density of the plume onto solar plane in each case (sum up to 5 000 h)

續(xù)圖6Fig.6 Continued

根據(jù)圖6，可以獲得霍爾推力器羽流對太陽翼的3個(gè)主要規(guī)律：

1) 隨著推力器與翼板距離的增加，濺射在整體分布和總量兩方面都有降低的趨勢，根據(jù)總濺射量進(jìn)行曲線耦合，服從二次函數(shù)規(guī)律，見圖7(a)；

2) 隨著推力器方位角的增加，濺射在整體分布和總量兩方面都有降低的趨勢，根據(jù)總濺射量進(jìn)行曲線耦合，服從指數(shù)函數(shù)規(guī)律，見圖7(b)；

3) 太陽翼受到各工況下的濺射影響都具有一個(gè)濺射集中區(qū)域，并且這個(gè)集中區(qū)域會隨著方位角和距離的增大而向外移動。

上述3種規(guī)律的產(chǎn)生主要與翼板壁面處的離子能量、數(shù)密度以及入射角度有關(guān)，隨著距離和方位角的增加，由離子能量和數(shù)密度所帶來的影響程度均有下降趨勢，該機(jī)制主要導(dǎo)致規(guī)律1)和規(guī)律2)的產(chǎn)生。入射角度的影響機(jī)制比較復(fù)雜，在Xe+對SiO2的濺射產(chǎn)額隨角度的變化曲線中，峰值點(diǎn)出現(xiàn)在65°～80°之間，角度太小或太大都會降低濺射量，于是該機(jī)制會導(dǎo)致距離、方位角發(fā)生變化時(shí)，濺射峰值區(qū)域發(fā)生漂移，主要對規(guī)律3)有較大影響。另外，由規(guī)律1)和規(guī)律2)揭示出：推力器方位角對濺射的影響程度要高于推力器與翼板的距離。

在獲得羽流對翼板濺射分布規(guī)律的基礎(chǔ)上，將進(jìn)一步討論該影響規(guī)律對透光率的影響情況。首先，將羽流對二氧化硅的濺射作用假設(shè)為一種玻璃磨砂處理。一般地，當(dāng)玻璃的磨砂厚度在1～10 μm時(shí)，每提高磨砂深度1 μm，透光率降低約6%。假設(shè)沒有經(jīng)過濺射的玻璃透光率為92%，將玻璃表面的濺射質(zhì)量分布推導(dǎo)出翼板入紙面方向微小矩形區(qū)域的濺射深度平均值，以濺射深度表示磨砂深度。通過對濺射深度的計(jì)算，可以計(jì)算所有工況下太陽翼玻璃表面的透光率變化的分布，見圖8。圖中20種工況的布置情況與圖4一致，以太陽翼表面某位置的濺射深度計(jì)算出該位置的透光率，并且以白、黃、棕、黑4種顏色的線來表示透光率降低到初始值的百分比。

圖7 濺射總量與距離、方位角的耦合曲線Fig.7 Fitting curves of total sputtering yield vs. installation angle and distance between thruster and panel

根據(jù)圖8，對各工況下翼板透光率影響惡劣的區(qū)域只占太陽翼的一部分，那么對于整個(gè)太陽翼來說，平均透光率實(shí)際上并不會降低太多。本文對各個(gè)工況下的翼板平均透光率進(jìn)行計(jì)算后，為了將影響程度高于5%的工況和其它工況區(qū)分開來，人為地在圖8中給出一條參考曲線(紅色虛線)。該曲線的意義為：只要太陽翼不穿過該曲線(方位角不小于0°)，則可認(rèn)為平均透光率下降不會超過5%。需要說明的是，該曲線是一條由苛刻的仿真條件給出的上限預(yù)估邊界，是一條保守邊界，即可能存在穿過該曲線的翼板，其透光率不會下降5%，但是，不穿過該曲線的翼板透光率一定不會降低5%。接著，給出仿真所涉及到的所有苛刻條件設(shè)定：

圖8 太陽翼各處透光率分布(累計(jì)5 000 h)Fig.8 Distribution of the light transmittance on the solar panel (sum up to 5 000 h)

1) 翼板壁面等離子體參數(shù)預(yù)估苛刻：由驗(yàn)證試驗(yàn)與計(jì)算對比可見(圖3)，只有在80°～90°時(shí)，試驗(yàn)值才高于計(jì)算值，其它0°～70°時(shí)均是計(jì)算值高于試驗(yàn)值，而太陽翼位于推力器0.2 m以外、0°角以上，其工況大部分位于探針?biāo)鶞y的0°～70°的區(qū)域，所計(jì)算得到的翼板壁面等離子體數(shù)量和能量都會高于實(shí)際情況。

2) 鞘層無碰撞假設(shè)預(yù)估苛刻：鞘層模型采用的是離子無碰撞假設(shè)，意味著只要能夠進(jìn)入鞘層，離子均會以無能量損失的狀態(tài)撞擊翼板壁面。

3) 濺射作用的等效預(yù)估苛刻：玻璃磨砂作用要比濺射作用更為嚴(yán)重，因?yàn)槟ド笆抢糜捕容^高的金剛砂等物質(zhì)進(jìn)行表面打磨以及用相關(guān)化學(xué)物質(zhì)進(jìn)行表面殘?jiān)謇?，以造成玻璃的不透明，這種處理方式明顯針對性較強(qiáng)，阻光機(jī)制更為有效，所以以濺射機(jī)制等效為磨砂機(jī)制是一種苛刻條件假設(shè)。

4 結(jié)束語

本文采用考慮離子擴(kuò)散作用的PIC/DSMC模型以及Yamamura濺射模型對霍爾推力器羽流在軌對太陽翼表面的濺射影響規(guī)律進(jìn)行了數(shù)值分析，獲得以下主要結(jié)論：

1)霍爾推力器與太陽翼的距離、方位角的增大對羽流產(chǎn)生的濺射有降低作用，分別服從二次函數(shù)衰減和指數(shù)衰減趨勢，其中方位角對濺射的衰減作用更加顯著；

2)霍爾推力器羽流對太陽翼透光率的影響僅在某些惡劣區(qū)域比較嚴(yán)重(見圖7)，而其他區(qū)域?qū)ν腹饴实挠绊戄^低，通?？梢院雎浴?/p>