基于射线追踪原理的微流星体撞击通量计算研究

2019-09-16焦子龙姜利祥李涛孙继鹏黄建国朱云飞

空间碎片研究 2019年2期

焦子龙,姜利祥,李涛,孙继鹏,黄建国,朱云飞

(1.可靠性与环境工程技术重点实验室,北京100094；2.北京卫星环境工程研究所,北京100094)

1 概述

微流星体一般是指直径小于1mm或质量小于1g的流星体，主要来自于小行星和彗星。微流星体在地球附近、行星际空间、星际空间等广泛存在。微流星体撞击航天器可能带来设备功能失效等多种影响。为确保航天器的微流星体撞击防护措施有效、适当，需要首先定量评估微流星体的撞击风险。早在20世纪60年代的阿波罗任务和礼炮号空间站任务中就开始重视并评估微流星体撞击影响[1]。目前已有 BUMPER、ESABASE/DEBRIS、 COLLO、 BUFFER、 MDPANTO、 MODAOST、ARMOR等国内外多种空间碎片/微流星体撞击风险评估软件[2-4]。

本文提出了一种基于射线追踪计算微流星体撞击通量的方法。该方法原理简单、易于实现，便于与其他空间环境因素效应分析模块集成[5]。本文首先介绍了基于射线追踪的微流星体撞击通量计算方法，然后通过算例对该方法进行了对比验证分析，最后对某空间站复杂构型进行了微流星体撞击通量计算。

2 微流星体撞击通量计算方法

2.1 微流星体环境模型

本文采用的微流星体环境模型为Grün模型，它定义了距太阳一个天文单位 (1AU)处相对静止的随机翻转的面积为1㎡的平板在各向同性的微流星体流中一年时间内受到的撞击数目F(m):

由于各环境模型中使用的均是微流星体的尺寸下限，所以它对应的通量值表示大于该尺寸的所有粒子的通量总和。

微流星体相对于地球的速率分布函数为:

速率单位为km/s。

由于航天器的运动使得其不同朝向的表面微流星撞击通量显著不同，例如迎风朝向的表面撞击通量较尾流方向的通量高出将近一个数量级。此外，航天器复杂的结构使得不同表面之间存在遮挡，造成微流星通量显著不同。因此不能直接采用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)完成初始化工作，包括航天器表面划分三角形网格，设置每个网格单元格需要追踪的射线数目等参数。

(2)设置航天器轨道和姿态参数，建立计算坐标系。本文中计算坐标系如图3所示。航天器位于+Z轴。

图3 微流星体通量计算坐标系Fig.3 Coordinate system for calculating flux of micrometeoroids

(3)针对每个网格单元:

1)生成新的射线。在立体角为4π的空间内均匀随机产生射线方向，其方法如下:

式中，r0和r1为 [0，1]区间的均匀分布随机数；θ和φ分别为天顶角和方位角；x，y，z为笛卡尔坐标系下的射线方向。

2)根据微流星体环境模型随机抽样微流星体速率值。随机抽样可采用舍选抽样法。

3)计算微流星体相对撞击速度vi。如果vi·n＜0，则微流星体撞击该表面单元，继续d。否则判断射线抽样数目是否足够。如果不够，继续执行a。

4)抽样射线的初始位置，其方法如下:设三角形单元三个顶点分别为A、B和C，其在系统坐标系中的坐标向量分别为，则射线初始位置为:

式中，r0和r1为 [0，1]区间的均匀分布随机数。若r0+r1＞1，则r0=1-r0，r1=1-r1。

5)判断射线是否与其他三角形单元相交。为加速计算是否存在交点，可采用八叉树等空间划分结构[6]。如果相交，判断射线抽样数目是否足够。如果不够，返回步骤1)。如果不相交，继续步骤6)。

6)根据式 (3)计算通量增量ΔFi。三角形单元受到的总通量为:

式中，N为追踪的射线数目；F0为根据式 (1)计算的微流星体通量；G为引力会聚因子。对于地球，表达式为:

式中，RE为地球半径；h为轨道高度。

7)循环步骤1) ～6)，直到达到指定的射线数目。

(4)循环步骤 (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软件的结果进行了比较。算例具体参数此处不再赘述。由于遮挡效应和速度方向效应对所有尺寸的微流星体影响相同，因此本文仅给出了直径小于0.1mm的微流星体通量计算结果。计算结果讨论如下。

图4 立方体构型卫星微流星体通量计算结果Fig.4 The calculation results of cube configuration satellite micrometeoroid flux

立方体构型卫星网格划分如图4(a)所示。总的网格数为1562个。计算结果如图4(b)所示。

与MDPANTO软件对比的结果如表4所示。

表4 立方体构型卫星微流星体通量计算结果对比Tab.4 A comparison of the micrometeoroid flux results of the cube configuration satellite

由表5可见，与MDPANTO比较，立方体工况误差3%。但背风面 (Back)和对地面(Earth)两个表面的误差较大。可能原因在于与MDPANTO采用的地球遮挡和微流星体速度分布抽样计算方法不同，该问题仍需进一步研究。

球形构型卫星的网格划分结果如图5(a)所示，共有5832个网格，计算结果为14.277，MDPANTO结果为14.57，误差2.01%。

空间站构型网格划分结果如图6(a)所示，共有9068个网格。计算得到总的通量为89.233，MDPANTO结果为92.47，相对误差分别为3.5%。

三平行平板遮挡效应算例计算结果如图7所示。相对与MDPANTO的误差分别为2.40%、4.34%和4.31%，符合较好。

4 复杂结构航天器微流星体撞击通量计算

利用上述计算方法对某复杂结构大型航天器微流星体撞击通量进行了计算，计算结果如图8所示。轨道高度为400km，坐标系如图3所示。共划分表面三角形网格24696个。每个网格射线数目为10万条。对于直径小于0.1mm的微流星体，计算得到的总撞击通量为256.37。

从图8中可以看出，由于地球遮挡和各舱段自身的遮挡，图中所示部位2附近撞击通量最小，仅为0.305。而图中所示部位1附近则通量最大，这是由于该处为迎风方向和天顶方向。此处单元最大通量为6.597。因此，本文计算方法能够正确捕捉表面单元之间复杂的遮挡关系，计算结果符合物理规律。

图5 球形构型卫星微流星体通量计算结果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%，符合较好。通量分布结果符合物理规律，但不同朝向的计算结果最大误差达16%，需进一步研究。利用该方法对某复杂构型大型航天器微流星体撞击通量进行了计算，计算结果符合物理规律。