破片飞散仿真SPH方法中的初始粒子间距

2018-08-30周宵灯崔村燕赵蓓蕾

兵器装备工程学报 2018年8期

周宵灯，崔村燕，赵蓓蕾，詹翔，王岩

(航天工程大学 a.研究生管理大队； b.宇航科学与技术系，北京 101400)

战斗部爆炸效应数值计算涉及高压力、高应变率以及材料的破坏等问题，需要模拟边界的畸变、跟踪裂纹的扩展等。在处理这类复杂问题时，传统方法因网格自身约束而造成的缺陷尤为明显。1977年 Lucy[1]和Monaghan[2]最早提出了光滑粒子流体动力学(SPH)方法，该方法最初在天体物理学领域中用来模拟三维无界空间中天体的演化。在此基础上Johnson与Cook通过建立粘塑性破裂模型,成功模拟了薄片在拉力作用下裂纹产生到破片形成的过程[3]，将该方法推广到材料动载响应领域[4-5]。SPH方法作为一种无网格数值方法，其离散化过程不需要划分单元，而是使用固定质量的运动点的积分方程进行求解，是目前求解爆炸冲击问题的主要方法之一。在使用SPH方法进行爆炸破片仿真时，初始粒子间距的选取对仿真结果有重要影响。若取值过小，则在搜索域内不能对质点提供足够的作用力，从而降低计算精度；若取值太大，则质点的详细特征和局部特质可能被抹平，同样会降低精度。为使SPH方法精确可靠，需要合理选取粒子间距。

目前，国内外对SPH方法开展了大量的理论和数值方面的分析研究，孙晓旺等[6]通过比较一维波动方程的SPH格式和有限差格式获得了衡量SPH方法模拟应力波的重要指标；赵亚洲等[7]针对传统的光滑长度自适应的缺陷，构造了一种避免“拖尾”现象的自适应准则；周杰等[8]提出了SPH方法进出口数值边界条件处理方法；徐金中等[9]采用SPH方法对高速碰撞问题做了数值模拟，分析了初始光滑长度和非一致粒子间距对计算结果的影响。这些研究对SPH方法的准确性、稳定性等进行了一定的优化，然而大多数是基于均匀分布的粒子，有的仅仅适用于一维情况，在实际数值仿真中对于初始粒子间距选取大多还是凭借经验，缺乏一个简单有效的选取规律，对于不均匀粒子间距对结果的影响研究较少。

本文以破片飞散仿真为例开展SPH方法中的初始粒子间距对数值模拟影响的研究，考虑均匀初始粒子间距以及不均匀粒子间距两个因素对仿真结果的影响，并与理论结果进行对比分析，为合理选取粒子间距提高计算精度提供参考。

1 SPH方法基本原理

SPH方法的理论基础来源于粒子方法，即通过一系列具有一定质量、坐标、速度、内能的离散粒子的插值近似描述连续的物理量。该方法通过“核函数”的积分核进行“核函数估值”近似，将流体力学基本方程组转换成数值计算的SPH方程。整个流场离散成一系列“粒子”，所有力学量由这些粒子负载。这些粒子可以按计算公式，即按流体力学流动的规律任意流动。可以假设每个粒子都有质量，拉格郎日插值点，速度以及内能，而其它的变量由插值或所构成的相互关系中求得。

常用的SPH离散守恒方程为[4-5]

(1)

(2)

(3)

其中:式(1)为质量守恒方程，式(2)为动量守恒方程，式(3)为能量守恒方程，ρ为密度,m为粒子质量，δ为粒子应力，v为粒子速度，E为比内能，w为核函数。下标i为粒子编号，下标j为i粒子的近邻粒子编号。

SPH方法中近似过程包括核近似和粒子近似。对函数的核近似是通过核函数和权函数积分实现，粒子近似是对某一区域内粒子加权求和得到。针对不同情况，研究人员提出了不同的核函数进行模拟。J.J.Monaghan等[10]在三次样条函数的基础上提出了B样条核函数，目前应用最广泛，其表达式为[4]

(4)

其中z=(xi-xj)/h，h为光滑长度，对于一维问题N=1.5h，对于二维问题N=0.7πh2，对于三维问题N=πh3。

SPH方法循环过程如图1所示。

2 数值仿真与分析

2.1 计算模型与材料参数

本文数值计算的对象为某型战斗部，其底部1/4切面模型如图2所示。战斗部壳体质量M=7 kg，装药TNT质量m=1.27 kg，采用中心点起爆，底部壁厚d=12 mm其余诸元见文献[11]。

mm-mg-ms单位制，材料包括TNT炸药和45钢，TNT炸药爆轰压力P、单位体积内能E、相对体积V之间的关系采用JWL状态方程进行描述：

(5)

式(5)中：A、B、R1、R2、ω为JWL状态方程参数；E为炸药内能；V为相对体积，具体数值[12]如表1所示。

ρ/(g·cm-3)E/(GJ·cm-3)V/(km·s-1)A/GPa1.6366.93374B/GPaR1R2ω3.754.150.90.35

45#钢材料的本构方程采用JOHNSON_COOK材料模型：

σ=(A1+A2εm1)(1+A3lnε*)(1-T*m2)

(6)

式(6)中：σ表示材料等效抗拉强度；ε表示等效塑性应变；ε*=ε/ε0表示无量纲的总塑性应变率；T*=(T-Troom)/(Tmelt-Troom)，Tmelt表示材料的熔点温度，Troom表示参考温度。A1、A2、A3、m1和m2表示与材料有关的实验参数，具体参数[12]如表2所示。

表2 45钢JOHNSON_COOK模型参数

2.2 不同均匀初始粒子间距下的破片飞散结果

定义无量纲参数λ为战斗部壁厚d与初始粒子间距h0之比,即λ=d/h0，固定战斗部壁厚d=12 mm不变，改变粒子间距。图3所示为λ=12、8、6、4、3、2的战斗部底部1/4切面粒子分布模型图，各工况粒子尺寸和数量见表3。

图4给出了不同λ值时战斗部破片生成情况，图5为不同λ值时战斗部爆炸生成破片数量随时间变化曲线，图6表示了仿真结果与实验结果对比。仿真计算结果表明，初始粒子间距选取不同λ值时，因搜索域内对质点提供的作用力不同，直接造成了爆炸生成碎片数量各不相同。生成碎片总数随着λ值增大而增大。已知该型战斗部爆炸后实测破片为1 082片[13]，因此合理的λ取值6和8之间。

λh0/mm粒子数量261 782346 0504314 0156249 40081.5115 710121393 000

由计算破片初速的能量模型[13-14]得

(7)

其中，D表示TNT爆速，β=m/M，m为TNT质量，M为壳体质量。

将战斗部相关参数代入能量模型得到该战斗部破片的理论初速为938.2 m/s，与不同λ值破片初速数值仿真结果对比情况如表5所示。从表中可以发现不同λ值的破片平均初速仿真结果与理论结果拟合较好，但普遍要略小于理论初速，因为在理论计算能量模型的能量守恒方程中，认为战斗部爆炸时空气介质吸收能量和壳体破坏能量占总能量比值小，人为忽略了这两部分能量，使得破片动能较实际情况有所增加，导致理论结果大于仿真结果。

表5 不同λ值破片平均初速数值仿真与经验结果

2.3 非均匀粒子间距的破片飞散结果

定义γ为战斗部壳体粒子间距h1与内部TNT粒子间距h2之比,即γ=h1/h2。为研究非均匀粒子间距对仿真结果的影响，保持壳体粒子间距h1=2 mm不变，共选择以下γ=2、5/3、4/3、1、1/2五种工况：壳体粒子间距h1=2 mm，TNT粒子间距h2=1 mm、1.2 mm、1.5 mm、2 mm、4 mm。不同γ值战斗部破片生成情况如图7所示，图8为五种工况下最终生成的碎片总数，可以看出非均匀粒子间距对计算结果影响较小。

图9粒子数量与计算精度计算时间的关系曲线中，曲线1表示数值模拟计算精度随粒子数量的变化，曲线2代表计算时间随粒子数量的变化。可以看出对于同一个模型，粒子数量较少时增加粒子数量可以使计算精度明显提高，而计算时间不会有太大变化。当粒子数量增加到一定数量后，再继续增加粒子数量计算精度不仅不会提高反而有所下降，而计算时间却大幅度增加。因此，在进行初始粒子间距选取最理想的情况是，对于重要的区域，粒子间距应该合理缩小，对于计算结果影响不大的区域，粒子间距应该有所增大，这样可以提高计算效率。