改进物质点法及其在爆炸与冲击动力学中的应用

陈卫东， 刘澜， 路胜卓， 吴世博， 马敬鑫， 吴培文

(哈尔滨工程大学 航天与建筑工程学院，黑龙江 哈尔滨 150001)

物质点法(material point method，MPM)是一种由粒子算法与网格算法相结合的数值计算方法。该算法最早是由Sulsky等[1]率先提出的，此后国内外许多学者都对该算法进行了深入研究，并且在时间积分[2]、形函数改进[3-5]、接触算法[6]以及多种算法耦合[7-9]等多个方面对该算法进行了改进，从而使该算法不断趋于成熟。因此在近年来，该算法已经在许多工程计算领域中得到了广泛应用。

但是在爆炸与冲击类工程计算领域中，由于计算对象往往会出现体积剧烈膨胀或结构大尺度拉伸变形等状况，例如炸药爆轰产物膨胀扩散以及破片飞溅等。所以在应用物质点法及其现有的改进算法求解此类问题时，当模型运行至一定阶段后，往往会出现背景网格区域容量不足以及物质点陷入拉伸不稳定状态等不利情况，这些情况会对计算结果产生严重影响，从而制约物质点法在爆炸与冲击类工程计算领域中的应用。因此现有必要对物质点法作进一步改进，以使该算法能够更加适用于求解爆炸与冲击类工程实际问题。

1 物质点法的显式时间积分原理

在物质点法中，计算对象被离散成一系列占据一定体积的物质点，并且被设置于背景网格之中，物质点储存着其所在区域的材料的全部物理信息，而背景网格节点则只用于求解动量方程和空间导数。因此基于应力后更新(update stress last，USL)求解格式的物质点法显式时间积分原理[10-11]的具体实现过程可表述如下。

1)将物质点所携带的质量、动量和速度映射到网格节点上，从而得到网格节点的质量、动量和速度。该步骤的具体表达式分别为：

(1)

(2)

(3)

2)计算网格节点力。网格节点合力的表达式为：

(4)

3)计算下一时间步的节点动量。该步骤的表达式为：

(5)

式中Δtk表示时间步长。

4)将下一时间步的节点动量映射回物质点，并更新物质点的速度和位置坐标。该步骤的表达式分别为：

(6)

(7)

5)计算物质点的旋率增量和应变增量。旋率增量和应变增量的表达式分别为：

(8)

(9)

6)更新物质点的密度和几何尺寸。该步骤的表达式分别为：

(10)

(11)

7)更新物质点的应力。该步骤的表达式为：

(12)

8)重复上述步骤1)～7)，直至达到终止计算时间。

2 物质点法的改进

2.1 构建物质点自由流出边界条件

2.1.1 物质点自由流出边界条件的作用与构建方式

根据物质点法的显式时间积分原理可知，背景网格在每次时间积分计算结束后都会自动复原，并继续参与下次时间积分计算。因此物质点法的背景网格区域始终处于非严格封闭状态，当物质点扩散至背景网格边界时，物质点及其所携带的物理信息都会受到边界节点的约束，而无法从背景网格边界自由流出，导致数值计算结果产生巨大误差。所以一般应设置容积较大的背景网格区域，以尽量减小边界节点约束对数值计算的影响。但是背景网格区域设置过大，又势必会加大数值计算的工作量，导致计算成本大幅度提高，这是物质点法所面临的一种局限。为了破除物质点法的这一局限，本文将基于物质点法显式时间积分原理来构建物质点自由流出边界条件，以使物质点及其物理信息能够根据实际情况从背景网格边界自由流出，从而在保证数值计算精确度的同时提高计算效率。下面对物质点自由流出边界条件的构建方式进行具体介绍。

在每次时间积分计算开始前，都要对模型中的全部物质点进行位置坐标的重新确定，然后按照

(13)

2.1.2 算例分析

为了检验物质点自由流出边界条件的有效性，本文将该边界条件应用到近物面水下爆炸算例中，并通过算例的对比分析结果来对该边界条件的实际效果以及钢板前后的正向冲击波(入射波)、反射波和透射波的形态进行量化描述。在近物面水下爆炸算例中，采用基于Fortran语言自编的物质点法程序来构建3个近物面水下爆炸二维模型，这3个模型在该算例中分别简称为模型A、模型B和模型C，这3个模型的初始形态如图1所示，其中，模型A和模型B的矩形背景网格区域的面积均为100.0 cm×100.0 cm，模型A施加物质点自由流出边界条件，模型B不施加该边界条件；而模型C的矩形背景网格区域的面积为200.0 cm×200.0 cm，且不施加物质点自由流出边界条件。除此之外，这3个模型的其他方面都是相同的，即：直角坐标系的坐标原点均设置在背景网格区域的形心位置，网格单元尺寸均为0.5 cm×0.5 cm，背景网格区域中均填加了面积为100.0 cm×100.0 cm的矩形物质点区域，物质点区域的中心部位构建了尺寸为8.0 cm×8.0 cm的矩形TNT炸药，炸药右侧构建了厚度为2.0 cm的钢板，钢板的前壁面(迎爆面)与炸药形心的距离为250 cm，物质点区域的其他部位由水填充，离散物质点的尺寸为0.25 cm×0.25 cm，在物质点区域中从左至右依次设置了4个观测点，这4个观测点的坐标见表1。

通过图1可知，模型C在其物质点区域的四周均预留出宽度为50 cm的背景网格扩充区域，从而可给物质点膨胀扩散提供充足的空间，这是应用物质点法构建数值模型的常规方式，所以模型C的计算结果可以认定是正确的，因此本算例采用模型C的计算结果来作为模型A和模型B的计算结果的正确性检验标准；而模型A和模型B的背景网格区域的尺寸与物质点区域的尺寸相同，并没有给物质点膨胀扩散提供额外的网格空间，因此可通过将这2个模型的计算结果与模型C的计算结果进行相互比较，以检验施加于模型A的物质点自由流出边界条件的有效性。在分别对这3个模型进行数值计算后，得到了这3个模型中4个观测点的压力时程数据以及这3个模型在各时刻的压力云图。

图1 3个模型的初始形态Fig.1 Initial form of 3 models

现分别提取距离背景网格左侧边界5 cm的观测点1和距离背景网格右侧边界5 cm的观测点4的压力时程数据，绘制出图2所示的这2个观测点的压力时程曲线。

表1 近物面水下爆炸模型的观测点坐标

通过观察图2可知，由于模型C给物质点膨胀扩散提供了网格空间，所以模型C中水下爆炸形成的正向冲击波以及钢板前后形成的反射波和透射波都能够跟随物质点的膨胀扩散而被自由释放，因此在450 μs的计算时间内，模型C中这2个观测点的压力时程曲线都只具有1个波峰，其中观测点1的波峰是由正向冲击波形成的，观测点4的波峰是由板后透射波形成的；而模型B既没有给物质点膨胀扩散提供网格空间，也没有施加物质点自由流出边界条件，这导致模型B中水下爆炸形成的正向冲击波以及钢板后方的透射波在传播至背景网格边界后不能被自由释放，并且在网格边界产生了反射波，因此模型B在这2个观测点位置的压力时程曲线都具有2个波峰，其中第2个波峰都是由网格边界产生的反射波所造成的，这与模型C的压力时程曲线存在显著差异；而模型A虽然没有给物质点膨胀扩散提供网格空间，但是由于在背景网格边界施加了物质点自由流出边界条件，所以模型A中水下爆炸形成的正向冲击波以及钢板后方的透射波都能够从网格边界自由释放，不会在网格边界产生反射波，因此模型A能够得到与模型C十分相近的数值计算结果，从而使模型A和模型C在这2个观测点位置的压力时程曲线高度重合。

图2 3个模型中观测点1和观测点4的压力时程曲线Fig.2 Pressure history of observation point 1 and 4 in 3 models

再通过观察图3所示的3个模型第300 μs的压力云图可知，由于模型C给物质点膨胀扩散提供了网格空间，因此在模型C的压力云图中，能够清楚的观察到物质点区域的膨胀变形过程和压力释放过程；而模型B既没有给物质点膨胀扩散提供网格空间，也没有施加物质点自由流出边界条件，因此模型B水下爆炸形成的正向冲击波以及钢板后方的透射波在传播至背景网格边界后会出现明显的反射现象，反射波压力峰值会达到0.3 GPa以上；而模型A虽然也没有给物质点膨胀扩散提供网格空间，但是由于在背景网格边界施加了物质点自由流出边界条件，致使模型A中的部分物质点及其物理信息能够从网格边界自由流出，所以模型A水下爆炸形成的正向冲击波以及钢板后方的透射波在传播至网格边界后依然能被成功释放，因此模型A的物质点区域内的压力分布与变化情况能够与模型C基本保持一致。

通过以上对比分析结果可知，由于模型B未施加物质点自由流出边界条件，因此模型B的背景网格边界是一种非严格封闭边界，物质点及其所携带的物理信息不能从网格边界自由流出。而模型A在背景网格边界施加了物质点自由流出边界条件，因此物质点及其物理信息能够从网格边界自由流出，从而使模型A能够得到与模型C非常相近的数值计算结果，并且由于模型A的网格节点数量较少，物质点又能够在计算过程中从网格边界大量流出，因此模型A的计算效率较模型C有明显提高，计算结果显示，模型A的计算耗时约为32 869 s，模型C的计算耗时约为50 316 s，模型A的计算耗时仅为模型C的0.653倍。由此可知，施加物质点自由流出边界条件的物质点模型能够在保证计算精确度的同时大幅度提高计算效率，因此物质点自由流出边界条件是十分有效的。

此外，为了进一步观察近物面水下爆炸模型中正向冲击波以及钢板前后的反射波和透射波的基本形态，现分别提取模型A中距离钢板前壁面8 cm的观测点2和距离钢板后壁面8 cm的观测点3的压力时程数据，绘制出图4所示的这2个观测点的压力时程曲线。

图3 3个模型第300 μs的压力云图Fig.3 Pressure contours regarding 3 models at the 300th μs

通过观察图4可知，观测点2的压力时程曲线具有2个峰值超过0.7 GPa的波峰，这2个波峰分别是由水下爆炸的正向冲击波和钢板前壁面反射波所造成的；而观测点3的压力时程曲线则只具有1个峰值超过0.25 GPa的波峰，该波峰是由钢板后壁面透射波所造成的。通过观察图4还可知，观测点2位置的正向冲击波和板前反射波的峰值较高，但正压作用时间较短；而观测点3位置的板后透射波的峰值较低，但正压作用时间较长，这与近物面水下爆炸的真实情况相符合。

下面分别选取图5所示的模型A第130 μs和第170 μs的压力云图来对该模型中正向冲击波以及钢板前后的反射波和透射波的基本形态做更为直观的观察。

通过观察图5可知，在模型A中，水下爆炸所形成的正向冲击波作用到钢板后，会在钢板的前壁面和后壁面分别产生反射波和透射波，且正向冲击波和反射波的压力峰值要高于透射波的压力峰值，同时钢板会发生明显变形，这也符合近物面水下爆炸的真实情况。

图5 模型A第130 μs和第170 μs的压力云图Fig.5 Pressure contours regarding model A at the 130th μs and the 170th μs

2.2 改进物质点自适应分裂方案

2.2.1 物质点自适应分裂改进方案的作用与原理

在标准物质点法中，物质点数量在数值计算过程中是保持不变的，这会导致在计算对象沿某一方向膨胀或拉伸变形的过程中，该方向的物质点也会随之变得稀疏，这种情况会对计算结果产生一定程度的影响，甚至会使数值计算产生严重的拉伸不稳定现象，这是物质点法所面临的另一种局限[10,12]。为了解决这一问题，Ma等[12]设计并提出了物质点自适应分裂方案(adaptive particle splitting scheme)，但是该方案只能使物质点在每次时间积分计算中沿某一个坐标轴方向进行自适应分裂，而不能使物质点在每次时间积分计算中沿多个坐标轴方向同时进行自适应分裂，所以基于该方案的物质点法在多维度物质点模型数值计算中仍然具有局限性。为了彻底破除物质点法的这一局限，本文对物质点自适应分裂方案进行了改进，使物质点在每次显式时间积分计算中能够沿多个坐标轴方向同时进行自适应分裂，以满足多维度物质点模型的数值计算要求，并由此形成物质点自适应分裂的改进方案。下面以二维问题为例，来对物质点自适应分裂改进方案的基本原理进行详细介绍。

在由直角坐标系所表征的二维问题中，计算对象具有沿x轴方向和y轴方向进行膨胀或拉伸变形的可能性，因此在物质点自适应分裂改进方案中，应将物质点自适应分裂分为以下3种情况来分别进行讨论。

1)只沿x轴方向进行分裂。

当物质点沿x轴方向发生膨胀或拉伸变形时，如果物质点在x轴方向上的特征尺寸达到分裂判定条件：

(14)

(15)

2)只沿y轴方向进行分裂。

当物质点沿y轴方向发生膨胀或拉伸变形时，如果物质点在y轴方向上的特征尺寸达到分裂判定条件：

(16)

3)同时沿x轴方向和y轴方向进行分裂。

当物质点沿x轴方向和y轴方向发生膨胀或拉伸变形时，如果物质点在这2个方向上的特征尺寸同时达到式(14)和式(16)所示的分裂判定条件，则使物质点沿这2个方向同时进行均等分裂，分裂以后将产生4个新物质点，这4个新物质点的各项广延量的数值为原物质点的0.25倍，而各项强度量的数值则保持不变，因此新物质点的各项物理量的取值方式为：

(17)

以上即为二维条件下的物质点自适应分裂改进方案的基本原理。而三维条件下的物质点自适应分裂改进方案的基本原理也与之类似，只是需要考虑物质点沿z轴方向进行自适应分裂的情况，因此存在7种分裂方式，此处不再逐一赘述。

2.2.2 算例分析

为了检验物质点自适应分裂改进方案的有效性，本文将该改进方案应用到聚能射流算例中，并通过算例的对比分析结果来对该改进方案的实际效果进行量化描述。在聚能射流算例中，同样采用自编的物质点法程序构建了图6所示的聚能射流模型(本算例中简称物质点聚能射流模型)，该模型的背景网格区域的面积为25.0 cm×10.0 cm，背景网格区域中构建了药型罩锥角为60°的射流弹，射流弹的药型罩的材料为铜，壳体材料为钢，内装药为PBX9501炸药，背景网格区域的其他部位用空气填充，射流弹各结构的物质点的尺寸为0.05 cm×0.05 cm，空气物质点的尺寸为0.1 cm×0.1 cm，网格单元的尺寸为0.1 cm×0.1 cm，整个模型共离散出29 440个物质点和25 000个网格单元，背景网格边界施加了物质点自由流出边界条件，在背景网格区域内沿着射流路径依次设置了6个观测点，这6个观测点的坐标见表2。

图6 聚能射流模型Fig.6 Shaped charge jet mode

表2 聚能射流模型的观测点坐标

此外，由于Autodyn软件是目前公认的数值模拟爆炸与冲击动力学问题效果最好的商业软件之一，因此本算例运用Autodyn软件的Euler-2D Multi-material求解器来构建与物质点聚能射流模型相对应的Autodyn二维轴对称聚能射流模型(本算例中简称Autodyn聚能射流模型)，并且采用Autodyn聚能射流模型的计算结果来作为物质点聚能射流模型的计算结果的正确性检验标准。Autodyn聚能射流模型的Euler网格单元尺寸为0.025 cm×0.025 cm，网格边界添加Autodyn自带的流出边界条件(Flow_Out边界条件)。

然后分别对物质点聚能射流模型和Autodyn聚能射流模型进行数值计算，其中，物质点聚能射流模型的数值计算分别在启动物质点自适应分裂子程序和不启动该子程序这2种条件下进行，以检验物质点自适应分裂改进方案的有效性。数值计算完成后，得到图7所示的第35 μs模型外观形态对比结果以及表3所示的聚能射流头部到达6个观测点位置时头部速度的对比结果，且为了便于观察聚能射流的外观形态，空气介质在图7的3张图片中均不予显示。

图7 3个模型第35 μs的外观形态对比结果Fig.7 Comparative results of the visual forms regarding 3 models at the 35th μs

首先通过观察图7可知，在未启动物质点自适应分裂子程序的物质点聚能射流模型数值计算中，随着内装药爆轰产物的持续扩散和铜质药型罩的持续挤压拉伸变形，内装药爆轰产物物质点和药型罩物质点在背景网格区域内的空间分布密度逐渐降低，数值计算过程产生了严重的拉伸不稳定现象，以至于当物质点聚能射流模型运行至35 μs时，内装药爆轰产物区显露出大面积空白区域，并且药型罩物质点之间也发生了明显的非物理性分离，因此所形成的聚能射流并不连续，这与Autodyn聚能射流模型的聚能射流外观形态存在本质性差别；而在启动物质点自适应分裂子程序的物质点聚能射流模型数值计算中，随着内装药爆轰产物的持续扩散和药型罩的持续挤压拉伸变形，内装药爆轰产物物质点和药型罩物质点在x轴方向和y轴方向上不断进行自适应分裂，从而可保证物质点在背景网格区域内的空间分布密度不会明显降低，有效避免了拉伸不稳定现象的产生，因此当物质点聚能射流模型运行至35 μs时，共产生出84 600个物质点(其中包括29 440个初始物质点)，此时的物质点总数约是初始物质点总数的2.87倍，内装药爆轰产物区没有出现大面积空白区域，并且药型罩物质点之间也没有发生非物理性分离，因此所形成的聚能射流的连续性非常好，这与Autodyn聚能射流模型的聚能射流外观形态具有较高的吻合度。

再通过表3可知，物质点聚能射流模型在未启动物质点自适应分裂子程序的条件下，其与Autodyn聚能射流模型只是在观测点1位置的射流头部速度相对差值的绝对值低于4%，而在其余5个观测点位置的射流头部速度相对差值的绝对值均高于4%；而物质点聚能射流模型在启动物质点自适应分裂子程序的条件下，其与Autodyn聚能射流模型在6个观测点位置的射流头部速度相对差值的绝对值均低于2.5%。由此可知，启动物质点自适应分裂子程序条件下的物质点聚能射流模型与Autodyn聚能射流模型的射流头部速度相差较小，即该条件下的物质点聚能射流模型的计算精确度更高。

表3 3个模型在6个观测点位置的射流头部速度对比结果Table 3 Comparative results of shaped charge jet tip velocity by the location of 6 observation points regarding 3 models

通过本算例的上述对比分析结果可知，物质点自适应分裂改进方案能够使物质点在每次显式时间积分计算中沿多个坐标轴方向同时进行自适应分裂，从而可避免物质点模型的数值计算产生拉伸不稳定现象，进而使数值计算的精确度得以显著提高。由此可知，物质点自适应分裂改进方案对物质点法的改进效果是十分显著的。

3 结论

1) 物质点自由流出边界条件能够使物质点及其所携带的物理信息在每次显式时间积分计算中根据实际情况从背景网格边界自由流出，所以构建背景网格时并不需要给物质点膨胀扩散提供额外的网格空间，因此可使物质点模型在保证计算精确度的同时大幅度提高计算效率；

2) 物质点自适应分裂改进方案能够使物质点在每次显式时间积分计算中根据分裂判定条件沿多个坐标轴方向同时进行自适应分裂，从而有效避免了数值计算过程中拉伸不稳定现象的发生，使多维度物质点模型的计算精确度得以显著提高；

3) 经过这2项改进的物质点法能够更加适用于求解爆炸与冲击类工程实际问题。