超高速撞击下碎片云的OTM 分析*

廖祜明，黎 波，樊 江，焦立新，于帅超，林健宇，裴晓阳

（1. 北京航空航天大学能源与动力工程学院，北京 102206；2. 云翼超算(北京)软件科技有限公司，北京 100027；3. 北京大学工学院，北京 100871；4. 中国工程物理研究院流体物理研究所，四川 绵阳 621900）

地球轨道上散布着大量的空间碎片和微流星体，若它们与在轨运行航天器相撞，相对撞击速度可达十几千米每秒，对航天器造成巨大威胁。大型碎片撞击可直接撞毁航天器；中型碎片撞击会造成航天器部件功能失效；而小型碎片撞击可在航天器表面造成累积损伤，导致部件功能下降[1]。因此，空间碎片超高速撞击防护是航天器结构设计须重点考虑的问题。厘清空间碎片超高速撞击的力学行为是防撞设计的重要基础。碎片超高速撞击的物理过程非常复杂，撞击过程中激波在材料中传播，同时激活各种相互耦合与竞争的能量耗散机制，引起材料的非弹性大变形、相变（液化、气化、凝固等）、裂纹萌生及扩展和层裂、碎裂等复杂的物理过程。采用数值方法模拟超高速撞击面临如何解决速度过大导致的网格畸变，以及如何准确描述超高速撞击过程中各种能量耗散模式之间的竞争与耦合关系等巨大挑战。

目前分析超高速撞击问题的数值计算方法主要包括基于网格的方法和无网格方法两大类。基于网格的数值方法中，传统拉氏有限元法受网格畸变困扰，难以处理超大变形问题；欧氏方法具有网格不变形的固有优势，但如何跟踪历史变量和物质界面仍然是亟待解决的问题，同时在模拟动态裂纹扩展、多体碰撞与应变率相关的材料热力学行为时还需进一步研究；任意拉格朗日-欧拉（arbitrary Lagrangian-Eulerian, ALE）[2-3]法在一定程度上缓解了网格畸变带来的困难，但ALE 继承了拉氏有限元的大部分缺点，同样面临大变形时网格畸变的问题[4]。无网格法由于不需要进行网格离散及采用高阶插值函数，有利于解决大变形问题，在超高速撞击问题上应用更为广泛。常用的无网格方法有光滑粒子流体动力学法（smoothed particle hydrodynamics, SPH）[5]、再生核粒子法（reproducing kernel particle method, RKPM）[6]以及物质点法（material point method, MPM）[7]等。有大量的研究应用以上方法[8-9]及其改进方法[10]对超高速撞击问题进行仿真分析并取得了丰富的研究成果，但在算法的精度、稳定性、健壮性、效率等方面仍然存在许多不足[11]。

最优运输无网格方法（optimal transportation meshfree, OTM）[12]是一种专门针对动态冲击问题提出的显式增量更新拉格朗日无网格方法。OTM 方法采用局部最大熵无网格插值函数（local maximum entropy,LME）[13]，在边界上满足Kronecker-Delta 属性，解决了传统无网格法在处理位移边界条件时的缺陷；同时在节点的基础上增加物质点，以物质点作为积分点有效避免了计算结果在拉伸载荷下的不稳定性；在时间离散上采用了最优运输理论（optimal transport theory）[14]，保证了其时间离散形式在理论上严格满足动量守恒，具有严格的收敛性数学证明；在并行计算上结合了共享内存多线程并行与分布式多进程并行实现了大规模多层混合并行[15]，显著提高了计算效率；在材料失效方面，采用本征侵蚀裂纹扩展算法（Eigenerosion），以材料固有属性能量释放率作为失效判据[16-18]，克服了传统裂纹扩展数值方法无法清晰表征材料真实变形和失效物理机制的缺点，理论证明收敛到Griffith 理论，与网格无关。OTM 方法的这些特性为超高速撞击模拟预测提供了高效精准的解决方案[19]。本文采用基于OTM 方法的极限力学仿真软件ESCAAS，对空间碎片超高速撞击问题进行数值模拟，并与实验结果进行对比验证。

1 OTM 最优运输无网格方法

1.1 时间与空间离散

假设参考构型中的连续介质 Ω0R3运动过程由变形映射（位移） φ:Ω0×[t0;t]!R3描述，其中 [t0;t] 为运动过程的总时间。X2 Ω0表示材料质点在参考构形中的位置，x=φ(X;t) 表示材料质点在现时构形Ωt=φ(Ω0;t)的位置，则拉格朗日描述下质量守恒、动量守恒与能量守恒方程为：

式中：下标k代表时间步长索引。式(5)为时间离散的变分公式，为了便于计算机求解，要对此半离散的变分公式进行空间离散得到时间-空间全离散格式。

在有限元法中，每个单元有特定的形状且各单元间具有固定的连接关系，每个单元携带了材料一定量的质量并具有相应的体积，材料的响应在单元的积分点上进行求解。与有限元法类似，OTM 方法中采用物质点xp与节点xa相结合的方式对半离散增量变分公式进行无网格空间离散，其中下标a表示节点索引，p表示物质点索引。如图1 所示，所有的场变量信息分别由节点xa和物质点xp携带。物质点可以理解为有限元中的积分点，而其携带的体积可以看作是有限元中的单元，而不同于有限元的是节点和物质点不再有固定的连接关系，即此“单元”没有固定的形状，其形状将根据物质的变形可以是任意的形状，由于解放了网格的束缚，因此对于大变形问题不存在网格畸变带来的问题[21]。

材料的动力学信息存储在节点上，包括位移、速度、加速度与温度等。在计算过程中tk+1时刻节点位置xa;k+1和节点速度 φ˙a;k+1的更新公式为：

从tk时刻到tk+1时刻的位移传输（变形）映射 φk!k+1为：

式中：Na为LME 局部最大熵插值函数，NH(xp;k代表物质点xp在tk时刻下的邻域。

材料的物理信息存储在物质点上，包括质量、体积、密度、变形、应力与材料内部参数等。材料的响应在物质点上求解，物质点携带的质量在计算过程中保持不变，而体积和密度在不断的变化，每个时间步都得到更新。物质点的运动通过节点的动力学信息插值获得[22]。在计算过程中，物质点在tk+1的位置更新，通过对邻域内的节点进行插值获得：

物质点邻域范围NHxp;k在每个时间步进行更新，同时更新各节点上插值函数值。LME 插值函数Na满足严格的非负性，以及零阶和一阶连续性要求，同时在边界满足Kronecker-Delta 属性。此外，在LME 插值函数中，通过调整 γ 值，插值函数作用范围可从局部有限元无缝过渡到全局无网格，如图2 所示，这种衰减特性建立了与高斯径向基函数的联系，其重要作用是只有少量的节点对待求函数有明显的贡献，大大降低计算成本。

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

1.2 捕捉式多体动态接触算法

OTM 方法采用物质点和节点进行离散，并且动态更新局部邻域，使其具备自动捕获接触（seizing contact）的能力。具体地，物体 ΩA上的节点同时也可以作为物体 ΩB中物质点的邻域，因此对两物体的接触可以不加特殊处理，如图3 所示，其中Ra为节点xa的插值函数的作用范围，Rb同理，vp;k代表物质点速度，vb;k代表节点速度。

在捕捉式多体动态接触算法中，当碰撞节点进入物质点邻域后，其线性动量会被自然地抵消，如碰撞节点xa的线性动量la;k为：

式中：M为节点邻域内的物质点数量，为mp为物质点质量。节点的动量为其邻域物质点动量的平均值，结合运算中采用的集中质量法，节点的速度本身代表了局部物质点速度的质量加权平均值，平均了速度的所有分量，接触时物体之间的法向速度和切向速度都被用于相互抵消，从而自动确定了碰撞体之间的接触条件。因此，这种捕捉式接触时的摩擦系数可视为无限大，适用于模拟碰撞终端抛射物与目标物混合及伴随的绝热加热过程。

1.3 本征侵蚀裂纹扩展算法

材料在极限条件下的裂纹扩展过程包含了多种能量传播与耗散方式的结合与竞争，准确地描述这些特征需要采用基于物理的裂纹扩展算法。OTM 方法中采用能量的方式来描述整个材料的响应过程，提出了本征侵蚀（Eigenerosion）裂纹扩展算法。在该方法中，每个物质点代表一小块物质，通过移除失效物质点的方式来模拟材料中裂纹扩展和碎裂的形成，物质点失效代表该小块物质内部产生裂纹并形成自由表面，而自由表面的产生则必然伴随着能量耗散，裂纹路径是由弹性能量（弹性能量作为能量释放机制促进断裂）、材料特定的断裂能量（与裂纹面积成正比的断裂比能）、单调性和接触约束之间的竞争引起的结果。物质点失效的判断准则为：

式中：Gp表示物质点xp的等效能量释放率，下标p代表计算域内物质点的索引；Cη表示局部邻域， η 决定了其范围，kCηk 表示局部邻域总体积；λ为几何校正因子；Vq为局部邻域内物质点xq体积，下标q代表邻域内物质点的索引。如图4 所示，黑色点为一组失效物质点组成的裂纹，圆圈内的红点为位于裂纹尖端的物质点的邻域 η ，通过平均该物质点邻域内的弹性能We得到物质点的等效能量释放率Gp，若Gp大于材料本身的临界能量释放率Gc时 ，该物质点将失效不再参与计算。

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

2 铜飞片超高速撞击铝合金靶板数的值模拟与验证

2.1 数值模拟模型

文献[25]开展了柱形OFHC 铜飞片超高速撞击Al6061-T6 铝合金靶板实验，并运用高速摄像系统拍摄了碎片云形貌。本节将采用基于OTM 方法开发的ESCAAS 软件，建立同样的超高速撞击模型（如图5 所示），对该超高速撞击过程进行模拟并开展与实验结果的对比验证。

工况1：撞击角度为5.4°，质量为3 g 的柱形OFHC 铜飞片，以5.55 km/s 的速度撞击Al6061-T6 铝合金靶板，几何模型如图5(a)所示，靶板外侧施加固定边界条件，取铜飞片和靶板的1/2 对称模型进行无网格离散，物质点最小特征尺寸为0.14 mm，共计31.3 万个物质点，6.26 万个节点；

工况2：撞击角度为11.7°，质量为10 g 的柱形OFHC 铜飞片，以5.12 km/s 的速度撞击Al6061-T6 铝合金靶板，几何模型如图5(b)所示，靶板外侧施加固定边界条件；取铜飞片和靶板的1/2 对称模型进行无网格离散，物质点最小特征尺寸为0.6 mm，共计40.6 万个物质点，22.9 万个节点，如图5(c)所示。

2.2 材料模型

超高速撞击过程涉及材料应变硬化、应变率硬化和高温软化等效应，J2 黏塑性模型（J2-power law viscoplasticity）能较好地模拟铜和铝合金材料塑性响应，其等效屈服应力为：

表1 材料物性参数Table 1 Material parameters

表2 材料本构模型参数Table 2 Parameters of material constitutive model

2.3 结果分析

工况1：在Linux 平台下采用ESCAAS 模拟 3 g 柱形铜飞片，以5.4°的撞击角度、5.55 km/s 的速度撞击铝合金靶板的过程。整个模拟采用16 核并行计算，10 µs 撞击过程模拟用时约为1.33 h。图6 为铜飞片穿透缓冲墙形成碎片的过程，模拟结果很好地呈现出碎片云的形态，包括碎片云呈现出头部椭圆形、内核碎片云的位置、外泡碎片云的形态以及反溅碎片云的形态等特征信息。

除了碎片云轮廓形貌的定性对比外，碎片云特征尺寸定量对比方面：铝合金靶板前端碎片云轮廓轴向最大位置点Point 1 距靶板的轴向距离的OTM 模拟结果为(50.74±0.66) mm，实验测量结果为51.3 mm，与实验的相对误差为(1.09±1.29)%；铜飞片前端碎片云轮廓轴向最大位置点Point 2 距靶板的轴向距离的OTM 模拟结果为(31.9±1.4) mm，实验测量结果为33.5 mm，与实验的相对误差为(4.8±4.18)%；铝合金靶板反溅碎片云的高度测量点Point 6 距靶板的轴向距离的OTM 模拟结果为(7.87±0.74) mm，实验测量结果为7.76 mm，与实验对比相对误差为(1.42±9.54)%；铜飞片和靶板材料碎片径向相交距离（测量点Point 4 和Point 5 的径向位置）的OTM 模拟结果为(27.6±1.3) mm，实验测量结果为27.9 mm，与实验对比相对误差为(1.07±4.66)%。

在超高速撞击过程中，弹丸厚度/板厚(D/H)、打击速度与角度的结合决定了冲击波的分段分区效应。本例中D/H=1.2，属于中厚板情形，其主要失效模式为塑性变形、裂纹扩展、熔化三种能量转换与耗散方式相结合。如图8 所示，OTM 方法模拟了飞片撞击下，靶板发生塑性变形引起靶板翻边撕裂的过程。测量得到OTM 模拟的穿孔直径为(20.9±0.4) mm，实验测试值约为20.1 mm，相对误差为(3.98±1.99)%。从上述分析可见仿真量化结果与实验吻合良好。

工况2：模拟了10 g 柱形OFHC 铜飞片，以11.7°的撞击角度，5.12 km/s 的速度撞击Al6061-T6 铝合金靶板的过程。整个模拟采用16 核并行计算，10 µs 撞击过程模拟时间约为2.4 h。图9 所示为不同时刻下碎片云的轮廓，由于铜飞片撞击角度11.7°较大，铜飞片碎片云与靶板碎片云前端呈现出显著偏转。

图10 为7.6 µs 时刻下，OTM 模拟与实验测量的碎片云形貌对比情况。同样地，轴向位置点Point 1 为靶板碎片云前端；轴向位置点Point 2 为铜飞片碎片云前端；轴向点Point 3 为撞击之前靶板的前端；径向位置点Point 4 和点Point 5 为铜飞片和靶板碎片径向相交距离测量点；径向位置点Point 6 为靶板后端反溅碎片云高度。此外，由于撞击角度11.7°较大，铜飞片碎片云与靶板碎片云前端呈现出显著偏转，增加了碎片云偏转角度与实验测量结果的对比情况，包括碎片云前端偏离中线的角度 α ，靶板上侧反溅碎片云偏离角 θ ，靶板下侧反溅碎片云偏离角 β 。

铜飞片速度为5.12 km/s 时，铝合金靶板前端碎片云轮廓轴向最大位置点Point 1 距靶板轴向距离的OTM 模拟结果为(47.2±1.2) mm，实验测量结果为47.5 mm，与实验的相对误差为(0.63±2.53)%；铜飞片前端碎片云轮廓轴向最大位置点Point 2 距靶板轴向距离的OTM 模拟结果为(27.1±2.1) mm，实验测量结果为26.0 mm，与实验的相对误差为(4.23±8.08)%；铝合金靶板反溅碎片云的高度测量点Point 6 距靶板轴向距离的OTM 模拟结果为(32.0±2.6) mm，实验测量结果为30.9 mm，与实验对比相对误差为(3.56±8.41)%；碎片云前端偏离中线的角度 α的OTM 模拟结果为(10.55±0.65)°，实验测量结果为11.5°，与实验对比相对误差约为(8.26±5.65)%；靶板上侧反溅碎片云偏离角 θ 的OTM 模拟结果为(109.5±3.5)°，实验测量结果为100.8°，与实验对比相对误差约为(8.63±3.47)%；靶板下侧反溅碎片云偏离角 β 的OTM 模拟结果为(109.75±3.25)°，实验测量结果为108.0°，与实验对比相对误差约为(1.62±3.00)%，模拟结果与实验结果吻合理想。

图11 所示为工况2 铜飞片撞击铝合金靶板，本例中D/H=1.17 属于中厚板，同样地，在撞击过程中可以观察到铝合金靶板具有显著的塑性变形和翻边效果。OTM 模拟的靶板穿孔直径为32.125 mm，实验测试值约为31.6 mm，相对误差大约1.66%，穿孔孔径与实验高度吻合。

3 结 论

超高速撞击涉及材料局部超大变形、裂纹扩展、层裂与碎裂、多体动态接触、固液气相变与多相混合等复杂的材料动态响应过程，对传统基于网格的数值方法提出了巨大挑战。本文引入最优运输无网格方法（OTM）对空间碎片超高速撞击过程进行数值分析，指出了OTM 在处理此类问题中有如下优势：

(1) OTM 方法基于变分原理，在统一的拉格朗日框架下求解物质运动方程，实现了对不同材料在高温、高压、高应变率等极限条件下多种能量耗散模式结合与竞争的现象进行了完善的考虑，在理论上保证了不同形式能量耗散的自主耦合分配；

(2) OTM 方法采用物质点和节点空间离散，在时间与空间上的离散被严格地限制在变分框架内，保证了质量，能量及动量的守恒；结合局部最大熵无网格插值函数，克服了一般无网格法在精准施加位移边界条件上的不足，稳定性和收敛性有完备的数学证明；

(3) OTM 方法采用了真实的裂纹扩展算法，考虑了材料内部能量耗散的各种机制，包括塑性变形、裂纹扩展和相变，采用材料的固有属性能量释放率作为失效判据，而非简单采用单元变形量、应变大小，或是基于实验的失效模型来作为材料的失效判据，克服了传统裂纹扩展数值方法无法清晰表征材料真实变形和失效物理机制，以及不收敛和网格相关的缺点；

(4) OTM 方法及ESCAAS 软件在计算弹孔直径、碎片云形态方面与实验结果吻合理想，验证了其在超高速撞击数值模拟方面的适用性，显示出OTM 方法及ESCAAS 软件可以作为超高速撞击的有力数值分析手段。