&nbsh1;

 

OTM/HOTM极限力学仿真在小行星防御中的应用

廖祜明，龚自正，宋光明，樊宗岳，杨宏涛，黎波*

(1.北京航空航天大学能源与动力工程学院，北京 102206；2.北京卫星环境工程研究所，北京 100094；3.北京大学工学院，北京 100871；4.云翼超算(北京)软件科技有限公司，北京 100027)

1 引言

地球自诞生以来就不断受到小行星(Asteroid)的撞击，这是人类生存面临的重大潜在威胁之一。如何减缓或避免小行星撞击地球已经成为学术界和国际社会关注的热点问题[1]。在天文学上，定义轨道在距离太阳1.3AU(1AU=1.496×108km)范围内，距离地球轨道最小距离在0.3AU范围内的小行星为近地小行星(Near-Earth Asteroid, NEA)，截至2022年9月26日，已发现NEA 29901颗[2]。当小行星与地球距离0.05AU时，有可能被地球引力俘获，改变运行轨道与地球相撞，而直径140m以上的小行星，撞击地球的威力足以造成区域性灾难，毁掉一个中等大小的国家，因此把距离地球轨道最小距离在0.05AU范围内、直径大于140m的小行星定义为具有潜在碰撞威胁的小行星(Potentially Hazardous Asteroid, PHA)[3]。目前已发现的PHA有2289颗[2]，是防御任务的重点目标。

在能够对PHA提前预警的前提下，根据预警时间长短以及目标小行星尺寸的不同，将小行星分裂成碎片或者改变小行星轨道是避免其撞击地球的两种基本方式。具体来说：一是对于预警时间充裕且尺寸较小的PHA，可采用长期作用力来改变小行星轨道(包括太空拖船、引力拖车、质量推进器、激光烧蚀、表面喷漆、离子束等)；二是对于尺寸较小、预警时间较短或尺寸较大且预警时间较长的PHA，可采用动能撞击技术直接改变小行星运行轨道；而对于预警时间很短且直径大于500m的PHA，目前研究认为摧毁或部分摧毁小行星的核爆炸(表面爆炸、对峙爆炸以及穿透爆炸)是进行防御的有效手段[4-8]。

2022年9月27日，在太空经历了长达10个月的飞行后，执行美国国家航空航天局(NASA)“双小行星重定向测试”(DART)任务的卫星成功撞击了直径约163m的Dimorphos小行星，这是全球首次行星防御技术的演示，也是NASA首次尝试在太空中对小行星的运行进行偏转，作为NASA整体行星防御战略的一部分。DART对小行星Dimorphos的撞击展示了一种可行的防御技术，观测DART任务撞击效果的“小行星撞击监视器”(HERA)与“小行星撞击偏转评估 计划 项目进度计划表范例计划下载计划下载计划下载课程教学计划下载 ”(AIDA)[9]是正在实施的针对动能撞击技术的在轨验证试验，将进一步评估动能撞击对于Dimorphos公转速度和轨道的影响[4]。预先计算显示此次撞击将给Dimorphos的运行速度造成0.2mm/s的改变，在长时间的飞行过程中其轨迹将达到一个可观的累积偏转量。由于试验研究存在周期长、成本高与测量技术受限等困难与挑战，数值计算成为了辅助预测防御措施有效性的高效手段之一。

在动能撞击数值仿真领域，张韵[10]等采用物质点法(Material Point Method，MPM)进行了铝弹超高速撞击碎石堆结构(rubble-pile structure)和完整结构(monolithic structure)两种不同小行星 材料 关于××同志的政审材料调查表环保先进个人材料国家普通话测试材料农民专业合作社注销四查四问剖析材料 模型的动能撞击效果评估；Jutzi[11]等提出了一种改进的SPH 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 ，模拟了小行星灾难性破坏的过程，结果展示了材料屈服极限、剪切强度和孔隙率在碰撞中的相对重要性；姜宇[12]等将小行星Bennu视作由球型颗粒聚成的碎石堆，采用软球离散元耦合引力N体模型对其进行了碰撞数值模拟；Raducan[13-15]等以DART任务作为研究案例，使用基于ALE方法的iSALE软件，对类DART卫星撞击小行星作了大量的数值模拟，总结了卫星几何形状、撞击速度与角度和小行星材料属性等参数对撞击坑坑径和撞击动量转换效率的影响；马鑫[16]等采用High-explosive-burn材料模型和JWL状态方程描述聚能爆炸成型弹丸(EFP)，采用Johnson-Cook材料模型和Gruneisen状态方程描述混凝土球体靶板，利用有限元方法进行深空撞击载荷仿真，形成了撞击速度、靶板强度、靶板密度、靶板体积不同参数条件下的撞击坑坑径变化规律。在核爆仿真领域，李毅[17]等在欧拉型冲击动力学仿真软件NTS中加入能量源，模拟了核爆装置在不同深度爆炸对小行星产生的偏转与破坏效应；汤文辉[18]等采用有限元方法以及能够描述气化反冲的PUFF物态方程进行了核爆炸偏转小行星数值模拟；Patrick K[19]等使用Spheral++[20]软件中的自适应SPH算法和N体引力算法模拟了核爆过程和碎片重组过程，研究结果表明在实际操作中需要合理预估核爆冲击能量，权衡风险系数。可见，国内外研究人员在小行星防御的数值仿真领域已开展相关研究工作并取得了一定的成果。

然而，由于小行星防御及与小行星撞击地球涉及的极高速、极高温以及强耦合的极端苛刻环境，如何高效高精度预测材料在高温、高压、高应变率等极限热-力-化学条件下的动态响应机制一直是数值仿真模拟的难点之一。目前围绕小行星防御的数值计算在基础理论、计算方法、材料模型与高性能计算方面都还面临许多困难与挑战。基础理论方面，由于在核爆、动能撞击、激光驱动等将小行星分裂成碎片或改变轨道的防御措施中，以及小行星进入地球大气与撞击地球表面的过程中，大量存在着烧蚀、解体、爆炸、火球、撞击成坑、反溅碎片云、地震等一系列复杂的物理化学和力学现象，需要理论模型能够综合考虑其中的热流固耦合、热力化学耦合等多物理场强耦合效应；计算方法方面，目前开展的工作主要采用基于网格和无网格方法两大类。基于网格的数值方法中，传统拉氏有限元法受网格畸变困扰，难以处理超大变形问题；欧氏方法具有网格不变形的固有优势，但如何精确模拟与加载速率和加载路径相关的材料响应、准确跟踪物质界面和固液气动态相变仍然是亟待解决的问题，同时在计算动态裂纹扩展与多体碰撞等非线性行为时还需进一步研究；任意拉格朗日-欧拉(Arbitrary Lagrangian Eulerian, ALE)[21-23]法在一定程度上缓解了网格畸变带来的困难，但ALE需要复杂的网格管理和映射技术，同时也面临超大变形时网格畸变的问题[24]。无网格法由于不需要进行网格离散及采用高阶插值函数，在大变形问题上应用更为广泛。常用的无网格方法有光滑粒子流体动力学法(Smoothed Particle Hydrodynamics, SPH)[25,26]、再生核粒子法(Reproducing Kernel Particle Method, RKPM)[27]、软球离散元法(Soft-sphere Discrete Element)[28]以及物质点法(Material Point Method, MPM)[29]等，这些方法在超高速动能毁伤和空间碎片防护研究中发挥了巨大作用，但在对多物理场强耦合的处理、稳定性和效率等方面仍然存在不足[30]。材料模型方面，与航天飞行器经过人类精心 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 的材料和结构不同，小行星的材料和结构是自然形成的，不仅呈现出各向异性的特征，而且充满孔隙和裂纹，因而对典型材料模型和高状态相变数据都提出了更多的需求[1]。

对小行星防御场景的高效高精度仿真是一个对基础理论、计算方法、材料模型不断丰富完善的过程，其对数值计算方法的基础理论框架提出了具有通用性和可扩展性的需求，从而实现不断融合新的物理化学现象、新的材料模型和更加真实的边界条件，形成小行星防御专用仿真平台，更加真实高效地仿真各种小行星防御场景。最优运输无网格方法(Optimal Transportation Meshfree, OTM)[31]是一种针对动态冲击问题提出的显式增量更新拉格朗日无网格方法。OTM采用基于变分原理的多物理场自主强耦合理论框架，为了求解热流固耦合问题，在OTM基础上通过引入能量守恒方程和热流固耦合变分本构，扩展形成热力强耦合最优运输无网格方法(Hot Optimal Transpor-tation Meshfree, HOTM)[32]，进一步地通过结合分布式多进程并行与共享内存多线程并行策略，实现了大规模双层混合并行最优运输无网格方法(massively parallel OTM method, pOTM)[33]，显著提高了计算效率。在材料失效方面，在OTM系列方法的基础上，采用基于能量准则的裂纹扩展算法(EigenFracture)[34]，以材料能量释放率作为失效判据，克服传统裂纹扩展数值方法裂纹扩展路径网格相关、不收敛、无法清晰表征材料真实变形和失效物理机制的缺点。通过有机融合这些方法形成的ESCAAS高度非线性多物理场强耦合极限力学仿真平台，如图 1所示，实现在统一的框架下求解材料超大变形、熔化、气化、冲击爆炸、流固与热流固耦合、自由表面、多体接触、碎裂、层裂与碎片云等复杂物理现象的多物理场、多尺度强耦合问题，为小行星防御仿真预测提供了具有潜力的解决方案。

图1 基于OTM的极限力学仿真理论框架

本文将对OTM/HOTM极限力学仿真理论框架进行介绍，并探讨其在陨石超高速撞击成坑、动能撞击、激光烧蚀驱动、核爆拦截等偏移小行星轨道以及摧毁小行星等防御措施中的潜在应用方向和相关案例，为小行星防御提供相关理论与技术支撑。

2 OTM/HOTM极限力学仿真理论框架

本节简要表述极限力学问题的理论基础与基于OTM/HOTM方法的数值求解框架。

2.1 拉格朗日描述控制方程

不失一般化，以一般带耗散热力耦合问题中Ω0⊂R3的连续物体进行讨论，其运动过程可表述为随时间变化的变形映射φ:Ω0×[t0,t]→R3，其中Ω0代表参考构型，[t0,t]为运动过程总时间。X∈Ω0表示材料质点在参考构形中的位置，x=φ(X,t)表示材料质点在现时构形Ωt=φ(Ω0,t)的位置，则拉格朗日描述下质量守恒、动量守恒与能量守恒方程为：

ρ(x)J=ρ0

(1)

(2)

(3)

2.2 多物理场强耦合变分框架

结合热力耦合变分原理[35]，建立与运动控制方程、初始条件与边界条件等效的变分结构。在连续介质力学中，与速率问题相对应的变分结构可表示为：

(4)

(5)

式中：A表示Helmholtz自由能，φ*为粘性能量耗散，ψ*为塑性变形能量耗散，χ为与热传导相关的能量耗散。矢量G=-(∇)0T)/T描述热能状态，∇0表示在参考坐标系下的偏导。Θ(F,T,Z)=∂U/∂N表示材料内部温度场由内能推导得到，并在热平衡状态下满足Θ=T。可见等效势能由系统存储的能量(Helmholtz自由能)和系统耗散的能量(熵、粘性耗散、塑性耗散、热传导耗散)所组成。其中，Helmholtz自由能可分解为：

(6)

2.3 最优运输无网格数值离散

在有限元法中，每个单元由若干节点组成特定形状，各节点与单元有固定的从属连接关系，每个单元携带了一定量的材料质量并具有相应的体积，材料的响应在单元积分点上求解，插值则在单元节点上进行。与有限元法类似，OTM方法中采用物质点xp与节点xa相结合的方式进行空间离散，如图 2所示，其中下标a表示节点索引，p表示物质点索引，下同。OTM中的节点即为有限元单元中的节点，而物质点则为有限元单元中的积分点，不同于有限元的是OTM中节点和物质点不再有固定的连接关系，物质点与节点的联系将在每一步计算过程中动态构建，即此“单元”没有固定的形状，其形状将根据物质的变形可以是任意的形状，由于解放了网格的束缚，因此克服了拉氏方法在处理大变形问题时因网格畸变带来的困难[36]。

图2 空间离散示意图[36]

材料的动力学信息存储在节点上，包括位移、速度、加速度与温度等。tk+1时刻的节点的位移、速度及当前的温度场更新公式为：

(7)

(8)

(9)

从tk时刻到tk+1时刻的位移传输(变形)映射φk→k+1为：

φk→k+1=∑a∈NH(xp,k)xa,k+1Na(xp,k)

(10)

材料的物理信息存储在物质点上，包括质量、体积、密度、变形、应力与材料内部参数等。物质点的运动通过节点的动力学信息插值获得。在计算过程中，物质点在tk+1的位置和温度更新，通过对邻域内的节点进行插值获得：

xp,k+1=∑a∈NH(xp,k)xa,k+1Na(xp,k)

(11)

(12)

式中：Na(xp,k)为局部最大熵无网格插值函数(Local Maximum Entropy, LME)[37,38]，NH(xp,k)代表物质点xp的邻域。LME插值函数满足严格的非负性，以及0阶和1阶连续性要求，同时在边界满足Kronecker-Delta属性。此外，在LME插值函数中，通过调整γ值，插值函数作用范围可从局部有限元无缝过渡到全局无网格，如图 3所示，这种衰减特性建立了与高斯径向基函数的联系，其重要作用是只有少量的节点对待求函数有明显的贡献，从而大大降低计算成本。

图3 局部最大熵插值函数[38]

通过将上述离散的参数代入半离散变分公式中，并实施驻点条件，可以获得完全离散的应力平衡方程和热平衡方程：

(13)

(14)

2.4 基于物理的裂纹扩展算法

材料在极限条件下的裂纹扩展过程包含了多种能量传播与耗散方式的结合与竞争，准确地描述这些特征需要采用基于物理的裂纹扩展算法。OTM方法中采用基于本征侵蚀的EigenFracture裂纹扩展方法，在该方法中以能量的方式来描述整个材料的响应过程，每个物质点代表一小块物质，允许物质点失效，物质点失效代表该小块物质内部产生裂纹并形成自由表面，而自由表面的产生将伴随能量耗散，裂纹路径则由各种能量耗散的竞争与耦合共同决定。在EigenFracture中物质点失效的判断准则为：

(15)

图4 黑色点为一组失效物质点组成的裂纹，圆圈内的红点为位于裂纹尖端的物质点的邻域η

由于EigenFracture算法采用材料固有参数能量释放Gc为失效判据，不需要对裂纹的方向、大小进行显式描述，具备简便的三维几何与拓扑结构的处理方式，通过平均化能量克服了由于网格分布形式而带来的收敛性问题和裂纹路径网格相关的问题，严格的数学证明收敛于Griffith解[39]。此外，由于EigenFracture方法从材料失效物理机制出发，考虑了材料内部能量耗散的各种机制，包括塑性变形、裂纹扩展和相变，使得精确预测脆性或者韧性材料在不同温度以及载荷下的裂纹扩展成为可能。

3 应用场景与结果 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析

本节主要探讨OTM/HOTM极限力学仿真方法在小行星防御技术中的应用方向。

3.1 动能撞击

动能撞击防御技术是指撞击器以一定的速度和角度撞击小行星，使其自旋状态和轨道发生改变。撞击器可以选用航天器、火箭甚至可操控的小行星。这一技术的关键在于掌握动能撞击过程的动态响应和动能撞击的能量传递规律。动能撞击防御方法技术简单、启动迅速、灵活性强、作用效果明显，是一种实际可行的成熟技术。由于动能撞击过程中涉及高温、高压、高应变率等极限条件下材料的高度非线性动力学响应，在发生超高速动能撞击的局部区域，小行星表面材料在瞬间经历高温、高压、高应变率极端状态后，将发生大变形、断裂、破碎及熔化、气化乃至等离