OTM/HOTM极限力学仿真在小行星防御中的应用
时间：2023-03-09 21:35:02 来源：柠檬阅读网 本文已影响 人 
廖祜明，龚自正，宋光明，樊宗岳，杨宏涛，黎波*
(1.北京航空航天大学能源与动力工程学院，北京 102206；
2.北京卫星环境工程研究所，北京 100094；
3.北京大学工学院，北京 100871；
4.云翼超算(北京)软件科技有限公司，北京 100027)
地球自诞生以来就不断受到小行星(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方法，模拟了小行星灾难性破坏的过程，结果展示了材料屈服极限、剪切强度和孔隙率在碰撞中的相对重要性；
姜宇[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]。材料模型方面，与航天飞行器经过人类精心设计的材料和结构不同，小行星的材料和结构是自然形成的，不仅呈现出各向异性的特征，而且充满孔隙和裂纹，因而对典型材料模型和高状态相变数据都提出了更多的需求[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极限力学仿真理论框架进行介绍，并探讨其在陨石超高速撞击成坑、动能撞击、激光烧蚀驱动、核爆拦截等偏移小行星轨道以及摧毁小行星等防御措施中的潜在应用方向和相关案例，为小行星防御提供相关理论与技术支撑。
本节简要表述极限力学问题的理论基础与基于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方法从材料失效物理机制出发，考虑了材料内部能量耗散的各种机制，包括塑性变形、裂纹扩展和相变，使得精确预测脆性或者韧性材料在不同温度以及载荷下的裂纹扩展成为可能。
本节主要探讨OTM/HOTM极限力学仿真方法在小行星防御技术中的应用方向。
3.1 动能撞击
动能撞击防御技术是指撞击器以一定的速度和角度撞击小行星，使其自旋状态和轨道发生改变。撞击器可以选用航天器、火箭甚至可操控的小行星。这一技术的关键在于掌握动能撞击过程的动态响应和动能撞击的能量传递规律。动能撞击防御方法技术简单、启动迅速、灵活性强、作用效果明显，是一种实际可行的成熟技术。由于动能撞击过程中涉及高温、高压、高应变率等极限条件下材料的高度非线性动力学响应，在发生超高速动能撞击的局部区域，小行星表面材料在瞬间经历高温、高压、高应变率极端状态后，将发生大变形、断裂、破碎及熔化、气化乃至等离子体化等复杂的力学、物理过程及其耦合作用。温度从300K增到104K，压强从0.1MPa增至1TPa，应变率高达107s-1。其中，熔化和气化是超高速撞击下材料动态响应的重要特征，是超高速撞击过程中能量转化的主要表现形式[3]。对这些现象的描述直接影响着动能撞击过程中小行星动态响应和动能撞击的能量传递规律的准确性。基于OTM通过有机结合物质点空间离散、最优运输理论时间离散、局部最大熵无网格近似，以及基于物理的裂纹扩展算法和固液气全域材料模型，为动能撞击仿真提供了解决方案，如图5所示。
图5 基于OTM的超高速撞击仿真模拟
在OTM超高速撞击仿真中，采用J2粘塑性模型(J2-power law viscoplasticity)描述材料应变硬化、应变率硬化和高温软化等效应。高温高压下材料变形与温度、压力的关系可采用SESAME状态数据库或Grüneisen状态方程。通过采用Lindemann修正模型[40]描述熔化温度的体积依赖性。目前该方案已成功应用于金属[41-43]、尼龙[44]、陶瓷[45]、冰[46]等材料冲击碰撞仿真，如图 6所示。图 7为利用OTM方法初步模拟动能撞击过程中的弹性波传播、塑性变形、裂纹扩展，以及由于塑性与状态方程引起的温度升高而导致的材料相变之间的相互耦合与竞争关系的过程。
图6 OTM超高速撞击仿真应用
图7 小行星动能撞击仿真
为了详细描述小行星内部结构在超高速撞击效应下的动态热力学响应，可对复杂的小行星体建立分层复合材料模型，基于OTM的多场耦合变分理论框架，将对小行星体撞击过程中的弹性波传播、塑性变形、裂纹扩展以及由于塑性与状态方程引起的温度升高而导致的材料相变之间的相互耦合与竞争关系进行精准的描述，量化状态方程、材料强度以及热效应对超高速撞击下行星体动态热力学响应的影响。
3.2 激光烧蚀
激光烧蚀驱动技术是指使用强激光照射小行星表面，利用表面烧蚀产生的等离子体喷射所带来的反作用力，使小行星自旋状态和轨道发生改变。大量研究认为，激光烧蚀驱动技术是一种高效的空间碎片清除技术[47-50]，基于同样的作用机理可用于小行星防御。激光驱动改变小行星轨道是一种非接触式防御方法，其关键问题在于烧蚀驱动小行星机理、对驱动效果的影响因素、驱动的动力学模型。
在HOTM方法中，热流固耦合整体求解框架结合固液气相敏全域材料模型为小行星的激光烧蚀仿真提供了理想解决方案，而施加在小行星表面的激光可通过热通量边界条件进行模拟，顶部平坦的超高斯分布热流可以较好地近似物理激光束[51]，该激光模型可以被描述为：
(16)
式中：A为材料的激光吸收率，P为激光功率，r为激光半径，xc(t)、yc(t)为激光中心坐标，可用任意随时间的变化的函数描述。通过这种模型可容易地控制激光轨迹、速度等参数。图 8(a)表示n取不同值时的超高斯热流的分布，其中横坐标为该点到激光圆心距离与激光半径之比，纵坐标为取值。在HOTM方法中使用光线追踪算法来动态获取材料表面施加热流条件的节点集。计算域内其他部分的加热则通过热传导来实现。除了激光热通量之外，还可以施加热对流边界条件和热辐射。
在HOTM方法中，小行星表面物质气化形成的等离子体对小行星形成的反冲压力可由施加在表面结点的力边界条件实现，具体为采用Anisimov[52]等提出的随温度变化的模型来模拟反冲压力，表达式如下：
(17)
图8 (a)高斯阶数n取值不同时的高斯热流分布；
(b)金属铝的反冲压力曲线
图9 (a)金属粉末床激光选区熔融工艺(SLM)仿真；
(b)超高速激光金属熔覆工艺(EHLA)仿真
激光烧蚀技术使用大功率激光照射到小行星表面使其表面物质温度急剧上升，熔化及气化并产生等离子体。剧烈气化的物质向外喷射产生很大的反冲压力推动小行星运动，使其速度发生变化并改变轨迹，如图 10所示。相比于通过爆炸改变PHA轨迹的方法，激光实际部署地点与PHA一般相距较远，实际照射到PHA表面的激光能量密度较低，所以不存在PHA表面碎片飞出的风险，因而采用这种技术更加安全。
图10 小行星激光烧蚀驱动仿真
3.3 核爆拦截
核爆是应对短预警时间、大尺寸小行星撞击的技术。核爆炸防御方式有两种：一是利用核爆装置直接炸毁PHA；
二是利用核爆直接炸毁小行星或改变小行星的轨道以避免其与地球相撞，是近地小行星防御最主要的手段之一。根据PHA尺寸、材质、结构的不同，可选择表面爆炸、对峙爆炸以及穿透爆炸三种方式[4]。①表面爆炸：针对小体积的PHA，可以采用作用能量较大的表面爆炸或浅地下爆炸的方式，使PHA分裂成数块碎片；
②对峙爆炸：对峙爆炸是在距离PHA表面一定距离时引爆核装置，爆炸产生的热中子、X射线以及γ射线辐射PHA表面，产生高温，引发PHA表面物质的喷射，喷射时产生的推力使PHA发生偏转。此外，爆炸产生的部分碎片与PHA发生撞击，传递动能。两种作用效果叠加，实现防御目的。对峙爆炸是规避爆炸碎片威胁的有效方法之一，适用于防御体积较大的PHA；
③穿透爆炸：穿透核爆炸是指核装置穿入PHA内部一定深度处发生爆炸。该方法的优势在于，除了核爆产生的爆炸能量外，爆炸引起的表面冲击波能够扩大作用威力，穿透深度很浅的爆炸就足以改变PHA的运行轨道。
利用OTM方法可从核弹侵彻开始模拟侵彻过程因高速冲击压力、塑性变形升温、摩擦生热等共同作用引起炮弹内部炸药点火起爆、爆轰波在小行星内部传播、爆炸产生的高能气体膨胀流固耦合结构破坏的全过程。其中，采用Lee-Tarver三项点火增长模型[54]描述炸药点火起爆，在此模型中，点火、生长和完成分为三个过程。当炸药完成点火增长起爆后，爆炸产生的爆轰波气体以JWL状态方程[55]的形式膨胀对结构发生破坏作用，JWL状态方程由Jones、Wikins及Lee提出的一种不显含化学反应、由实验方法确定参数的半经验状态方程，能较精确地描述爆轰产物膨胀驱动做功过程。JWL状态方程及其等熵方程由三项组成，第一项主要在高压区起作用，第二项在中压区起作用，第三项在低压区起作用。爆炸产生的高温气体先冲击结构表面，然后冲击波在结构内部传播，可引发结构的大范围碎裂效应。通过结合EigenFracture裂纹扩展算法，可描述爆轰波在结构内部的流固耦合破坏作用。图 11初步展示了核弹穿入PHA内部发生爆炸引起的表面冲击波和核爆产生的爆炸能量共同作用使PHA分裂成许多碎片的过程。
图11 小行星核爆拦截仿真
针对小行星防御数值仿真在基础理论、计算方法、材料模型与高性能计算方面存在的困难与挑战，本文介绍了OTM/HOTM极限力学仿真理论。OTM/HOTM有机融合了多物理场变分原理、基于物理的本构理论和基于能量的断裂力学理论，为小行星防御中综合考虑热流固耦合、热力化学耦合等多物理场效应强耦合复杂的物理化学和力学现象提供了理论基础；
在数值离散方面，OTM采用了最优运输时间积分、物质点空间离散、局部最大熵无网格近似技术，并采用多线程多进行混合并行策略实现大规模并行计算，克服了基于网格的数值方法和无网格方法在精度、稳定性、健壮性、效率等方面的不足；
基于这些特性所发展的ESCAAS极限力学仿真平台实现了在统一的框架下求解材料超大变形、熔化、气化、冲击爆炸、流固与热流固耦合、自由表面、多体接触、碎裂、层裂与碎片云等复杂物理现象的多物理场、多尺度强耦合问题；
最后，展示了OTM/HOTM方法及ESCAAS平台在陨石超高速撞击成坑、动能撞击、激光烧蚀驱动、核爆拦截等小行星防御措施仿真中的有效性。由于OTM理论框架具有良好的可扩展性，后续可在小行星材料模型、动能撞击喷出的岩石羽流、破碎的石块和尘埃、激光烧蚀驱动小行星机理、对驱动效果的影响因素、驱动的动力学模型、高能量密度的核爆炸药材料模型、核爆机理及对小行星的破坏效果等一系列融合新的物理化学现象、新的材料模型和更加真实的边界条件方面进一步开展研究工作，形成小行星防御专用仿真平台，更加真实高效地仿真各种小行星防御场景。
