说说超高速碰撞仿真
本文摘要（由ai生成）：
这篇文章主要介绍了超高速碰撞的概念、需求、模拟难点以及仿真方法和软件，还介绍了 ESCAAS 软件的特点、算法和验证结果。文章指出，超高速碰撞是速度在 5km/s 以上的碰撞现象，典型应用包括航天器防护、反弹道导弹技术等，模拟难点在于物理现象的复杂性和算法、离散精度、材料模型及参数的选取对模拟结果的影响。目前常用的仿真方法包括基于网格的拉格朗日法、基于网格的欧拉法、ALE 任意拉格朗日-欧拉法和 SPH 无网格法等，商业非线性动力学仿真软件在超高速碰撞领域的应用较多，但也存在一些局限性。ESCAAS 软件是创新 OTM 方法的新一代软件，具有超高速碰撞仿真的特点，结合了最优输运理论、局部最大熵插值函数和物质点几何离散等核心算法，采用多进程多线程混合并行方案，实现了高效稳定、精确收敛的计算方法，还提供了用户自定义功能。通过与实验数据的对比，ESCAAS 计算结果的准确性得到了充分验证。
什么是超高速碰撞？
超高速碰撞一般是指速度在5km/s以上的碰撞现象，此时加载速度已接近或者超过了材料内部弹性波的传播速度。典型应用如卫星、飞船、航天飞机、空间站等航天器的防护，反弹道导弹技术、陨石撞击下的宇宙飞船防护技术，小行星撞击、类地行星表面陨石坑研究、核反应堆外安全壳安全防护设计、小行星采矿等等。航天领域是超高速碰撞应用的主要领域。
航天宇航活动为什么必须考虑超高速碰撞？
伴随人类越来越多的航天发射活动，地球轨道内（高度200~36000km）产生了大量的太空垃圾即空间碎片，另外还有微流星体等。空间碎片平均相对碰撞速度约为10km/s，而微流星体相对碰撞速度更高，可达20km/s。空间碎片或微流星体与航天器发生超高速碰撞时，巨大的碰撞动能足以破坏航天器结构，使其功能失效，甚至导致灾难性的后果。大的空间碎片碰撞可以导致航天器彻底毁坏，且几乎无法防护，唯一有效的方法是在轨机动规避；小空间碎片可以通过设置防护结构来加以防护。超高速碰撞试验代价十分昂贵，需要二级或三级气炮、甚至电磁炮、定向聚能加速器等，动态测量，高速照相（>1,000,000帧每秒）等技术，试验成本高、周期长，进行大量试验相对困难。因此航天领域对能够精细进行超高速碰撞现象仿真的软件技术需求十分迫切。
超高速碰撞模拟难在何处？
超高速碰撞模拟与低速碰撞、高速碰撞的区别在于物理现象不同。低速碰撞时，属于结构动力学问题，局部侵彻与结构的总体变形效应紧密耦合；高速碰撞，碰撞点附近区域靶材料性质（密度、强度和断裂韧性）起主要作用，结构效应退居次要地位；超高速碰撞，材料的惯性甚至可压缩性效应或相变效应会起重要作用；极高的超高速碰撞（大于12km/s）时，撞击区能量沉积速度很快，以致发生气化爆炸现象。超高速撞击后激波在材料中传播会同时激活各种相互耦合与竞争的能量耗散机制，热力强耦合，引起材料的非弹性大变形、相变（液化、气化、凝固等）、裂纹萌生及扩展和层裂、碎裂等，更高的能级也会造成等离子体的产生。
算法、离散精度、材料模型及参数的选取会在很大程度影响数值模拟结果的准确性和可信性。如何准确将各种能量耗散机制及其相互耦合作用准确考虑进去，对仿真精度具有重要的影响。
空间碎片对航天器的破坏效果是随碎片质量、碎片材料、碎片形态、撞击角度和撞击速度而变化的，航天器的防护层设计及其防护能力需要考虑大量不同碰撞场景，进行大量的仿真计算，从而优化防护设计。超高速碰撞过程中，弹丸会形成固、液、气三种状态混合的碎片云。同时将热流固强耦合、塑性变形、脆性与韧性裂纹扩展以及破碎、熔化、气化等现象很好的模拟出来，对商业软件是一个考验。
超高速碰撞仿真方法及软件
超高速碰撞仿真本质上属于高度非线性动力学问题，常用方法包括基于网格的拉格朗日法（比如有限元方法FEM）、基于网格的欧拉法（例如有限体法FVM与有限差分法FDM）、任意拉格朗日-欧拉法（ALE）以及无网格法（meshfree）等，其中无网格方法又有光滑粒子流体动力学法（SPH）、物质点方法（MPM）、最优输运无网格方法（OTM）等。
图1 采用基于网格的拉格朗日法、基于网格的欧拉法、ALE任意拉格朗日-欧拉法和SPH无网格法模拟碰撞问题的例子[1]。
基于网格的拉格朗日法在整个计算过程中网格附着于材料上，网格随着材料的运动而运动，对于超高速碰撞这类典型的超大变形问题，由于巨大的网格畸变、扭曲、缠绕或单元分裂等问题，将造成有限元求解困难甚至导致求解失败。为了解决这些问题，往往在有限元计算中需要频繁的进行网格重分、在新旧网格之间交换数据，使得计算过程繁琐复杂，大大地增加了计算时间，同时导致可观的数值误差。或者删除畸变网格，引入数值裂纹，破坏各种守恒关系。对于有些问题仅仅重新划分网格并不能完全得以解决，最后往往导致计算时间步长过小或者数值误差过大从而程序提前终止，因此这种方法不适用于超高速碰撞的仿真。
基于网格的欧拉法，网格固定在空间中，材料在网格中流动，适于处理材料大变形问题，处理流体问题较为理想。欧拉法不存在网格问题，但是其空间描述特性在追踪材料界面等具有固有的缺陷，需要对较大的空间进行离散，计算量大，需要复杂的算法与极大的计算成本来确定不同材料和材料不同相之间的交界面，此外无法精确模拟历史与应变率相关的材料热力学行为（比如塑性变形），无法精确模拟裂纹扩展、多体碰撞与物质的自由表面。因此这种方法模拟超高速碰撞也有难以克服的局限性。
任意拉格朗日-欧拉法（ALE），通过引入一个参考域并采用网格离散，参考网格能够独立于拉格朗日坐标系和欧拉坐标系以任意的形式运动，可以用来处理界面协调及自由表面的问题，并且适用于有限元方法以及边界元方法，ALE结合了基于网格的拉格朗日与欧拉法的优势，但同时也继承了拉格朗日的大部分缺点，即便是当前最先进的ALE方法，在面对大变形问题时也达到了该方法本身的极限，在大变形和多相无约束流问题中因网格运动而导致的单元拉伸和扭曲，ALE必须采用重分网和状态变量的重映射，每个计算步骤需对新网格进行规整化，意味着将带来较大的数值扩散误差。而且在大变形问题中，网格重新划分也会导致物理参数在新旧网格之间的传递而引入插值误差，此外由于ALE方法的网格定义于不同的坐标系中大规模并行化也是ALE方法面临的一大困难。
光滑粒子流体动力学法（SPH），是近30多年发展起来的一种纯拉格朗日方法，与有限元或有限差分等基于网格划分技术的传统数值算法相比，SPH方法通过在计算域中填充具有独立材料性质的SPH粒子来替代网格划分，粒子遵守质量、动量和能量守恒方程。SPH方法的无网格粒子（节点）不只是作为插值近似点，同时还携带有材料性质，在内、外力作用下运动与拉格朗日算法相结合。SPH方法守恒方程的构建不受SPH粒子分布随意性的影响，可以很自然处理材料发生极大变形（如破坏、碎裂）下的非线性动力学问题。为此SPH在超高速碰撞领域有了非常广泛的应用。但随着应用的深入，其缺点也逐渐显现。主要问题在于插值函数不满足Kronecker-delta属性，无法直接施加位移边界条件，插值函数不满足严格的非负性，求解热力学问题不收敛，缺乏有效的数值积分方式，缺乏基于物理的裂纹扩展算法，严重的零能模式及拉应力不稳定性等。物质点方法（MPM），结合了基于网格的拉格朗日法与欧拉法的优点，具有一定的进步，但在拉格朗日物质点与欧拉背景网格间的数据交换引入数值误差，缺乏基于物理的裂纹扩展算法，计算结果与质量点的分布有关，计算成本较高等。
最优输运无网格方法（OTM），是近十年发展起来的创新方法，集成了拉格朗日方法和无网格方法的所有优点，具有严格的数学证明，包括时间积分的动量守恒与辛守恒，质量守恒自动满足，具有严格的误差分析；采用高阶插值函数，满足Kronecker-delta特性，可直接施加位移边界条件。缺点是计算成本较高。ESCAAS软件作为第一款采用OTM方法的商业软件，其创新的多层混合并行化方法pOTM，可线性甚至超线性加速到个计算核心，具有卓越的处理超大规模问题的计算能力，这对如今硬件成本大大降低，超算云计算遍地开花的时代来说，这种计算方法的应用可以说潜力无限。
图 2：(a) OTM物质点与节点相结合的方式对几何模型进行空间离散的示意图[2]，系统的运动学信息存储在节点，而力学信息存储在物质点。物质点离散为无网格法中的数值积分提供了快捷的解决方案。同时物质点与节点动态联系，无需像自适应有限元方法或者动态网格细化一样需要在新旧网格中进行信息的转换，避免了传统自适应求解器中引入的极大的数值误差。由于插值与积分在不同的离散点进行，完全避免了SPH方法中存在的拉应力不稳定性（Tensile Instability）[3,4]。在每个物质点采用高阶高斯积分可有效的规避零能模式的出现。; (b) pOTM多进程多线程动态负载平衡混合并行方案[5]。
如今商业非线性动力学仿真软件有LS-DYNA，AUTODYN,ANSYS-DYNA,DYTRAN,PAM-CRASH, ABAQUS,RADIOSS等等。这些商业软件基本都是采用传统的模拟方法，如显式有限元，欧拉法，ALE，SPH等。这些商业软件，尤其是LS-DYNA，AUTODYN在空间碎片领域的仿真应用较多，也取得了一定的成果。但随着对仿真精细化的要求越来越高，对物理机理的研究要更多细节更精确的展现，真正做到对试验的预测而不是拟合，这些商业软件的局限就非常明显了。例如算法本身有很多固有的缺陷，很多现象尤其是能量的各种耗散方式的强耦合与竞争关系，裂纹扩展模式等不能得到很好的解决，另外很难进行二次开发，有些软件根本不支持二次开发。ESCAAS作为创新OTM方法的新一代软件，还需要用户更多的使用和验证，应该说为广大的相关领域工作者提供了一个创新的工具，可以用它做很多之前想做而做不了的工作，对ESCAAS和广大科技工作者都是难得的机遇。
ESCAAS软件适于超高速碰撞仿真的特点
为了高效稳定的求解极大变形、高速冲击及几何畸变、金属材料成型及多相耦合等问题。OTM方法采用物质点与节点对初始问题域进行离散，采用局部最大熵插值函数来构造连续的运动函数，避免了有限元方法处理极大变形网格畸变，无网格法无法直接添加Dirichlet边界条件，以及计算不收敛等问题。同时由于插值与积分在不同的离散点进行，提供了有效的无网格数值积分模式，解决了拉应力不稳定性的问题。OTM方法作为增量更新拉格朗日方法，质量守恒自动满足，无需求解。通过采用最优运输理论进行时间离散，保证了离散系统的动量守恒以及辛守恒，极大的提高了运算的效率与精度。
经过十几年的研发以及理论与大量实验验证，OTM方法为当今世界上唯一能够精确稳定的模拟不同材料在高温（>103~104K），高压（>10~100Gpa）、高应变率（>104~108/s）等极限条件下复杂热力学响应的计算方法。能够高精度模拟从一般碰撞到超高速碰撞问题（0~10km）。
能够很好考虑起决定性作用的复杂的热力学现象，其中包含了复杂三维结构的超大变形、热力强耦合过程以及固液气相变与多相混合、高度非线性以及应变率相关的材料热力学响应、三维裂纹扩展与碎片云，以及多体动态接触。
采用了一种全新的裂纹扩展算法，物质点失效算法（EigenErosion）。EigenErosion方法是基于能量与变分原理的裂纹尖端扩展算法。数学上可以严格的证明EigenErosion的计算结果是收敛的，即离散越细，结果越精确。是无网格方法中难得的基于物理的，收敛的，与离散无关的裂纹扩展算法。
OTM方法中，节点包含了运动学信息与温度场。通过热力耦合变分原理，对能量守恒与热力学第二定律的求解，OTM方法可以精确的预测材料内部的温度分布以及系统熵的变化。尤其是金属材料在高速撞击与高应变率加载条件下，能量会通过粘性以及塑性等内部耗散过程转换为热能。ESCAAS是现今唯一可以精确模拟带能量耗散的热力强耦合系统的软件。拉格朗日无网格法与热力耦合变分原理的完美结合，为热流固耦合，固液气动态相变与多相混合提供了高效稳定、精确收敛的计算方法。
ESCAAS采用多进程多线程（MPI/PThread）混合的方式对OTM方法进行了高度并行化（pOTM方法）来降低计算成本。pOTM方法已成功地将OTM计算的效率超线性的提升到个计算核心。
ESCAAS采用基于设计模式程序架构和高效的C++实现有机结合，保证了商业化产品级的健壮性与扩展性，用户可以充分利用ESCAAS的开放性架构，将自己的创新成果（算法、材料模型、公式等）加入软件中，形成自己的专用工具。
ESCAAS碰撞仿真的例子
图 3 ESCAAS中先进的核心算法：OTM最优运输无网格法方法通过有机结合最优输运理论、局部最大熵无网格插值函数与物质点几何离散，并采用多进程与多线程双层混合大规模并行方案。
图 4展示了ESCAAS对碰撞冲击试验的数值模拟，在弹丸以不同的速度撞击下，金属靶板的破坏与失效呈现出了截然不同的物理机制：在低速弹道冲击试验中，靶板以常见的冲塞机制被穿孔；在高速冲击下，靶板产生了严重的塑性变形与局部的熔化现象。然而，超高速撞击过程中，急剧增加的压力使材料内部的温度远远高于其熔点，固液气相变随之产生。但是，随着板厚的增加，超高速撞击实验中能量的耗散机制同时也伴随着严重的塑性变形与动态裂纹扩展（比如层裂）。ESCAAS计算结果的精确性通过与几百次低、中、高以及超高速碰撞冲击试验（撞击速度从100m/s到10km/s）进行比较得到了充分的验证。
图4 冲击速度从100m/s到10km/s的金属靶板冲击试验的OTM仿真
图 5 ESCAAS模拟航天器防护层超高速碰撞：Al6061-T6铝盘被尼龙圆柱体高速撞击。铝盘厚度范围 [0.5,3.0]mm, 入射速度范围[5, 7]km/s，入射角范围[0, 60]。铝盘在不同的加载条件下，展现了完全不同的失效机制，其中包括粘弹塑性，材料的熔化与气化以及液固气多相混合，裂纹扩展与碎片云等。
图 6 ESCAAS模拟航天器防护层超高速碰撞，基于模型的不可靠性算法（Model-based UQ）通过结合ESCAAS模拟与实验，严格的量化不同板厚的铝合金标靶在不同速度、不同角度的超高速打击下的失效概率：每个点对应了4个ESCAAS模拟分别采用了不同的盘旋角度（假设在碰撞试验中非圆形弹丸的盘旋是不可控变量）；每个OTM计算分别采用150万个物质点模拟标靶，20万物质点模拟弹丸每个ESCAAS模拟采用512 mpi 进程，可在10小时内完成计算（20微秒）。2500 次OTM数值模拟与45次实验数据对比: Error=11%。
参考文献
[1]Quan, X.; Birnbaum, N.K.; Cowler, M.S.; Gerber,B.I.; Clegg, R.A.; Hayhurst, C.J. Numerical Simulation of StructuralDeformation under Shock and Impact Loads Using a Coupled Multi-Solver Approach.In Proceedings of the 5th Asia-Pacific Conference on Shock and Impact Loads onStructures, Hunan, China, November 12–14, 2003.
[2]B. Li, F. Habbal, and M. Ortiz, “Optimaltransportation meshfree approximation schemes for fluid and plastic flows,”Int. J. Numer. Methods Eng., vol. 83, no. 12, pp. 1541–1579, 2010.
[3]Li, Shaofan, and Wing Kam Liu. "Meshfree andparticle methods and their applications." Applied Mechanics Reviews 55.1(2002): 1-34.
[4]Fernández-Méndez, Sonia, and Antonio Huerta."Imposing essential boundary conditions in mesh-free methods."Computer methods in applied mechanics and engineering 193.12 (2004): 1257-1275.
[5]B. Li, M. Stalzer, and M. Ortiz, “A massivelyparallel implementation of the Optimal Transportation Meshfree method forexplicit solid dynamics,” Int. J. Numer. Methods Eng., vol. 100, no. 1, pp.40–61, 2014.
[6]贾光辉，《航天器结构超高速碰撞数值模拟》
[7]林俊德等，超高速撞击实验的三级压缩气炮技术，爆炸与冲击，2012年9月
[8]邢柏阳等，强爆轰驱动超高速碎片发射装置设计因素分析
[9]ESCAAS软件介绍手册及相关论文
