1.研究背景

冲击动力学问题在工业装备和军事科技上有着广泛而重要的应用，如鸟撞击飞机、导弹穿透建筑物、震动致管道碰撞等。冲击动力学行为的研究对于重大装备和工程的结构设计和安全预测至关重要，该方向一直是学术界和工程界关注的热点。数值方法作为一种研究冲击动力学问题的有效方法，相较于实验手段，具有成本低和周期短等优势，逐渐受到研究人员青睐。有限元方法作为当前发展最为成熟的数值模拟方法，在冲击动力学领域也有广泛应用。但针对冲击/爆炸工况下的结构变形分析，有限元方法常因网格畸变而影响计算精度。为了克服这个缺点，有限元引入了网格重剖分技术，但这也造成了计算量的大幅增加。

鉴于有限元方法的这些缺陷，无网格粒子类方法应运而生，其中物质点法作为一种可便捷处理大变形和接触问题的无网格质点类方法，自提出以来，因其易于编程实现和精度高等优点在计算力学领域得到了快速发展。物质点法拥有两套描述方式，即基于欧拉描述的背景网格和基于拉格朗日描述的物质点，可以有效避免网格畸变和网格重剖分。该方法目前已成为处理高速冲击、爆炸和极端大变形等问题的有效数值模拟手段。

2.物质点法介绍

图1给出了物质点法拥有的两种离散形式，一是用于定义物体几何构型及携带物理量的物质点离散，二是用于求解离散控制方程的背景网格离散。该方法结合了欧拉法和拉格朗日法的优点，物理信息存储在物质点上，拉格朗日的物质点便于跟踪材料界面，引入与变形历史相关的本构模型。物质点上携带的物理信息通过插值函数映射到背景节点上，在背景网格上建立平衡方程进行计算，由于背景网格保持固定不变，因此可以有效避免网格畸变。

图1 物质点方法中背景网格和物质点示意图

物质点法求解格式与有限元法类似，也可分为隐式求解和显式求解两大类。其中显式方法是条件稳定算法，不需要求解线性代数方程组，单步计算量小，但其稳定时间积分步长也很小，常用于分析冲击/爆炸等关注应力波传播过程的高频响应问题。隐式算法一般是无条件稳定算法，需要求解大型线性代数方程组，单步计算量大，但其积分步长可比显式算法高2~4个数量级，适用于求解长时低频响应问题。Newmark法和中心差分法分别是物质点法中常用的隐式积分和显式积分算法。

物质点法大部分研究和应用均基于显式求解格式，其主要更新格式包括以下三类：USF(update stress first) 格式、USL(update stress last)和MUSL(modified update stress last)格式。研究表明，USL具有较强的能量数值耗散，而USF和MUSL格式具有较好的能量守恒性质，其中MUSL格式不直接利用更新的结点动量来计算结点速度, 而是将更新的质点动量映射回背景网格后重新计算结点速度，应力更新时先基于速度场计算应变率, 计算精度相对更高。

图2给出了标准显式物质点法的基本求解过程，其主要包括以下步骤：

图2 标准物质点方法的四步计算流程 (引自de Vaucorbeil, et al.AAM, 2020(53))

Step 1：物质点到背景网格的信息映射。定义背景网格, 将所有质点的质量、动量和内力项映射到背景网格结点上;

Step 2：背景网格的信息更新。更新背景网格节点的动量和位移信息;

Step 3：背景网格到物质点的信息映射。采用背景网格结点的速度场和加速度场更新质点的位置和速度，并利用节点速度场计算应变增量和旋量增量，根据本构方程更新应力;

Step 4：背景网格重设置。丢弃已经变形的背景网格，重新设定未变形背景网格。

上述步骤中，物质点与背景网格间的两次信息映射至关重要，其直接影响着物质点法的计算精度和计算效率。

3.物质点法插值技术

物质点法中一个重要的研究内容是关于映射插值函数的构造，其架起了物质点与背景网格节点间的信息传递，且对物质点法的计算精度至关重要。早期，Bardenhagen等（Bardenhagen, et al. CMES, 2004(5)）以变分形式和Petrov-Galerkin离散格式来推导物质点法，提出了广义插值物质点法（GIMP），用于解决基于delta函数插值方式的标准物质点法存在的严重跨网格误差，并指出这种跨网格误差来源于插值函数缺乏足够的光滑性。

随后受GIMP思想的启发，涌现出许多改进的算法，例如相邻粒子GIMP、双域物质点法（DDMP）、一阶和二阶对流粒子域插值法（CPDI和CPDI2）、B样条物质点法（BSMPM）和交错网格物质点法（SGMPM）等。这些方法都有一个共同的特征，即构造更为光滑的插值函数梯度。

这里，我们主要介绍CPDI型插值函数构造方式，其插值函数及其梯度的计算公式如下：

其中代表物质点p控制域的节点坐标,X代表物质点坐标,和分别代表背景网格和质点控制域的形函数,为质点控制域节点数，详见图3。

图3 对流粒子域型插值方式

4.数值算例

目前已初步实现基于标准插值和对流粒子域型插值的三维显式物质点法计算功能，对三维Taylor杆冲击和子弹高速冲击靶板算例进行测试，验证计算精度。

4.1算例1 Taylor杆冲击

该算例模拟Taylor杆冲击问题，图4给出了Taylor杆冲击前后示意图。几何模型尺寸根据Johnson和Holmquist实验，杆初始长度 =25.4mm，初始直径 =7.6mm，初始速度 =190m/s，本构采用Johnson-Cook本构模型。Taylor杆底部设置为z方向约束。整个模型离散为13,413个物质点，背景网格采用均匀六面体网格，尺寸为0.76×0.76×0.76mm³，分析总时间为8×10s。其中，插值方式采用CPDI2，质点控制域形状为八节点六面体。计算结果和实验结果、标准物质点法和有限元法解进行对比。

图4 Taylor杆冲击前后示意图

图5 最终状态变形位移图（左，x方向，右，y方向）

图6 最终状态等效塑性应变云图

表1 最终状态位移结果对比

4.2 算例2 弹丸高速冲击靶板

该算例模拟弹丸高速冲击靶板问题，图7给出了模型示意图。弹丸初始速度 =6,600m/s，本构采用带失效模式的Johnson-Cook本构模型。整个模型离散为184,058个物质点，背景网格采用均匀六面体网格，尺寸为0.5×0.5×0.5mm³，分析总时间为1×10s。其中，插值方式采用CPDI2和标准MPM插值，质点控制域形状为八节点六面体，此算例用到了质点控制域动态转化为标准MPM技术。计算结果和实验结果进行对比。

图7 弹丸高速冲击靶板示意图

图8 最终状态碎片云示意图（上 物质点法结果，下 实验结果）

图9 冲击过程动画

图10 冲击过程中CPDI2动态转换为标准MPM

5.结论

冲击、爆炸等问题是工程数值模拟中的重要分析内容，也是当前技术中的难点问题。显式物质点法作为一种有效的新型数值模拟技术，可以处理极端大变形、侵彻、破坏等位移强间断型物理问题。目前实现的三维显式物质点法计算功能，通过典型算例测试，计算精度和计算效率良好。后续将在此基础上，进一步开发新的材料本构模型、多物理场耦合等分析功能，并逐步完善插值技术和MPI并行功能。