霍爾推力器羽流對太陽翼的濺射影響規律

，，，，*

1. 上海空間推進研究所，上海 201112 2. 電子科技大學 電子學院，成都 611731 3. 上海空間發動機工程技術研究中心，上海 201112

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

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

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

1 模型

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

1)以PIC算法求解計算粒子在每個時間步長內的運動狀態(力、速度、位移等物理參數)；

2)完成各個計算粒子的物理參數對周圍節點的權重分配；

3)在完成運動狀態的更新后，以DSMC碰撞模型對粒子間的碰撞過程發生概率及碰撞類型進行判斷；

4)更新本時間步長內，部分發生碰撞的粒子的速度和加速度，以此來實現羽流中等離子體粒子的運動及碰撞過程的模擬。

在羽流-翼板相互作用過程中，以濺射模型來判斷粒子在碰撞太陽翼后的運動狀態以及計算濺射產額，以此統計羽流粒子對太陽翼表面作用的分布數據。具體模型細節將在第1.1和1.2節中介紹。

1.1 等離子體輸運過程

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

FTotal=FE+FE,AF+Fd

(1)

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

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

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

(2)

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

綜上，需要對2個參量進行節點分配計算：電荷和數密度。為此，引入單元粒子法。

1.1.1 單元粒子法(PIC)

在PIC模擬中，通常要用節點的統計參數來代表網格內的粒子狀態，因此需要以某種原則，將網格內所有粒子的關注參數分配到節點上。一般采用權重分配，常用的處理方法是：對于二維問題，采用面積權重法；對于三維問題，采用體積權重法。以二維問題為例，電量為q的離子分配給網格節點i的電量為：

(3)

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

圖1 物理參數在節點分配的方法Fig.1 Assignment strategy of the physical parameters to mesh node

可以利用節點上各項物理參數(如電荷、數密度)等求解相應的運動狀態參數(電場力、濃度梯度力)。

1.1.2 直接蒙特卡洛法(DSMC)

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

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

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

3) Xe-Xe+CEX碰撞，該類型碰撞是影響羽流發散最為重要的碰撞，自由程在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碰撞截面公式中含有經驗參數，本文根據已有推力器試驗件進行了修正，所以該公式與先前文獻不完全一致。

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

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

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

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

離子通量密度為：

Γi=npci

(4)

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

(5)

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

Γi=Γe

(6)

于是有：

(7)

因此，Xe+離子入射能量為進入鞘層時的動能加上鞘層電勢能，即

(8)

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

1.2.2 濺射產額計算

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

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

(9)

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

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

(10)

該公式中的各項參數取值依然可以參考文獻[10-11]。

2 驗證試驗

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

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

1)工況1(標準工況)，放電電壓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半徑內的10個角度位置進行電流密度的測量，具體見圖2(b)。

各算法下的計算結果與試驗結果的對比見圖3。在沒有修正離子擴散模型前，離子在軸向輸運過程中只有受到碰撞才會向兩側運動，但在羽流遠場中，濃度擴散作用與電場遷移作用相當，如果采用原有擴散模型會有很大的誤差，所以圖3的黑色曲線數據會偏離試驗結果較多。

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

圖2 驗證試驗的系統布置Fig.2 Systematic diagram of the rectification test

圖3 各工況下試驗與計算結果對比Fig.3 Comparisons of test and calculation results in each case

3 計算結果與討論

針對衛星在軌工作情況，進行在多種推力器布置下的太陽翼濺射產額分布的數值計算，具體推力器的工作參數和布置方位見圖4。

針對上述20組工況，進行霍爾推力器羽流的相關參數計算，羽流在穩態下的各物理參數分布結果見圖5，該結果作為濺射計算的輸入條件。濺射質量密度的分布計算結果見圖6。

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

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

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

續圖6Fig.6 Continued

根據圖6，可以獲得霍爾推力器羽流對太陽翼的3個主要規律：

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

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

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

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

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

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

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

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

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

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

3) 濺射作用的等效預估苛刻：玻璃磨砂作用要比濺射作用更為嚴重，因為磨砂是利用硬度較高的金剛砂等物質進行表面打磨以及用相關化學物質進行表面殘渣清理，以造成玻璃的不透明，這種處理方式明顯針對性較強，阻光機制更為有效，所以以濺射機制等效為磨砂機制是一種苛刻條件假設。

4 結束語

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

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

2)霍爾推力器羽流對太陽翼透光率的影響僅在某些惡劣區域比較嚴重(見圖7)，而其他區域對透光率的影響較低，通常可以忽略。