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

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

1 概述

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

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

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

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

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

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

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

速率單位為km/s。

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

2.2 射線追蹤計(jì)算方法

如圖1所示為微流星體撞擊平板時(shí)通量計(jì)算的示意圖。平板法線方向?yàn)閚，運(yùn)動(dòng)速度為vS，微流星體速度為vm。微流星體撞擊平板的速度為vi，vi=vm-vS。

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

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

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

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

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

對(duì)應(yīng)圖2，計(jì)算流程如下:

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

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

圖3 微流星體通量計(jì)算坐標(biāo)系Fig.3 Coordinate system for calculating flux of micrometeoroids

(3)針對(duì)每個(gè)網(wǎng)格單元:

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

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

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

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

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

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

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

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

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

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

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

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

3 計(jì)算方法驗(yàn)證

3.1 定向運(yùn)動(dòng)平板

利用上述方法，對(duì)定向運(yùn)動(dòng)平板的微流星體撞擊通量進(jìn)行了計(jì)算，并與文獻(xiàn) [7]進(jìn)行了比對(duì)，其結(jié)果如表1—表3所示。表中β為平板運(yùn)動(dòng)速度Vs與平板法線的夾角，如圖1所示。Vm為微流星體運(yùn)動(dòng)速度。

表1 定向運(yùn)動(dòng)平板微流星體通量計(jì)算結(jié)果，Vs=VmTab.1 The calculation results of micrometeoroid flux for directional movement plates，Vs=Vm

表2 定向運(yùn)動(dòng)平板微流星體通量計(jì)算結(jié)果，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 定向運(yùn)動(dòng)平板微流星體通量計(jì)算結(jié)果，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

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

3.2 典型結(jié)構(gòu)微流星體通量

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

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

立方體構(gòu)型衛(wèi)星網(wǎng)格劃分如圖4(a)所示?？偟木W(wǎng)格數(shù)為1562個(gè)。計(jì)算結(jié)果如圖4(b)所示。

與MDPANTO軟件對(duì)比的結(jié)果如表4所示。

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

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

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

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

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

4 復(fù)雜結(jié)構(gòu)航天器微流星體撞擊通量計(jì)算

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

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

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

圖6 空間站構(gòu)型微流星體通量計(jì)算結(jié)果Fig.6 The calculation results of space station configuration micrometeoroid flux

5 總結(jié)

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

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

圖7 三平行平板遮擋效應(yīng)計(jì)算Fig.7 The calculation results of the shielding effect of three parallel plates

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