基于蒙特卡羅射線追蹤的航天器空間環(huán)境因素暴露通量統(tǒng)一計(jì)算方法

焦子龍，胡松林，姜利祥，黃建國(guó)，孫繼鵬，朱云飛

（1. 可靠性與環(huán)境工程技術(shù)重點(diǎn)實(shí)驗(yàn)室； 2. 北京衛(wèi)星環(huán)境工程研究所：北京 100094）

0 引言

航天器在軌運(yùn)行過(guò)程中，受到太陽(yáng)輻射、軌道殘余大氣、帶電粒子輻射、微流星體及空間碎片等空間環(huán)境因素的影響，導(dǎo)致分系統(tǒng)或航天器性能下降，甚至任務(wù)無(wú)法完成[1]。雖然人們對(duì)空間環(huán)境因素及其效應(yīng)的認(rèn)識(shí)有了長(zhǎng)足進(jìn)步，但隨著航天器復(fù)雜度不斷提高，新技術(shù)應(yīng)用比例不斷增長(zhǎng)，受空間環(huán)境影響而發(fā)生的故障仍不斷涌現(xiàn)。對(duì)國(guó)外航天器故障進(jìn)行統(tǒng)計(jì)分析發(fā)現(xiàn)，主要分系統(tǒng)故障原因都包含空間環(huán)境因素的影響[2-3]。同樣，對(duì)國(guó)內(nèi)遙感衛(wèi)星和CAST968平臺(tái)衛(wèi)星的故障統(tǒng)計(jì)發(fā)現(xiàn)，歸因于空間環(huán)境效應(yīng)的故障分別占其故障總數(shù)的38%[4]和46%[5]。因此，開展空間環(huán)境效應(yīng)仿真分析，了解故障機(jī)理，是極為必要的。并且，空間環(huán)境效應(yīng)仿真可減少物理試驗(yàn)成本，提高航天器研制效率。

美國(guó)、歐空局等開展了大規(guī)模的空間環(huán)境研究項(xiàng)目，并將研究成果及數(shù)據(jù)集成為軟件工具系統(tǒng)，包括SPENVIS[6]、The Environment WorkBench（EWB）[7]、Space Radiation[8]以及 Geospace Environment Data Analysis System（GEDAS）[9]等。航天器空間環(huán)境效應(yīng)分析必須結(jié)合空間環(huán)境模式和航天器自身參數(shù)，因此各分析系統(tǒng)均集成有航天器軌道計(jì)算、環(huán)境計(jì)算、效應(yīng)計(jì)算等模塊，用戶只需要輸入相關(guān)參數(shù)，即可計(jì)算評(píng)估空間環(huán)境因素對(duì)航天器造成的危害程度。

國(guó)內(nèi)也研制出多種單一空間環(huán)境因素效應(yīng)分析工具或方法，如原子氧通量分析[10-11]、污染分析[12-13]、微流星體及空間碎片分析[14-16]等。這些軟件均與國(guó)外同類軟件進(jìn)行了對(duì)比驗(yàn)證，結(jié)果符合較好。但國(guó)內(nèi)目前尚無(wú)與國(guó)外類似的集成分析工具。

空間環(huán)境效應(yīng)分析的首要條件是獲取空間環(huán)境因素暴露通量。本文對(duì)基于蒙特卡羅方法的原子氧、微流星體、太陽(yáng)光壓、氣動(dòng)力、輻射劑量、材料放氣污染等空間環(huán)境因素暴露通量的計(jì)算方法進(jìn)行分析總結(jié)，提出一種基于蒙特卡羅射線追蹤的空間環(huán)境因素暴露通量統(tǒng)一計(jì)算方法；介紹該算法的計(jì)算流程，并針對(duì)原子氧通量、大氣阻力系數(shù)、太陽(yáng)光壓、等效輻射屏蔽厚度等進(jìn)行計(jì)算，與理論值或文獻(xiàn)值進(jìn)行對(duì)比驗(yàn)證。

1 空間環(huán)境因素暴露通量計(jì)算方法分析

航天器空間環(huán)境效應(yīng)包括原子氧剝蝕、大氣阻力、總劑量效應(yīng)、單粒子效應(yīng)、顆粒高速撞擊、光壓攝動(dòng)、表面充電等。為對(duì)這些效應(yīng)進(jìn)行精確分析，首先需要獲得航天器表面遭受的空間環(huán)境因素暴露通量。

1.1 原子氧通量計(jì)算方法

波音公司建立了基于射線追蹤的航天器原子氧通量計(jì)算方法及軟件[17-20]。計(jì)算中，考慮了航天器表面幾何形狀及相對(duì)位置造成的遮擋影響，以及原子氧與表面相互作用發(fā)生的鏡面反射、漫反射、表面復(fù)合、剝蝕等行為的影響。

1.2 微流星體撞擊通量計(jì)算方法

ESA建立了微流星體撞擊通量的計(jì)算方法和軟件ESABASE/Debris[21]。計(jì)算中，將微流星體環(huán)境模型根據(jù)微流星體速度、尺寸等進(jìn)行離散化，同時(shí)將航天器表面離散化，對(duì)每個(gè)表面單元，采用射線追蹤方法計(jì)算其遭受的一定速率、一定尺寸的微流星撞擊通量，然后累加求和得到航天器整體的受撞擊通量。

1.3 氣動(dòng)力計(jì)算方法

軌道殘余大氣分子與航天器發(fā)生動(dòng)量交換產(chǎn)生氣動(dòng)力，是低軌航天器在軌運(yùn)行除地球非理想球形引力外的其他主要攝動(dòng)力之一[22]。Fritsche等人建立了氣動(dòng)力計(jì)算方法[23-24]。計(jì)算中，將大氣視為稀薄氣體，其分子運(yùn)動(dòng)符合自由分子流描述，則氣動(dòng)力計(jì)算可以采用試驗(yàn)粒子蒙特卡羅方法：首先根據(jù)航天器幾何構(gòu)型設(shè)置計(jì)算區(qū)域（長(zhǎng)方體），模擬試驗(yàn)粒子在區(qū)域內(nèi)的運(yùn)動(dòng)，當(dāng)試驗(yàn)粒子與航天器表面碰撞時(shí)，與表面按照一定的表面作用模型進(jìn)行動(dòng)量交換，并從表面反射；計(jì)算粒子與表面的動(dòng)量交換，然后以反射點(diǎn)為起點(diǎn)，繼續(xù)追蹤粒子的運(yùn)動(dòng)，直到粒子逸出計(jì)算區(qū)域?yàn)橹梗煌ㄟ^(guò)模擬大量的粒子運(yùn)動(dòng)，可以統(tǒng)計(jì)得到航天器所受大氣氣動(dòng)力參數(shù)。

1.4 太陽(yáng)光壓計(jì)算方法

太陽(yáng)光壓計(jì)算可采用與氣動(dòng)力計(jì)算類似的方法[23-24]：首先設(shè)置計(jì)算區(qū)域（長(zhǎng)方體），根據(jù)太陽(yáng)入射方向得到計(jì)算區(qū)域各面的太陽(yáng)輻射通量；然后從計(jì)算區(qū)域各面隨機(jī)選取位置作為入射光束的起點(diǎn)，光束方向?yàn)樘?yáng)入射方向，光束代表的太陽(yáng)輻射通量由蒙特卡羅抽樣總數(shù)計(jì)算，計(jì)算光束與衛(wèi)星表面部件的交點(diǎn)，并判斷距離最近的交點(diǎn)，其所在面元即為實(shí)際受太陽(yáng)光照射的面元；計(jì)算在該面元上產(chǎn)生的太陽(yáng)光壓，如果該面元的鏡面反射率大于0，則按照上述流程繼續(xù)追蹤反射光線的運(yùn)動(dòng)，并計(jì)算其產(chǎn)生的太陽(yáng)光壓；循環(huán)上述過(guò)程，直到反射光線強(qiáng)度小于規(guī)定的閾值時(shí)停止。

1.5 輻射劑量計(jì)算方法

國(guó)內(nèi)外開發(fā)了多種三維輻射劑量評(píng)估方法和軟件[25-29]。由于航天器結(jié)構(gòu)為薄壁結(jié)構(gòu)，吸收劑量計(jì)算中采用了近似直線傳輸原理[30]，即入射高能粒子在材料中沿直線運(yùn)動(dòng)，因此屏蔽厚度采用沿入射方向直線上的結(jié)構(gòu)厚度，即若高能粒子入射方向與結(jié)構(gòu)表面法線間的夾角為θ，則屏蔽厚度Leあ=L/cosθ，其中L為結(jié)構(gòu)厚度?；诖嗽?，吸收劑量的計(jì)算過(guò)程為：在星內(nèi)選定輻射劑量分析點(diǎn)，從該點(diǎn)向星外全向空間均勻分布引出若干條射線；然后計(jì)算出每條射線穿透的衛(wèi)星中各儀器設(shè)備、結(jié)構(gòu)件等單元的累積屏蔽厚度；再結(jié)合劑量-深度關(guān)系曲線，計(jì)算得到每條射線所代表的小區(qū)域中經(jīng)屏蔽后的空間輻射劑量，從而完成整星輻射劑量三維分析計(jì)算。

1.6 放氣污染計(jì)算方法

在軌航天器用有機(jī)材料造成的表面污染量及性能退化程度計(jì)算方法為[31-32]：認(rèn)為放氣分子發(fā)射及分子從表面反射均符合余弦規(guī)律，放氣點(diǎn)位置在放氣表面元均勻分布，利用蒙特卡羅方法抽樣放氣分子運(yùn)動(dòng)的起始位置和方向；然后抽樣分子運(yùn)動(dòng)的自由程，判斷該自由程代表的分子運(yùn)動(dòng)路程是否與物體表面相交，沒有經(jīng)歷碰撞的情形可以認(rèn)為自由程為無(wú)窮大；從大量的分子運(yùn)動(dòng)軌跡中抽樣計(jì)算得到污染源i和目標(biāo)表面j之間的質(zhì)量傳遞系數(shù)，將其代入矩陣方程即可求解得到表面j的入射流，則該表面的沉積速率為，其中為該表面污染分子的沉積系數(shù)。

2 空間環(huán)境因素暴露通量統(tǒng)一計(jì)算方法的提出

空間環(huán)境效應(yīng)分析是衛(wèi)星設(shè)計(jì)研制中的一項(xiàng)重要工作，呈現(xiàn)出不同空間環(huán)境因素效應(yīng)協(xié)同分析以及與衛(wèi)星任務(wù)級(jí)、系統(tǒng)級(jí)分析集成等發(fā)展趨勢(shì)。為此，本文提出構(gòu)建統(tǒng)一的基于蒙特卡羅射線追蹤的計(jì)算方法，將原子氧、微流星體、帶電粒子、放氣分子等視作模擬粒子，用射線代表粒子的運(yùn)動(dòng)軌跡，追蹤大量粒子與航天器表面的相互作用，從而得到原子氧通量、微流星體撞擊通量、氣動(dòng)力、太陽(yáng)光壓、輻射劑量、放氣污染量等。相比于單空間環(huán)境因素效應(yīng)分析工具，這種計(jì)算方法物理圖像清晰，無(wú)須航天器幾何模型在不同分析工具之間的轉(zhuǎn)換，可顯著提高分析評(píng)估的效率，有助于未來(lái)實(shí)現(xiàn)任務(wù)級(jí)、系統(tǒng)級(jí)集成分析。

3 空間環(huán)境因素暴露通量統(tǒng)一計(jì)算方法流程

本文提出的基于蒙特卡羅射線追蹤的統(tǒng)一計(jì)算方法流程如圖1所示。主要步驟如下：

圖1 航天器空間環(huán)境因素暴露通量計(jì)算流程Fig.1 Flowchart for exposure flux computation of space environmental factors

1）構(gòu)建航天器三維模型，劃分表面網(wǎng)格。根據(jù)三維模型設(shè)置長(zhǎng)方體包圍盒，作為模擬邊界。

2）設(shè)置環(huán)境因素模型參數(shù)。根據(jù)計(jì)算精度等要求設(shè)置模擬粒子數(shù)等其他全局參數(shù)。

3）對(duì)于原子氧、大氣阻力、太陽(yáng)光壓暴露通量，基于空間環(huán)境模型計(jì)算邊界處通量。例如，對(duì)于原子氧，采用NRLMSISE-00模型[33]得到某軌道位置的大氣原子氧數(shù)密度nO、溫度T等參數(shù)值，則模擬邊界處單位面積上的原子氧通量為

4）從模擬邊界處生成代表粒子的射線。對(duì)于原子氧、大氣阻力、太陽(yáng)光壓計(jì)算，模擬邊界為人為設(shè)置的包圍盒表面；對(duì)于輻射劑量、微流星體、材料放氣計(jì)算，模擬邊界為航天器表面單元。射線起點(diǎn)位置P0在邊界上均勻隨機(jī)分布，射線方向D按照環(huán)境因素特性隨機(jī)抽樣。射線可表示為

5）追蹤射線，判斷其軌跡與表面是否相交。如果相交，根據(jù)其與航天器表面的相互作用模型計(jì)算暴露通量；然后根據(jù)表面特性重新設(shè)置射線的方向，并以交點(diǎn)為新的起點(diǎn)；繼續(xù)追蹤，直到射線逃逸出限定計(jì)算區(qū)域，或者滿足一定的閾值條件后拋棄該軌跡。若不相交，則進(jìn)行下一個(gè)粒子的模擬。

6）對(duì)所有邊界進(jìn)行模擬后，統(tǒng)計(jì)得到暴露通量，并輸出計(jì)算結(jié)果。

4 空間環(huán)境因素暴露通量統(tǒng)一計(jì)算方法的仿真校驗(yàn)

空間環(huán)境因素暴露通量統(tǒng)一算法的校驗(yàn)可通過(guò)將統(tǒng)一算法計(jì)算結(jié)果與理論值（或文獻(xiàn)值）進(jìn)行對(duì)比完成。復(fù)雜幾何構(gòu)型航天器難以獲得理論值，且已有文獻(xiàn)中航天器構(gòu)型復(fù)雜、難以精確復(fù)現(xiàn)。為此，本文采用典型簡(jiǎn)單幾何模型作為模擬對(duì)象，將統(tǒng)一算法的計(jì)算結(jié)果與理論值進(jìn)行比對(duì)：以航天器幾何模型劃分的表面網(wǎng)格作為基本模擬單元，考慮幾何模型的代表性，選取文獻(xiàn)[32]中的圓盤算例作為模擬對(duì)象，圓盤直徑0.564 m，劃分三角形網(wǎng)格；統(tǒng)一算法的模擬誤差與模擬粒子數(shù)目N的平方根成反比，由于圓盤為簡(jiǎn)單幾何模型，計(jì)算量較小，所以選取模擬粒子數(shù)目為100萬(wàn)個(gè)。

4.1 原子氧暴露通量

計(jì)算圓盤不同攻角下的原子氧暴露通量。理論分析值取自文獻(xiàn)[17]，軌道原子氧數(shù)密度8.2×1011/m3，圓盤速度8000 m/s。粒子撞擊速度分布采用舍選抽樣法[34]，撞擊位置在包圍盒表面上均勻隨機(jī)分布。原子氧暴露通量的計(jì)算值和理論值對(duì)比如圖2所示。由圖可見：攻角為0°～80°時(shí)，計(jì)算值與理論值間的相對(duì)誤差小于0.1%；攻角為90°時(shí)，相對(duì)誤差為1%，計(jì)算值和理論值符合很好。但誤差隨攻角的變化規(guī)律有待進(jìn)一步研究。

圖2 圓盤算例原子氧暴露通量理論值與計(jì)算值比較Fig.2 Comparison of theoretical and calculated AO fluxes for circular disk sample

4.2 大氣阻力

大氣阻力可表示為[35]

式中：CD為大氣阻力系數(shù)；A為迎風(fēng)面面積；ρ為大氣密度；vr為航天器相對(duì)于大氣的運(yùn)動(dòng)速度；ev為航天器相對(duì)大氣運(yùn)動(dòng)方向單位矢量。航天器軌道、姿態(tài)一定時(shí)，其大氣阻力取值主要與大氣阻力系數(shù)CD有關(guān)，而影響CD的主要因素是氣體分子與航天器表面相互作用模型。本文分別采用鏡面反射模型和漫反射模型計(jì)算3.1節(jié)中圓盤算例的大氣阻力系數(shù)。理論分析值取自文獻(xiàn)[35]，圓盤表面溫度為300 K。大氣阻力系數(shù)的計(jì)算值與理論值對(duì)比如表1所示，2種反射模型下計(jì)算值與理論值間的相對(duì)誤差均僅有0.1%，二者符合較好。

表1 大氣阻力系數(shù)計(jì)算值與理論值對(duì)比Table 1 Comparison between theoretical and calculated values of drag coefficients

4.3 太陽(yáng)光壓

太陽(yáng)光壓可表示為[36]

式中：c為真空中的光速；Φ為單位時(shí)間單位面積接受的太陽(yáng)輻照度，其在地球表面的平均值約為1367 W/m2，隨衛(wèi)星與太陽(yáng)間距離不同而變化，變化幅度為每年±3.3%。不同太陽(yáng)光入射角下的太陽(yáng)光壓計(jì)算值與理論值對(duì)比如圖3所示。從圖中可以看出，計(jì)算值與理論值間的最大誤差約為0.4%，二者符合較好。

圖3 不同太陽(yáng)光入射角下的歸一化太陽(yáng)光壓計(jì)算值與理論值對(duì)比Fig.3 Comparison between theoretical and calculated value of normalized solar pressure for different incidence angles

4.4 微流星體撞擊通量

微流星體撞擊航天器表面的撞擊通量可以表示為

式中：vi為微流星體撞擊速度，是微流星體相對(duì)于航天器的速度，vi=vm-vs，vm為微流星體速度矢量，vs為航天器速度矢量；n(vm)為微流星體速度分布。文獻(xiàn)[21]給出了當(dāng)所有微流星體具有相同速度，即n(vm)=1時(shí)，航天器速度矢量與表面法線呈不同角度時(shí)，其微流星體歸一化撞擊通量的變化。圖4給出了統(tǒng)一算法計(jì)算結(jié)果與文獻(xiàn)值的對(duì)比。由圖中可以看出，計(jì)算值與文獻(xiàn)值間的誤差最大約0.4%，二者符合較好。

4.5 等效輻射屏蔽厚度

用統(tǒng)一算法對(duì)文獻(xiàn)[27]中邊長(zhǎng)400 mm、厚度3 mm的空心鋁質(zhì)立方體的屏蔽厚度進(jìn)行計(jì)算，得到其等效屏蔽厚度為3.55 mm，與文獻(xiàn)值一致?？招匿X質(zhì)立方體等效屏蔽厚度分布如圖5所示。從圖中可以看出，正方體8個(gè)頂點(diǎn)處等效屏蔽厚度較大，接近5.2 mm；而中心處等效屏蔽厚度較小，為3 mm，符合幾何形態(tài)規(guī)律。

圖5 空心鋁質(zhì)立方體等效屏蔽厚度分布Fig.5 Distribution of equivalent shielding thickness of hallow aluminum cube

4.6 分子污染返回流

利用統(tǒng)一算法進(jìn)行圓盤放氣污染返回流的計(jì)算，結(jié)果為 1.33×10-6，與文獻(xiàn)值[32]相比，相對(duì)誤差為0.3%，二者符合較好。

5 結(jié)論

本文提出了一種基于蒙特卡羅方法的空間環(huán)境因素暴露通量統(tǒng)一計(jì)算方法：將空間環(huán)境因素等價(jià)為粒子束或能量束集合，用射線代表粒子束或能量束，跟蹤射線在空間中的運(yùn)動(dòng)以及與航天器表面的相互作用，得到空間環(huán)境因素暴露通量。針對(duì)簡(jiǎn)單幾何構(gòu)型，通過(guò)典型算例（圓盤）計(jì)算了原子氧暴露通量、大氣阻力系數(shù)、太陽(yáng)光壓、微流星體撞擊通量、等效輻射屏蔽厚度、分子污染返回流等，與理論值進(jìn)行對(duì)比，均符合較好，證明本文提出的統(tǒng)一算法能夠正確模擬相關(guān)物理過(guò)程，可用于近地軌道、深空等各種航天器的空間環(huán)境因素暴露通量計(jì)算。

利用統(tǒng)一算法的計(jì)算過(guò)程中，僅對(duì)表面劃分一次網(wǎng)格，即可用于多種空間環(huán)境因素暴露通量的計(jì)算，因此對(duì)于復(fù)雜構(gòu)型航天器，可顯著節(jié)省幾何模型建模時(shí)間，提高空間環(huán)境效應(yīng)分析效率，便于與其他系統(tǒng)級(jí)、任務(wù)級(jí)分析軟件集成。

下一步將利用真實(shí)在軌或地面試驗(yàn)數(shù)據(jù)對(duì)該方法進(jìn)行對(duì)比驗(yàn)證分析，以檢驗(yàn)其有效性。