基于蒙特卡羅射線追蹤的航天器空間環境因素暴露通量統一計算方法

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

（1. 可靠性與環境工程技術重點實驗室； 2. 北京衛星環境工程研究所：北京 100094）

0 引言

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

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

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

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

1 空間環境因素暴露通量計算方法分析

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

1.1 原子氧通量計算方法

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

1.2 微流星體撞擊通量計算方法

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

1.3 氣動力計算方法

軌道殘余大氣分子與航天器發生動量交換產生氣動力，是低軌航天器在軌運行除地球非理想球形引力外的其他主要攝動力之一[22]。Fritsche等人建立了氣動力計算方法[23-24]。計算中，將大氣視為稀薄氣體，其分子運動符合自由分子流描述，則氣動力計算可以采用試驗粒子蒙特卡羅方法：首先根據航天器幾何構型設置計算區域（長方體），模擬試驗粒子在區域內的運動，當試驗粒子與航天器表面碰撞時，與表面按照一定的表面作用模型進行動量交換，并從表面反射；計算粒子與表面的動量交換，然后以反射點為起點，繼續追蹤粒子的運動，直到粒子逸出計算區域為止；通過模擬大量的粒子運動，可以統計得到航天器所受大氣氣動力參數。

1.4 太陽光壓計算方法

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

1.5 輻射劑量計算方法

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

1.6 放氣污染計算方法

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

2 空間環境因素暴露通量統一計算方法的提出

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

3 空間環境因素暴露通量統一計算方法流程

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

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

1）構建航天器三維模型，劃分表面網格。根據三維模型設置長方體包圍盒，作為模擬邊界。

2）設置環境因素模型參數。根據計算精度等要求設置模擬粒子數等其他全局參數。

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

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

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

6）對所有邊界進行模擬后，統計得到暴露通量，并輸出計算結果。

4 空間環境因素暴露通量統一計算方法的仿真校驗

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

4.1 原子氧暴露通量

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

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

4.2 大氣阻力

大氣阻力可表示為[35]

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

表1 大氣阻力系數計算值與理論值對比Table 1 Comparison between theoretical and calculated values of drag coefficients

4.3 太陽光壓

太陽光壓可表示為[36]

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

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

4.4 微流星體撞擊通量

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

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

4.5 等效輻射屏蔽厚度

用統一算法對文獻[27]中邊長400 mm、厚度3 mm的空心鋁質立方體的屏蔽厚度進行計算，得到其等效屏蔽厚度為3.55 mm，與文獻值一致。空心鋁質立方體等效屏蔽厚度分布如圖5所示。從圖中可以看出，正方體8個頂點處等效屏蔽厚度較大，接近5.2 mm；而中心處等效屏蔽厚度較小，為3 mm，符合幾何形態規律。

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

4.6 分子污染返回流

利用統一算法進行圓盤放氣污染返回流的計算，結果為 1.33×10-6，與文獻值[32]相比，相對誤差為0.3%，二者符合較好。

5 結論

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

利用統一算法的計算過程中，僅對表面劃分一次網格，即可用于多種空間環境因素暴露通量的計算，因此對于復雜構型航天器，可顯著節省幾何模型建模時間，提高空間環境效應分析效率，便于與其他系統級、任務級分析軟件集成。

下一步將利用真實在軌或地面試驗數據對該方法進行對比驗證分析，以檢驗其有效性。