小天体高速撞击溅射的接近观测风险分析

焦艺菲，程彬 *，宝音贺西, 2*

(1.清华大学航天航空学院，北京100084；2.内蒙古工业大学，呼和浩特010051)

1 引言

太阳系内已发现的小天体超过数百万颗，其中近地小天体(包括近地小行星、彗星等)超过3万颗，据估计，直径50m以上的近地小天体的真实数量可能超过10万颗，但目前发现的仅有2%，由于轨道与地球接近，这些小天体可能与地球相撞并引发区域级的灾难[1-3]。为了应对近地小行星的潜在威胁，国内外学者提出了众多防御手段，包括动能撞击、核爆炸、长期作用力偏转等，其中动能撞击方案旨在利用航天器的高速撞击来改偏小天体的轨道，进而避免与地球相撞，是目前最有效可行的防御手段之一[4-10]。

在2022年中国航天日大会上，国家航天局表示，我国将着手组建近地小行星防御系统，并计划在2025年左右实施一次对某颗有威胁的小行星既进行抵近观测，又实施就近撞击改变其轨道的技术实验[11]。小行星撞击防御任务需要对高速撞击过程进行尽可能精确的观测，地基望远镜和绕地的观测卫星难以实现，而使用探测器接近观测能够对撞击溅射过程和撞击的直接结果进行记录，但是高速撞击溅射的碎片可能对近距离的探测器造成威胁。此前国外航天机构也开展过多次深空的高速撞击和接近观测任务，如2005年的“深度撞击”(Deep Impact)探测器高速撞击彗星9P/Tempel，采用近距离飞越观测的方式，最小飞越距离为500km，并成功观测到撞击产生的大量发光溅射尘埃，但由于观测时间窗口较短且观测距离较远，成像质量较差且无法观测到最终的撞击坑[12]；2009年的“月坑观测与感知卫星”(LCROSS)月球撞击任务使用火箭上面级撞击月球表面，在撞击后4min内，使用飞船对喷射羽流进行近距离观测并传输数据，随后飞船再次撞击月球[13]；2019年“隼鸟”2(Hayabusa 2)探测器在小行星“龙宫”(Ryugu)表面进行了小型人工撞击试验，在撞击侧向约1km位置放置小型相机进行观测，探测器则逃离到小行星背面约20km处以免被溅射碎片撞击[14]；2022年9月，“双小行星重定向测试”(DART)探测器成功撞击Didymos双小行星系统的子星Dimorphos，携带的立方星LICIACube在撞击前分离，以55km的距离飞越观测撞击溅射情况，但目前发布的撞击图像无法用于精确评估轨道偏转效果[15-17]；此外，计划于2024年发射的Hera任务将再次造访Didymos双星并对DART的撞击结果进行全面评估[18]。

相比于飞越观测的方式，伴飞、绕飞的方式可以对小天体的撞击过程进行长时间、近距离的观测[19,20]，因此能够有效观测到高速撞击溅射、撞击坑形成、低速碎片在不规则的弱引力场中演化的整个过程，既是对撞击防御结果的直接评估，也对研究小天体的材料性质、内部结构和理解小天体的形成与演化过程有重要价值。针对接近观测任务的安全性需求，本文使用完全自主知识产权的小天体高速撞击数值模拟软件THU-SPHSOL，通过数值仿真研究了均质岩石小天体的高速撞击溅射过程，以及不同撞击角度对溅射结果的影响，分析了高速溅射碎片的空间分布情况和观测相机的安全区域，研究结果将为我国即将开展的小行星撞击防御任务提供重要参考。

2 数值仿真模型

小天体的高速撞击和溅射过程涉及强非线性的力-热耦合、流-固耦合，材料在高压高应变率的极端条件下发生大变形、断裂和破碎，传统的基于网格的数值方法难以求解，因此，基于无网格的光滑粒子流体动力学(Smoothed Particle Hydrodynamics, SPH)方法[21]，我们开发了完全自主知识产权的THU-SPHSOL软件[22]，软件面向航天动力学和行星科学等领域的高速撞击问题，能够实现千万粒子量级的大规模高效计算，本节主要对软件采用的物理模型进行介绍。

2.1 拉格朗日描述的控制方程

由于高速撞击产生的冲击压力远远大于材料的强度，材料的行为特性接近无粘的可压缩流体，使用拉格朗日描述的连续介质守恒方程作为控制方程[23]，包括连续性方程、动量方程、能量方程和运动方程：

(1)

(2)

(3)

(4)

2.2 状态方程

状态方程是高速撞击过程中控制材料行为特性的关键部分，描述了材料密度ρ、压力p和比内能e之间的关系，在小天体高速撞击场景中，Tillotson状态方程是最常使用的状态方程之一，在低压/低能量时满足Hugoniot关系，在高压(高能量)时趋于Thomas-Fermi极限[24]。当材料处于压缩(ρ>ρ0)或冷膨胀(ρ<ρ0且e

(5)

当材料处于热膨胀状态时，

(6)

对于中间状态，则根据比内能e进行线性插值，

(7)

状态方程还决定了材料的声速或纵波波速c

(8)

式中：G为剪切模量，对于使用Tillotson状态方程描述的材料，根据式(5)～(8)计算材料声速[25]。

2.3 强度模型

高速撞击的冲击压力随着时间和距离发生衰减，此时压力项和偏应力项将同时主导固体材料的响应特性，因此屈服失效和断裂等与强度相关的特性必须加以考虑，通常按照理想弹塑性模型进行处理。在弹性阶段，假设偏应力率和偏应变率成比例，即满足(增量型)胡克定律，并考虑刚体转动的影响，此时偏应力率的更新形式为：

(9)

材料是否进入塑性阶段则由屈服条件决定，对于小行星等地质材料，采用Lund模型，考虑了内聚力和内摩擦力对屈服强度的影响[26,27]，

(10)

材料的损伤程度描述了材料抗剪切和抗拉伸能力的降低，通常使用一维标量D作为损伤的度量，且0≤D≤1，其中D=0对应完好无损的材料，D=1对应完全失效的材料。对于小行星等脆性地质材料，损伤的形式主要是内部缺陷导致的拉伸损伤，因此采用基于Weibull形式缺陷分布的损伤累积模型[28,29]。根据Grady Kipp模型[30]，材料损伤对抗拉伸能力的削弱通过压力项的修正实现，

p′=(1-D)p,p<0

(11)

损伤对抗剪切能力的削弱则通过修正屈服强度来实现，在Lund模型中，材料完全失效时屈服强度与压力成正比，即Yd=μdp，此时损伤修正的Lund强度为

Y=(1-D)Yi+DYd

(12)

3 仿真结果与讨论

3.1 数值仿真算例

面向我国即将开展的小行星撞击防御任务，典型的撞击场景为质量500kg的航天器以5000m/s的相对速度垂直撞击目标小行星，但考虑到航天器的撞击位置偏差以及小行星的表面地形，真实撞击场景下撞击方向可能偏离撞击表面的法向，撞击溅射分布也可能完全不同，因此分别设置撞击方向偏离表面法向0°、15°和30°三组工况，对应算例(1)～(3)并假设撞击局部的表面平坦，研究不同撞击角度对结果的影响。航天器采用密度2700kg/m3的铝材料，使用Von Mises模型和Tillotson状态方程描述，小行星则采用密度2700kg/m3的玄武岩材料，不考虑孔隙度，使用Lund强度模型、损伤累积模型和Tillotson状态方程描述，材料参数见表1。为提高离散精度并减小计算量，将计算域限制为撞击点周围直径30m的范围，SPH粒子初始间距取0.132m，此时撞击器和小行星分别表示为81个和300万个SPH粒子，每个粒子的质量为6.21kg，撞击仿真总时长取20ms。

表1 材料参数

3.2 结果与讨论

以垂直撞击为例，算例1的数值仿真结果如图1所示，显然，20ms的瞬时撞击坑并不是最终形态，整个成坑过程可能持续几秒到几十秒，但20ms对于研究撞击溅射的较高速碎片是足够的，此时图1中大部分区域都已经完全失效，冲击压力也相比撞击前期大幅衰减，此时粒子的运动速度对后续的演化过程起主导作用，具有向外速度分量的粒子(图中黄色箭头)将继续溅射，而具有向内速度分量的粒子(图中蓝色箭头)在冲击作用下向内压缩，因此不会发生溅射，这也与Z-model的模型预测相符[32,33]。

图1 垂直撞击成坑的溅射过程，箭头方向为粒子速度方向，箭头颜色为沿撞击反向(x轴方向)的速度分量，其中黄色为向外溅射，蓝色为向内压缩；图中标注了撞击后20ms瞬时的撞击坑，以及根据速度场预测的最终撞击坑(直径约16m，深度约6m)

尽管撞击溅射的粒子仍然受到小天体、太阳和其他行星的引力作用，太阳光压、太阳磁场产生的电磁力甚至粒子之间的接触碰撞也都将对粒子的运动产生影响[34,35]，但对于撞击过程的短期(约几十到数百秒)接近观测，仍然可以假设所有溅射粒子都保持撞击仿真结束时的速度做匀速直线运动。我们筛选所有速度大小超过某个阈值且具有向外的速度分量的粒子作为溅射粒子集合，分别取1m/s和10m/s作为速度阈值(远大于百米量级小天体的表面逃逸速度)，为了对比不同撞击角度的溅射结果，根据溅射速度在yz平面的投影，进一步筛选速度投影方向与y轴夹角在±15°以内的粒子进行分析，如图2所示，对于斜撞击情况，筛选区域粒子的溅射速度具有明显的不对称性，并对应了不同撞击角度的溅射分布特征。根据不同的筛选条件，表2统计了溅射粒子的数量分布，可以看出，垂直撞击的溅射粒子最多，随着撞击角度的偏离，溅射粒子数量减少，撞击坑也会更小；速度10m/s以上的粒子数量约为速度1m/s以上粒子数量的6%～7%，因此绝大部分溅射粒子都是低速的，但对于接近观测的探测器，10m/s以上的溅射粒子具有较大的动能且覆盖范围更远，因此相比于低速粒子更具威胁性。

表2 溅射粒子统计

图2 撞击坐标系的定义和溅射粒子的筛选区域

对于溅射速度方向在所筛选区间的粒子，分别考虑1m/s和10m/s两个速度阈值，不同撞击角度(算例1～3)的溅射方向和空间分布如图3所示。可以看出，对于垂直于小行星表面的撞击，速度1m/s以上的溅射粒子与撞击反向的夹角主要分布在20°～70°之间，而速度10m/s以上的溅射粒子角度主要分布在30°～60°，其中分布在40°～50°的粒子数量最多，此前的文献研究和撞击实验[36,37]同样指出，垂直撞击成坑过程中，不同位置发射的溅射物速度与水平方向角度分布在40°～55°之间，但整体的溅射角度接近45°，并使用标度律来描述溅射位置、速度大小和累积质量之间的关系，当然斜撞击情况下溅射角度还会随着方位角发生变化[38,39]；如果只考虑对接近观测造成威胁较大的、速度10m/s以上的溅射粒子，并假设这些粒子在短期的演化过程中做匀速直线运动，在随后的1～100s时间内，溅射粒子的空间位置也基本分布在以撞击点为中心的锥形范围内，锥角区间与速度角度区间基本相同，都是30°～60°，且溅射分布关于撞击面法向(对于垂直撞击也是撞击方向)轴对称。随着撞击方向的偏离，溅射粒子在撞击点两侧的分布呈现出明显的不对称性，下坡侧的溅射粒子数量远多于上坡侧，因此在随后的1～100s时间内，溅射粒子的空间位置关于撞击方向不再旋转对称，但基本上仍是关于撞击面法向的轴对称分布，即使两侧的溅射粒子数量不同。

图3 不同撞击场景的溅射角度分布，其中(a)～(c)为速度超过1m/s的部分，(d)～(f)为速度超过10m/s的部分，极坐标系的横向坐标为溅射速度与撞击反向(x轴方向)的夹角，径向坐标为每个角度区间的粒子数量，红色箭头表示撞击方向，绿色线表示撞击表面；不同撞击场景的溅射空间分布，包括(g)～(i)，只考虑速度超过10m/s的溅射粒子，假设粒子从20ms(撞击仿真的结束时刻)开始做匀速直线运动，极坐标系的横向坐标为粒子的空间位置与x轴反向的夹角，径向坐标为溅射粒子远离撞击点的距离

在真实的小行星撞击防御任务中，航天器的撞击角度不仅取决于撞击点的选取和小行星的形状，还会受到导航和制导精度的影响[40,41]，因此真实撞击角度可能分布在某个范围内。考虑直径50m的球形小行星，假设标称撞击点指向球心，那么最终撞击点可能是以标称撞击点为中心、某个距离范围内的任意位置：如果撞击位置误差在半径6.5m范围内，对应的撞击角度总是小于15°，此时绝大部分撞击溅射粒子分布在与撞击反向(x轴方向)夹角20°～80°的锥形范围内，并且将在100s的时间内到达1～10km外的区域，观测器应避开该角度区间，或者在更远处观测一段时间后重新寻找安全区域；如果撞击位置误差在半径12.5m范围内，对应的撞击角度不超过30°，此时溅射粒子可能分布在与撞击反向夹角0°～90°的半空间范围内，但如果观测相机躲避到撞击的背面，将难以拍摄到撞击溅射过程，此时建议观测相机放置在约10km外的区域，并在约100s后离开该区域。

4 结论

本文结合我国即将开展的小行星撞击防御任务和接近观测的安全性需求，使用小天体高速撞击数值模拟软件THU-SPHSOL对撞击过程进行数值仿真研究，介绍了软件主要采用的物理模型和材料参数，数值仿真研究了均质岩石小天体的高速撞击溅射过程，以及不同撞击角度对溅射结果的影响。研究表明，垂直撞击的高速溅射粒子主要分布在与撞击反向夹角为30°～60°的锥形范围内，但撞击方向的偏离会引起溅射粒子数量和角度分布的不对称性，在真实撞击场景下，观测相机的安全区域范围还将取决于撞击位置的误差范围。本文的研究结果将为我国即将开展的小行星撞击防御任务提供重要参考。