程 彬，宝音贺西

(清华大学 航天航空学院·北京·100084)

0 引 言

人类自古以来就对浩瀚的宇宙充满好奇，渴望向未知领域探索求知。自20世纪50年代以来，人类探索宇宙的脚步不断深入，深空探测项目在近几十年来已成为各航天大国的重要发展方向之一，而小行星采样探测更是其中不可或缺的一部分。小行星的表面保留着太阳系形成初期的原始成分，是研究太阳系起源和演化历史的“活化石”，其表面的撞击坑更是研究太阳系爆炸时期环境变化的最佳材料[1]；同时，其上有丰富的稀有矿物资源。例如，某颗直径仅为500m的小行星蕴含的铂元素就相当于迄今为止人类铂开采量的总和，而其表面的水资源还可作为人类深空探测旅行的补给站[2]；另外，部分近地小行星有撞击地球的可能，通过深空探测了解其结构组成，对于未来可能的小行星防御计划的成功具有非常重要的作用[3]。

目前，各航天大国先后提出并实施了小行星探测规划，以占据航天技术战略制高点，从而在未来的太空资源争夺中保持优势。日本宇航局(JAXA)于2003年发射了隼鸟号(Hayabusa)，在近10年间针对S-类小行星Itokawa进行了近距离抵近探测与采样返回，成为了世界上第一个对小行星成功进行采样返回的国家[4]，而隼鸟二号(Hayabusa II)也已经于2018年6月到达了小行星Ryugu，并成功实施了表面巡游与撞击采样。美国宇航局(NASA)于2016年9月发射了OSIRIS-REx深空探测器，探访了富含碳质生命元素的C-类小行星Bennu[5]。OSIRIS-REx深空探测器在2018年10月到达了小行星Bennu，将在2019年实施表面附着采样。欧洲航天局(ESA)计划于2022年实施Hera小行星探测任务。我国小行星环绕与采样返回任务也正在规划论证中，国内相关研究部门已经积极开展了针对小行星探测的相关研究工作。

考虑到小行星表面引力非常低，且表面覆盖着颗粒状风化层，不易于实现探测器的着陆、锚定，因此在微重力条件下，通常采用与传统采样方案不同的策略，即接触-分离(Touch-And-Go)式采样，利用探测器接触小行星表层的短暂时间，通过收集溅射物质实施采样。接触时间往往仅有几秒，之后需立刻启动反冲发动机以离开小行星表面，从而不需要着陆、锚定机构。如探访小行星Itokawa的隼鸟号，在采样装置接触到小行星表层时发射高速射弹，撞击其浮土层，向上飞溅的岩石碎片在挡板的作用下反弹入返回舱，整个收集过程仅持续了1s，并成功收集到了1534个风化层颗粒[6]。对于小行星采样方案而言，最重要的是收集效率。在隼鸟号发射前，JAXA曾进行了大量实验，以研究射弹形状、射弹速度等参数的影响[7]。但受限于实验条件，JAXA并没有进行更接近小行星真实环境的超微重力实验。得益于数值方法的特性，我们可以更为方便地模拟小行星的微重力环境，研究更为广泛的参数空间的采样过程，同时跟踪风化层每个颗粒的运动，更深入地理解溅射过程，从而优化采样收集效率，为高效采样仪器的设计提供理论基础。针对这一问题，Wada等[8]首先将软球离散元法引入小行星表面风化层的撞击仿真中，研究了撞击导致的碎片速度、角度分布及撞击坑尺寸与撞击速度的关系；Nakamura等[9]在此工作的基础上进行了微重力实验，研究了射弹在低速穿过颗粒层时的速度衰减过程，发现离散元法具有较高的准确性；Schwartz等[10]首先使用离散元法进行了采样过程仿真，研究了溅射质量与射弹形状的关系，并将其与地球环境下的实验进行了对比，发现离散元方法可以较为准确地复现实验结果。但是，由于微重力环境下的仿真耗时较长，Schwartz等并没有进行接近小行星环境的采样过程仿真，也没有考虑采样器的复杂几何结构。

本文采用数值方法研究了小行星表面的撞击采样过程，考虑了小行星的微重力环境与颗粒风化层等特殊地质条件，基于初步的溅射采样装置原型机，研究了射弹形状对撞击坑形状的影响，以及其与采样效率的依赖关系，为我国未来的小行星采样装置设计提供了理论基础。

1 数值模型

现有天文观测与采样结果表明，大部分小行星表面覆盖着颗粒状的风化层，在微重力环境下其动力学行为复杂。目前针对颗粒物质的研究主要采用离散单元法，课题组基于此方法开发了三维大规模颗粒模拟软件DEMBody[11-12]。该软件以牛顿第二定律为基础，所研究的对象是大量具有特定形状的离散颗粒，其采用颗粒领域通用的软球模型来描述散体的动力学行为，如图1所示。其包括弹簧阻尼、Hertz接触力、Cundall切向力等模型，根据粒子交叠的程度刻画颗粒之间的正向力和摩擦力，即有

FN=-knxn+Cnun
FT=min{μ|FN|t;ktS+Ct|ut|t}

(1)

式(1)中，FN和FT分别为颗粒所受到的正向弹性力与摩擦力；kn和kt为弹性系数，Cn和Ct为阻尼系数，与颗粒碰撞恢复系数相关，在颗粒系统仿真中一般取经验值0.5；无量纲参数μ为颗粒间的摩擦系数；x为颗粒的正向侵入量，S为接触过程中颗粒的切向形变矢量；n和t分别为正向与切向的单位矢量；正向力与摩擦力中均包含阻尼力项，其分别与正向和切向相对速度(即un和ut)成正比，用于描述颗粒接触过程中的黏滞耗散行为。

在之前的颗粒系统仿真过程中，为了计算简便，往往会忽略滚动摩擦力矩的作用。已有研究[13]表明，这种简化将会导致颗粒系统行为与真实行为出现较大差异。因此，引入了滚动摩擦系数来代表该动能损耗。当滚动摩擦系数不为0时，将会抑制颗粒对在接触点的滚动相对速度，从而阻碍颗粒对的相对滚动，即有

MR=-μr|FN|Riωi

(2)

式(2)中，MR表示颗粒所受到的滚动力矩；Ri为颗粒半径，ωi为颗粒角速度；无量纲参数μr为颗粒间的滚动摩擦系数。

(a)颗粒接触正向模型

(b)颗粒接触切向模型

(c)颗粒接触滚动模型图1 颗粒接触的离散元模型Fig.1 The contact model in Discrete Element Model simulation

为了模拟数百万尺度的颗粒动力学问题，采用并行计算的方式解决大规模计算问题。DEMBody采用OpenMP与MPI进行多核并行，可在天河二号超级计算机上进行大规模计算，在实际测试过程中可数百倍地提高效率，基本满足本研究的需求。

2 数值程序测试

本部分通过与典型的颗粒动力学实验进行对比，确定了程序的正确性与适用性。LI等[14]通过实验，研究了在飞船与天体表面接触时，产生溅射物质的质量与撞击速度的关系。通过仿真重现了该实验的过程，即使用直径为19.71mm的玻璃射弹撞击由丙烯酸颗粒填充的颗粒床，撞击速度范围为0.9m/s～3.6m/s，颗粒床的厚度范围为46.1mm～138.3mm，直径保持为75.6mm。通过记录撞击前后颗粒床的质量之差，可确定由撞击引起的溅射质量。

图2所示为仿真与实验结果的比较，其中圆点表示仿真结果，叉号表示实验结果，误差线表示由多次重复实验得到的溅射质量范围；蓝色数据表示46.1mm颗粒床，黄色数据表示92.2mm颗粒床，红色数据表示138.3mm颗粒床。仿真结果表明，溅射质量随着撞击速度的增加而逐渐增加，这是由于撞击动能的增加导致了颗粒床所获得的动能增加；同时，溅射质量随着颗粒床厚度的增加而逐渐减小，这是由于颗粒床厚度的增加使得颗粒系统的摩擦能量损耗增加，从而导致颗粒所获得的动能减小。值得注意的是，在每一组仿真参数中，仿真得到的溅射质量都基本在实验误差范围内，这表明程序可以很好地复现颗粒床的低速撞击过程，具有较高的准确性。

图2 仿真溅射质量与实验结果的对比Fig.2 Comparison of ejected mass between simulation and experiment

3 仿真试验及结果分析

3.1 仿真设置

日本隼鸟号对小行星Itokawa风化层的观测表明，小行星风化层的力学性质接近碎石，颗粒半径满足0.2mm～0.25mm的立方幂律分布。由于隼鸟号的探测数据是目前小行星风化层最为精细的数据，因此仿真中采用的类似S-类小行星风化层的参数，如表1所示。

对于射弹，采用了与隼鸟号类似的SUS304铁的物理参数，即密度为8.0g/cm3；同时，为了研究射弹形状对溅射采样过程的影响，采用了相同质量不同形状的射弹，包括球形、60°锥形、90°锥形、120°锥形、半球形、1∶2圆柱、1∶1圆柱、2∶1圆柱，具体几何参数如表2所示。其中，R表示射弹的半径，H表示射弹的总厚度，M表示射弹的质量。

表1 溅射采样仿真中所使用的接触参数

表2 溅射采样仿真中所使用的射弹形状

对于采样装置，采用了与隼鸟号类似的结构，包括可伸缩采样桶(Extendable Sampler Horn)、锥形角(Conical Horn)、捕集器(Sampler Catcher)与收集仓(Container)，分别使用了圆柱、圆锥、圆柱、平板边界模拟。如图3为采样装置的尺寸说明。

图3 采样装置、部分射弹形状及采样原理示意图Fig.3 Sampling equipment, several projectile shapes and sampling principle

在仿真中，首先将13万颗粒随机撒入直径为28cm的圆柱容器中，在小行星重力作用下形成模拟风化层。随后，将采样装置放置于颗粒层上方，同时将射弹以25m/s的速度垂直向下射出，设定采样时间为1s。

3.2 仿真结果

不同形状的射弹侵入风化层后，均在其后形成不断扩大的撞击空腔，同时撞击物将携带的动能通过力链的方式传递给颗粒层，导致部分颗粒从风化层表面溅起，形成了倒圆锥状的溅射云。如图3(右)所示，这些溅射物质将撞击采样桶，随后反弹并沿采样桶不断上升，直至到达收集仓。需要指出的是，由于采样时间设置为1s，即若颗粒在1s内到达收集仓，才能认为该颗粒被成功采样。

如图4所示，对于头部尖锐的射弹，其形成的撞击坑近似为锥形，撞击坑更深，但宽度小，如60°锥形、1：2圆柱形；对于头部圆钝的射弹，其形成的撞击坑近似为U形，撞击坑更浅，但宽度大，如120°锥形、2：1圆柱形；对于介于其中的射弹，其形成的撞击坑近似为碗形，撞击坑的深浅与宽度也介于两者之间，如90°锥形。其原因在于，不同形状射弹的截面积与接触风化层面积均不同。侵入时，所受接触力不同，对风化层的冲击作用也不同，60°锥形射弹所受阻力较小，侵入深度大，易形成锥形撞击坑；120°锥形射弹所受阻力较大，侵入深度小，易形成U形撞击坑；90°锥形射弹所受的阻力介于两者之间，因而易形成碗形撞击坑。

图4 不同形状射弹撞击得到的撞击空腔。从上至下、从左至右分别为球形、60°锥形、90°锥形、120°锥形、半球形、1∶2圆柱、1∶1圆柱、2∶1圆柱Fig.4 The impact cavity generated by projectiles with different shapes. From left to right then top to bottom: sphere, 60° cone, 90° cone, 120° cone, hemisphere, 1∶2 cylinder, 1∶1 cylinder, 2∶1 cylinder

不同形状射弹在撞击后所形成的撞击坑形状、颗粒速度分布、溅射区域均有不同，从而影响了收集仓采样的效率。如图5所示，对于圆柱形射弹而言，其采样质量与射弹的尖锐程度近似成反比，即射弹越尖锐(如1：2圆柱形射弹)，采样质量越小；但对于锥形射弹而言，其采样质量与射弹的尖锐程度近似成正比，即射弹越圆钝(如120°锥形射弹)，采样质量越小。在所有被研究的射弹形状中，90°锥形射弹的采样质量最高，60°锥形射弹相比略低。因此，我国未来的小行星采样任务可以将90°锥形射弹作为首选。

图5 不同形状射弹撞击收集到的风化层质量 Fig.5 The sampling mass collected by projectiles with different shapes

4 结 论

本文将散体动力学仿真方法应用到了适合于微重力小行星环境的射弹溅射采样优化过程中，重点研究了同等质量的不同形状射弹对弹坑形状、采样效率的影响，涉及到了不同角度的圆锥形、半球形、圆柱形等形状，比较得出的结论是90°锥形射弹兼顾了提高溅射速度与溅射质量两方面的因素，采样效率最高。因此，我国在未来的小行星探测器中可首先考虑采用90°锥形射弹。但由于本文采用的离散元法未能考虑颗粒破碎、热量传递能问题，在后续的研究中需结合光滑粒子动力学法(SPH)等冲击动力学模型对这一问题继续进行研究。