基于射線追蹤原理的微流星體撞擊通量計算研究

焦子龍,姜利祥,李濤,孫繼鵬,黃建國,朱云飛

(1.可靠性與環(huán)境工程技術重點實驗室,北京100094；2.北京衛(wèi)星環(huán)境工程研究所,北京100094)

1 概述

微流星體一般是指直徑小于1mm或質量小于1g的流星體，主要來自于小行星和彗星。微流星體在地球附近、行星際空間、星際空間等廣泛存在。微流星體撞擊航天器可能帶來設備功能失效等多種影響。為確保航天器的微流星體撞擊防護措施有效、適當，需要首先定量評估微流星體的撞擊風險。早在20世紀60年代的阿波羅任務和禮炮號空間站任務中就開始重視并評估微流星體撞擊影響[1]。目前已有 BUMPER、ESABASE/DEBRIS、 COLLO、 BUFFER、 MDPANTO、 MODAOST、ARMOR等國內外多種空間碎片/微流星體撞擊風險評估軟件[2-4]。

本文提出了一種基于射線追蹤計算微流星體撞擊通量的方法。該方法原理簡單、易于實現(xiàn)，便于與其他空間環(huán)境因素效應分析模塊集成[5]。本文首先介紹了基于射線追蹤的微流星體撞擊通量計算方法，然后通過算例對該方法進行了對比驗證分析，最后對某空間站復雜構型進行了微流星體撞擊通量計算。

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

2.1 微流星體環(huán)境模型

本文采用的微流星體環(huán)境模型為Grün模型，它定義了距太陽一個天文單位 (1AU)處相對靜止的隨機翻轉的面積為1㎡的平板在各向同性的微流星體流中一年時間內受到的撞擊數(shù)目F(m):

由于各環(huán)境模型中使用的均是微流星體的尺寸下限，所以它對應的通量值表示大于該尺寸的所有粒子的通量總和。

微流星體相對于地球的速率分布函數(shù)為:

速率單位為km/s。

由于航天器的運動使得其不同朝向的表面微流星撞擊通量顯著不同，例如迎風朝向的表面撞擊通量較尾流方向的通量高出將近一個數(shù)量級。此外，航天器復雜的結構使得不同表面之間存在遮擋，造成微流星通量顯著不同。因此不能直接采用Grün模型計算表面的微流星通量，而應建立通用的計算方法。

2.2 射線追蹤計算方法

如圖1所示為微流星體撞擊平板時通量計算的示意圖。平板法線方向為n，運動速度為vS，微流星體速度為vm。微流星體撞擊平板的速度為vi，vi=vm-vS。

圖1 撞擊通量計算示意圖Fig.1 Schematic diagram for calculating impact flux

由圖1可知，撞擊通量可以表示為:

n(vm)為一定速度分布下微流星體速度為vm的概率。因此，有:

基于射線追蹤的微流星體撞擊通量計算流程如圖2所示。

圖2 微流星體通量計算流程圖Fig.2 Flowchart of micrometeoroid flux computation

對應圖2，計算流程如下:

(1)完成初始化工作，包括航天器表面劃分三角形網(wǎng)格，設置每個網(wǎng)格單元格需要追蹤的射線數(shù)目等參數(shù)。

(2)設置航天器軌道和姿態(tài)參數(shù)，建立計算坐標系。本文中計算坐標系如圖3所示。航天器位于+Z軸。

圖3 微流星體通量計算坐標系Fig.3 Coordinate system for calculating flux of micrometeoroids

(3)針對每個網(wǎng)格單元:

1)生成新的射線。在立體角為4π的空間內均勻隨機產(chǎn)生射線方向，其方法如下:

式中，r0和r1為 [0，1]區(qū)間的均勻分布隨機數(shù)；θ和φ分別為天頂角和方位角；x，y，z為笛卡爾坐標系下的射線方向。

2)根據(jù)微流星體環(huán)境模型隨機抽樣微流星體速率值。隨機抽樣可采用舍選抽樣法。

3)計算微流星體相對撞擊速度vi。如果vi·n＜0，則微流星體撞擊該表面單元，繼續(xù)d。否則判斷射線抽樣數(shù)目是否足夠。如果不夠，繼續(xù)執(zhí)行a。

4)抽樣射線的初始位置，其方法如下:設三角形單元三個頂點分別為A、B和C，其在系統(tǒng)坐標系中的坐標向量分別為，則射線初始位置為:

式中，r0和r1為 [0，1]區(qū)間的均勻分布隨機數(shù)。若r0+r1＞1，則r0=1-r0，r1=1-r1。

5)判斷射線是否與其他三角形單元相交。為加速計算是否存在交點，可采用八叉樹等空間劃分結構[6]。如果相交，判斷射線抽樣數(shù)目是否足夠。如果不夠，返回步驟1)。如果不相交，繼續(xù)步驟6)。

6)根據(jù)式 (3)計算通量增量ΔFi。三角形單元受到的總通量為:

式中，N為追蹤的射線數(shù)目；F0為根據(jù)式 (1)計算的微流星體通量；G為引力會聚因子。對于地球，表達式為:

式中，RE為地球半徑；h為軌道高度。

7)循環(huán)步驟1) ～6)，直到達到指定的射線數(shù)目。

(4)循環(huán)步驟 (3)，直到完成所有三角形單元計算。

3 計算方法驗證

3.1 定向運動平板

利用上述方法，對定向運動平板的微流星體撞擊通量進行了計算，并與文獻 [7]進行了比對，其結果如表1—表3所示。表中β為平板運動速度Vs與平板法線的夾角，如圖1所示。Vm為微流星體運動速度。

表1 定向運動平板微流星體通量計算結果，Vs=VmTab.1 The calculation results of micrometeoroid flux for directional movement plates，Vs=Vm

表2 定向運動平板微流星體通量計算結果，Vs=7.6km/s，Vm=16.8km/sTab.2 The calculation results of micrometeoroid flux for directional movement plates，Vs=7.6km/s，Vm=16.8km/s

表3 定向運動平板微流星體通量計算結果，Vm速度分布采用式 (2)，Vs=7.6km/sTab.3 The calculation results of micrometeoroid flux for directional movement plates，the velocity distribution of Vmusing formula(2)，Vs=7.6km/s

可見計算結果符合較好，最大相對誤差僅1.5%，證明本文計算方法能夠正確描述微流星體的速度效應。

3.2 典型結構微流星體通量

針對防護手冊[2]中用于軟件對比驗證的算例進行了計算，并與MDPANTO軟件的結果進行了比較。算例具體參數(shù)此處不再贅述。由于遮擋效應和速度方向效應對所有尺寸的微流星體影響相同，因此本文僅給出了直徑小于0.1mm的微流星體通量計算結果。計算結果討論如下。

圖4 立方體構型衛(wèi)星微流星體通量計算結果Fig.4 The calculation results of cube configuration satellite micrometeoroid flux

立方體構型衛(wèi)星網(wǎng)格劃分如圖4(a)所示。總的網(wǎng)格數(shù)為1562個。計算結果如圖4(b)所示。

與MDPANTO軟件對比的結果如表4所示。

表4 立方體構型衛(wèi)星微流星體通量計算結果對比Tab.4 A comparison of the micrometeoroid flux results of the cube configuration satellite

由表5可見，與MDPANTO比較，立方體工況誤差3%。但背風面 (Back)和對地面(Earth)兩個表面的誤差較大。可能原因在于與MDPANTO采用的地球遮擋和微流星體速度分布抽樣計算方法不同，該問題仍需進一步研究。

球形構型衛(wèi)星的網(wǎng)格劃分結果如圖5(a)所示，共有5832個網(wǎng)格，計算結果為14.277，MDPANTO結果為14.57，誤差2.01%。

空間站構型網(wǎng)格劃分結果如圖6(a)所示，共有9068個網(wǎng)格。計算得到總的通量為89.233，MDPANTO結果為92.47，相對誤差分別為3.5%。

三平行平板遮擋效應算例計算結果如圖7所示。相對與MDPANTO的誤差分別為2.40%、4.34%和4.31%，符合較好。

4 復雜結構航天器微流星體撞擊通量計算

利用上述計算方法對某復雜結構大型航天器微流星體撞擊通量進行了計算，計算結果如圖8所示。軌道高度為400km，坐標系如圖3所示。共劃分表面三角形網(wǎng)格24696個。每個網(wǎng)格射線數(shù)目為10萬條。對于直徑小于0.1mm的微流星體，計算得到的總撞擊通量為256.37。

從圖8中可以看出，由于地球遮擋和各艙段自身的遮擋，圖中所示部位2附近撞擊通量最小，僅為0.305。而圖中所示部位1附近則通量最大，這是由于該處為迎風方向和天頂方向。此處單元最大通量為6.597。因此，本文計算方法能夠正確捕捉表面單元之間復雜的遮擋關系，計算結果符合物理規(guī)律。

圖5 球形構型衛(wèi)星微流星體通量計算結果Fig.5 The calculation results of spherical configuration satellite micrometeoroid flux

圖6 空間站構型微流星體通量計算結果Fig.6 The calculation results of space station configuration micrometeoroid flux

5 總結

本文基于射線追蹤原理提出了微流星體撞擊通量的計算方法。該方法將微流星體以定向的射線表示，將航天器表面劃分表面三角形單元，然后隨機抽樣選取射線的方向及其代表的微流星體速度，判斷其是否撞擊表面，并計算與其他表面元和地球的交點。如果存在交點，則說明存在遮擋。通過計算大量射線的軌跡，并求其平均值，可得微流星體通量及其分布。

針對定向運動的平板進行了計算并與理論解進行了比較，最大誤差1.5%，符合較好。針對IADC防護手冊中給出的立方體、球體、簡化空間站等構型算例和三平板遮擋算例進行了計算，并與MDPANTO的結果進行了對比，總通量誤差分別為3%、2%、3.5%和4.3%，符合較好。通量分布結果符合物理規(guī)律，但不同朝向的計算結果最大誤差達16%，需進一步研究。利用該方法對某復雜構型大型航天器微流星體撞擊通量進行了計算，計算結果符合物理規(guī)律。

圖7 三平行平板遮擋效應計算Fig.7 The calculation results of the shielding effect of three parallel plates

圖8 大型航天器微流星體通量計算結果Fig.8 The calculation results of large spacecraft micrometeoroid flux