基于SPH方法的聚能射流侵彻混凝土靶板数值模拟*

强洪夫,范树佳,陈福振,刘 虎

(火箭军工程大学动力工程系,陕西 西安 710025)

基于SPH方法的聚能射流侵彻混凝土靶板数值模拟*

强洪夫,范树佳,陈福振,刘 虎

(火箭军工程大学动力工程系,陕西 西安 710025)

在完全变光滑长度SPH(smoothed particle hydrodynamics)方法的基础上，利用F.Ott等提出的修正SPH方法处理在求解多介质大密度问题时的数值不稳定性问题，运用Holmquist-Johnson-Cook本构模型处理混凝土在冲击载荷下的变形和损伤问题，对聚能装药射流侵彻混凝土靶板的过程进行了数值模拟，同时利用LS-DYNA非线性有限元程序进行对比，分析了2种方法得到的混凝土von Mises应力变化、射流头部特定节点处的速度变化及裂纹演变，验证了SPH方法的准确性。分析了另外2种不同尺寸的靶板在射流侵彻作用下的破坏形式，结果符合射流侵彻物理规律，表明该方法适合模拟聚爆炸与冲击等大变形破坏等问题。

爆炸力学;侵彻;光滑粒子动力学;混凝土;破坏

二级串联随进战斗部[1]作为攻击跑道目标的有效手段，在现在战争中发挥着重要作用。其第1级战斗部即采用聚能装药爆破后形成高速射流，对混凝土目标侵彻形成孔道，以便于第2级战斗部跟进，继续侵彻一定深度后爆破，从而对目标造成不可恢复的破坏。聚能射流侵彻混凝土的主要特点是：炸药能量高，聚能效应明显，形成高速射流,混凝土在高速冲击作用下由于其拉压不等的特性其破坏的模式复杂等。因此，对于该过程的数值模拟对于串联随进战斗部的研究与设计具有重要意义。

传统用于计算该过程的方法主要是基于现有商业软件中的成熟算法，如多物质Euler法、Lagrange法以及ALE(arbitrary Lagrangian-Eulerian)法等，这些方法大都采用有限差分法(finite difference method,FDM)、有限元法(finite element method,FEM)或有限体积法(finite volume method,FVM)等基于网格的数值方法进行离散求解，计算中不可避免地出现网格的扭曲和缠绕(针对Lagrange网格离散)或界面追踪复杂、精度较低(针对Euler网格离散)的问题。光滑粒子流体动力学(smoothed particle hydrodynamics, SPH)方法作为一种无网格粒子方法，在对射流进行弹塑性流体动力学计算，对混凝土进行大应变、高应变率的变形计算时,可避免网格重分及算法耦合，因此非常适合此类问题的求解。最早采用SPH方法模拟爆炸的可追溯到J.W.Swegle等[2]，而L.D.Libersky等[3-4]最先将SPH方法运用于高速冲击领域，随后，M.B.Liu等[5-6]采用SPH方法模拟了聚能装药的爆轰过程，Qiang Hongfu等[7]运用F.Ott等[8]提出的修正SPH方法对聚能装药射流过程进行模拟，分析了不同的起爆方式对射流的影响。本文中，拟在强洪夫等[9-10]提出的完全变光滑长度SPH方法的基础上，结合修正SPH方法[8]，处理爆炸和冲击过程中密度和光滑长度变化剧烈的问题以及多介质界面问题。采用SPH方法计算混凝土损伤大多引入Johnson-Cook[11]、SCG[12]等本构模型，与混凝土的特性不符，未能捕捉混凝土在高速冲击下出现的裂纹扩展等细节，且采用SPH方法计算包括炸药的爆轰、药罩的挤压、射流的形成及发展、射流高速侵彻混凝土、混凝土损伤破碎等过程在内的复杂全过程问题未见报道。因此，本文中将Holmquist-Johnson-Cook本构模型[13]引入SPH方法，处理混凝土在高速冲击作用下的变形损伤问题，对聚能炸药爆炸挤压形成射流，进而侵彻混凝土靶板的整个过程进行数值模拟，并对侵彻过程中射流头部特定点处的速度变化进行分析，采用传统网格方法进行对比验证。同时为更充分地描述混凝土靶板在高速射流侵彻作用下的破坏损伤效应，对另外2种不同厚度的靶板进行设计，并开展数值实验，探讨SPH方法用于模拟诸如聚能射流侵彻混凝土等涉及爆炸与冲击大变形多介质问题的适用性。

1 SPH基本方程

在SPH方法中，通过对问题域的粒子离散来求解系统的状态，系统中物理量f(r)及其导数·f(r)的SPH形式为：

f(ri)

(1)

(2)

式中：m、ρ、r分别为粒子的质量、密度和位置矢量;Wij=W(ri-rj,h)为核函数，h为光滑长度，通常选用三次样条核函数[14]。

为更好地模拟爆炸与冲击、大变形大扭曲等密度和光滑长度变化剧烈的问题，本文中采用完全变光滑长度SPH方法[7]，同时为了很好地解决药型罩和爆轰气体间密度梯度较大带来的间断面不稳定的问题，采用修正SPH方法[10]对聚能装药射流过程进行模拟，结合后的方程组如下：

(3)

式中:Wij=W(xi-xj,h)为核函数，它的选取直接影响计算的误差和稳定性，通常选用三次样条核函数；iWij为核函数对粒子i坐标的空间导数，vij=vi-vj；σ为总应力张量；Πij为人工黏度；h为插值核宽度的一种度量，称为光滑长度，表示W不显著为零的取值范围，控制着SPH粒子的影响域，通常设定时，W=0,这里h为空间和时间的函数，在连续方程中将其对时间求导，即：

(4)

式中:d为空间维数。对于光滑长度变化率dhi/dt与密度变化率dρi/dt之间相互耦合而带来的难以显式求解的问题，本文中采用迭代法[9]求解密度方程和光滑长度。

在计算过程中，为使粒子分布更有序，消除由于分布不均匀引来的粒子非物理聚集的问题，采用J.J.Monaghan[14]提出的XSPH对速度场进行光滑:

(5)

式中:ε(0≤ε≤1 )为常数，通过施加邻近粒子的影响使自身的运动速度与邻近粒子的平均速度相近，一般ε取0.3。

2 炸药及药型罩材料模型

本文模拟中，炸药选用TNT，其爆轰速度由实验[15]测得为6 930 m/s。对于爆轰气体,使用标准JWL状态方程[16]得到气体压力:

(6)

式中:ρ0=1 630 kg/m3为TNT炸药的初始密度,v=ρ0/ρ,ρ为爆轰气体的密度,e为爆轰气体的比内能;A、B、R1、R2和w为由实验结果拟合得到的系数[17],且A=371.2 GPa,B=3.231 GPa,R1=4.15,R2=0.95,w=0.3。

对于药型罩材料，选用铜，状态方程选用Mie-Grüneisen状态方程[16]:

(7)

(8)

式中:ρ0=8 530 kg/m3为铜的初始密度，Γ=1.99，cs=3 940 m/s,Ss=1.489。

3 混凝土本构模型及状态方程

图1 混凝土屈服强度模型曲线Fig.1 Curves for yield strength model of concrete

图2 混凝土累计损伤模型曲线Fig.2 Curve for cumulative damage model of concrete

HJC本构模型[13]分为强度模型、累计损伤模型和状态方程3部分，模型曲线分别如图1～3所示。

3.1 屈服强度模型

归一化的屈服强度函数为：

(9)

3.2 累计损伤模型

损伤值为由等效塑性应变和塑性体积应变引起的损伤累积起来的形式：

(10)

(11)

式中:D1和D2为混凝土的损伤常数。累计损伤模型曲线如图2所示。

3.3 状态方程

图3 混凝土状态方程曲线Fig.3 Curves for equation of state of concrete

混凝土靶板在冲击载荷下的压缩变形分为弹性阶段、过渡阶段和压实阶段3个阶段，状态方程曲线如图3所示，图中的plock和pcrush分别为下式中的pl和pc,μlock、μcrush和μp,lock分别为下式中的μl、μc和μpl。

(1)弹性阶段(0≤μ≤μc)

在这一阶段，混凝土发生可以恢复的弹性变形，加载与卸载的方程相同，为：

p=Keμ

(12)

式中:弹性体积模量Ke=pc/μc，pc和μc分别为单轴强度抗压实验中得到的压碎压力和压碎体积应变，μ=ρ/ρ0-1为单元体积应变，ρ和ρ0分别为单元的密度和初始密度。

(2)过渡阶段(μc<μ≤μpl)

在这一阶段中混凝土中的气体逐渐被挤压出去，混凝土产生破碎性断裂，加载时:

p=pc+Kp(μ-μc)

(13)

式中:Kp=(pl-pc)/(μpl-μc)为塑性体积模量，pl为压实压力;μpl=ρg/ρ0-1为压实体积应变，对应于pl处的体积应变,ρg为颗粒密度。

该区的卸载是通过相邻区域间插值的一条路径进行的，卸载的状态方程为：

p=[(1-F)Ke+FKp]μ

(14)

式中:插值因子F=(μmax-μc)/(μpl-μc) ，μmax为卸载前达到的最大体积应变。

(3)压实阶段(μ≥μpl)

该阶段混凝土完全被压碎，加载的状态方程为:

(15)

卸载的状态方程为：

p=K1μ

(16)

混凝土所能承受的最大拉伸强度为T(1-D),状态方程曲线如图3所示。

4 聚能射流侵彻混凝土靶板的数值模拟

数值模拟中采用的聚能装药模型及混凝土模型如图4(a)所示，药柱宽度为40.00 mm，装药头长度为20.75 mm，药孔长度为42.00 mm，张角为60 °,药型罩厚度为3.00 mm，炸高为40.00 mm，混凝土靶板尺寸为100 mm×30 mm。炸药采用点起爆的方式，起爆点为(0,0.07),具体粒子配置如图4(b)所示,粒子间距均为0.5 mm,其中TNT炸药粒子数为26 766,药型罩粒子数为4 480,混凝土粒子数为12 000。光滑长度取1.5倍粒子间距，时间积分采用蛙跳格式，时间步长为0.1 μs。

本文中利用LS-DYNA非线性有限元程序进行对比验证。计算模型由炸药、药型罩、空气和混凝土靶板4部分组成。其中，对炸药、药型罩和空气采用Euler网格建模，对混凝土靶板采用Lagrange网格建模，在混凝土靶板与空气和药型罩材料间采用耦合算法。聚能装药是线性的，考虑到模型的对称性，为了减少计算时间，将模型简化为二维平面对称问题并建立其1/2模型。在对称面上施加对称边界约束，在混凝土靶板边界施加固定约束，把靶板侧面设定为无反射边界。对炸药、药型罩、空气和混凝土靶板均选用与SPH数值模拟中参数相同的本构模型和状态方程。

图4 算例模型结构Fig.4 Construction of the model

图5 聚能射流形成过程对比Fig.5 Comparison of the formation of jet flows

图5给出了采用SPH方法和LS-DYNA程序计算过程中炸药爆轰及药罩被挤压、射流形成过程的对比。图6给出了采用SPH方法和LS-DYNA程序计算过程中射流侵彻混凝土形成开孔、混凝土损伤及混凝土完全破坏的过程的对比。可以看出，炸药爆轰波逐渐将药型罩向中轴线方向挤压形成射流和杵体，随后混凝土在射流的侵彻下发生压缩和剪切变形，产生剥落，靶板弹坑的上部呈漏斗状。随着射流的不断侵彻，靶板背面首先由于反射波到达壁面后的反射拉伸作用，导致背面出现拉伸损伤，部分脱落，开孔周围则逐渐出现不同程度的层裂，并且裂纹逐渐向外扩展，最终形成一个漏斗状损伤区域，中间部分出现粗细均匀的孔道，裂纹形态与图7中的实验结果基本相符。通过对比验证，发现利用2种数值模拟方法得到的结果均符合射流侵彻物理规律，但LS-DYNA程序通过定义失效删除单元产生的裂纹非常规则，数量少，基本没有在裂纹上产生小裂纹的状况，且不能模拟侵彻后碎片的物理特性，裂纹模拟效果受到网格划分大小的影响，SPH方法则恰恰克服了以上困难。

由图6可知，在侵彻过程中，前2段射流发挥了主要作用。因此，分别选取SPH方法及LS-DYNA程序计算结果中前2段射流的头部顶点粒子A、B进行分析。图8为2种方法的计算结果中A点和B点粒子速度随时间的变化曲线，由此可知，药型罩在0.02 ms时受到炸药爆轰波的挤压开始逐渐加速，在约0.10 ms时，在SPH方法的计算结果中A点粒子速度达到最大值3 247.9 m/s，在LS-DYNA程序的计算结果中A点粒子速度达到最大值3 457.4 m/s。射流头部在0.14 ms时开始与靶板接触，速度大幅降低，造成混凝土的损伤，随着侵彻的深入，速度逐渐减小为零。2种方法的不同之处在于，SPH方法的计算结果中第1段射流在0.20 ms时速度基本降为零的同时第2段射流开始侵彻，而在LS-DYNA程序的计算结果中，第1段射流在0.14 ms时与混凝土靶板接触，速度逐步降低，在0.17 ms左右时，第2段追上第1段射流，2段射流合为一段继续侵彻，直到速度基本降至零。

图6 聚能射流侵彻过程中不同时刻，混凝土靶板中Von Mises应力的分布Fig.6 Von Mises stress distribution in concrete target penetrated by shaped charge jet at different times

图7 聚能射流侵彻的混凝土靶Fig.7 Concrete target penetrated by shaped charge jet

图8 射流头部特定点处速度-时间曲线Fig.8 Velocity-time curves at special points of jet head

相对于LS-DYNA中通过删除网格单元形成裂纹的处理方法，SPH不删除粒子单元而是通过显示混凝土靶板的损伤值D来表述材料失效形成的裂纹，损伤后的SPH粒子按无黏性的流体粒子进行计算，可以完整地再现裂纹的损伤演化过程，损伤状态及裂纹扩展分布情况如图9所示，计算得到的靶板最终损伤状态如图10所示。由图10可知：采用LS-DYNA得到的计算结果中弹坑直径与采用SPH方法得到的计算结果中弹坑直径差别较大，其主要原因在于LS-DYNA中理想空气与混凝土靶板的重叠区域小；采用SPH方法得到的计算结果中靶板左右壁面出现裂纹，而采用LS-DYNA得到的计算结果中没有出现，其主要原因在于SPH中靶板的左右边界是自由边界，冲击波到达左右壁面时由于反射波的作用开始出现拉伸损伤，而LS-DYNA靶板定义的是无反射边界，反射波不能到达左右壁面。

图9 在聚能射流的侵彻下不同时刻混凝土靶板的损伤情况Fig.9 Damage in concrete target penetrated by shaped charge jet at different times

图10 不同方法模拟得到的裂纹长度及分布情况Fig.10 Length and distribution of crack simulated by different methods

为了更充分地描述混凝土靶板在高速射流侵彻作用下的破坏损伤效应，采用上述SPH方法对另2种尺寸的靶板进行了数值实验，粒子数分别为4 000和40 000，计算结果见图11，图12为选取的射流头部特定点进行追踪得到的速度-时间曲线。从图11～12可以看出，当靶板尺寸较小时，在较短的时间内只需射流头部一段即可将靶板完全击穿，射流孔处损伤均匀；靶板较厚时，需要3段或更多段射流才可将靶板击穿，而且在射流孔处不会形成大面积的损伤，而是形成不同程度的损伤裂纹延伸至靶板内。

图11 另外2种不同尺寸的混凝土靶板在聚能射流侵彻作用下的损伤Fig.11 Damage in other two concrete targets with different sizes penetrated by shaped charge jets

图12 另外2种不同尺寸的靶板聚能射流侵彻作用下特定节点的速度-时间曲线Fig.12 Velocity-time curves at special points of other two concrete targets with different sizespenetrated by shaped charge jets

5 结 论

运用无网格SPH方法计算得到了聚能装药射流形成、发展及高速冲击混凝土整个过程,并与传统网格进行了对比验证，得知：

(1)采用SPH方法能有效解决在高速冲击、大变形问题的计算中不可避免地出现Euler网格界面追踪复杂或Lagrange网格扭曲和缠绕的问题。基于完全变光滑长度SPH方法与修正SPH方法相结合的方法，有效解决了爆炸与冲击中光滑长度变化剧烈的问题以及多介质界面上大密度间断带来的计算不稳定问题。

(2)运用HJC本构模型和状态方程描述混凝土在高速冲击下的变形损伤特性，不仅得到了破孔的形状，而且对由于壁面反射波造成的靶板背面拉伸损伤及破孔周围裂纹的扩展等变形损伤细节描述精确，与实际物理现象符合。采用2种方法得到的计算结果基本一致，验证了SPH方法的准确性。同时，采用SPH方法计算得到的裂纹形态与实验结果吻合更好，更合理。

(3)随着靶板尺寸的增大，射流孔径逐渐减小，裂纹数目逐渐增加，裂纹的细度及延伸的长度逐渐增加，破坏的能量逐渐向周围扩散。

[1] 段建,杨黔龙,周刚,等.串联随进战斗部侵彻混凝土靶实验研究[J].爆炸与冲击,2007,27(4):364-369. Duan Jian, Yang Qianlong, Zhou Gang, et al. Experimental studies of a tandem follow-through warhead penetrating concrete target[J]. Explosion and Shock Waves, 2007,27(4):364-369.

[2] Swegle J W, Attaway S W. On the feasibility of using smoothed particle hydrodynamics for underwater explosion calculations[J]. Computational Mechanics, 1995,17(3):151-168.

[3] Libersky L D, Petscheck A G. Smoothed particle hydrodynamics with strength of materials[C]∥Trease H, Fritts J, Crowley W. Proceedings of the Next Free Lagrange Conference. NY: SpringerVerlag, 1991,395:248-257.

[4] Libersky L D, Petscheck A G, Carney T C, et al. High strain Lagrangian hydrodynamics: A three-dimensional SPH code for dynamic material response[J]. Journal of Computational Physics, 1993,109(1):67-75.

[5] Liu M B, Liu G R, Zong Z, et al. Computer simulation of high explosive explosion using smoothed particle hydrodynamics methodology[J]. Computers & Fluids, 2003,32(3):305-322.

[6] Liu M B, Liu G R, Lam K Y, et al. Meshfree particle simulation of the detonation process for high explosives in shaped charge unlined cavity configurations[J]. Shock Waves, 2003,12(6):509-520.

[7] Qiang Hongfu, Wang Kunpeng, Gao Weiran. Numerical simulation of shaped charge jet using multi-phase SPH method[J]. Transactions of Tianjin University, 2008,14(1):495-499.

[8] Ott F, Schnetter E. A modified SPH approach for fluids with large density differences[J]. Arxiv Physics E-prints, 2003:3112.

[9] 强洪夫,高巍然.完全变光滑长度SPH法及其实现[J].计算物理,2008,25(5):569-575. Qiang Hongfu, Gao Weiran. SPH method with fully variable smoothing lengths and implementation[J]. Chinese Journal of Computational Physics, 2008,25(5):569-575.

[10] 强洪夫,王坤鹏,高巍然.基于完全变光滑长度SPH方法的HE爆轰过程的数值试验[J].含能材料,2009,17(1):27-31. Qiang Hongfu, Wang Kunpeng, Gao Weiran. Numerical study of high explosive detonation process using SPH method with fully variable smoothing lengths[J]. Chinese Journal of Energetic Materials, 2009,17(1):27-31.

[11] Johnson G R, Cook W H. A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[C]∥Proceedings of the Seventh International Symposium on Ballistics. Hague, Netherlands, 1983:571-574.

[12] Steinberg D J, Cochran S G, Guinan M W. A constitutive model for metals applicable at high strain rate[J]. Journal of Applied Physics, 1980,51(3):1498-1504.

[13] Holmquist T J, Johnson G R, Cook W H. A computational constitutive model for concrete subjected to large strains, high strain rates, and high pressures[C]∥Proceedings of the 14th International Symposium on Ballistics. Quebec, Canada, 1993:591-600.

[14] Monaghan J J. Smoothed particle hydrodynamics[J]. Reports on Progress in Physics, 2005,68(8):1703-1759.

[15] Liu G R, Liu M B.光滑粒子流体动力学:一种无网格粒子法[M].韩旭,杨刚,强洪夫,译.长沙:湖南大学出版社,2005:195-197.

[16] Livermore Software Technology Corporation. LS-DYNA keyword user’s manual[M]. Livermore: Livermore Software Technology Corporation, 2012:17-45.

[17] Liu M B, Liu G R, Zong Z, et al. Computer simulation of high explosive explosion using smoothed particle hydrodynamics methodology[J]. Computers and Fluids, 2003,32(3):305-322.

(责任编辑 张凌云)

Numerical simulation on penetration of concrete target by shaped charge jet with SPH method

Qiang Hongfu, Fan Shujia, Chen Fuzhen, Liu Hu

(PowerEngineeringDepartment,Xi’anHi-TechInstitute,Xi’an710025,Shaanxi,China)

On the basis of the smoothed particle hydrodynamics (SPH) method with fully variable smoothing lengths, the modified Ott-Schnetter SPH method was used to cope with the computational instability in the large density gradient problems. The Holmquist-Johnson-Cook constitutive model was used to cope with the deformation and damage of concrete under impact. By comparing with the corresponding ones simulated by the non-linear finite element program LS-DYNA, the changing courses of the von-Mises stress and cracks in the concrete targets, and the velocities at the some special points of the jet head were analyzed, which proved the feasibility and accuracy of the SPH method. Other two concrete plates with different sizes penetrated by shaped charge jets were also simulated. The results are in good agreement with the physical principle. So this method can be used to deal with the multi-material and large deformation problems such as detonation and impact.

mechanics of explosion; penetration; smoothed particle hydrodynamics; concrete; damage

10.11883/1001-1455(2016)04-0516-09

2014-12-03；

2015-03-24

国家自然科学基金项目(51276192);国家重点基础研究发展计划(973计划)基金项目(61338); 火箭军工程大学创新型基金项目(EPXY0806)

强洪夫(1963— )，男，博士，教授，博士生导师;

范树佳，fan_shu_jia@163.com。

O389国标学科代码：1303530

A