炸药的安全性一直是军事科学重点关注的问题。在处理、制造、存储和运输的过程中, 即使是受到较低速度的撞击, 固体高能炸药也可能被引燃, 甚至升级为更危险的爆炸事故。国内外发展了多种预测炸药受意外撞击的安全性试验, 如Steven^[1]、Susan^[2-3]、落锤^[4]和跌落^[5]等试验。

由于武器装置易受意外撞击, 以实验手段全面预测含能材料的安全性成本较高。随着高科技信息技术的发展, 数值模拟由于成本低, 不受测试技术限制, 成为有效模拟炸药撞击响应特性的重要手段。

炸药的撞击响应特性模拟在20世纪80年代就有相关研究, Boyle等^[6]通过动压剪切装置对炸药的点火阈值做了一维预测。Partom^[7]在Boyle预测的基础上使用了二维模拟, 采用压力和塑性剪切应变率的瞬时乘积值来预测临界点火情况, 通过比较Susan试验和Steven试验得到LX-04炸药撞击点火的阈值速度, 但是压力和应变率的瞬时乘积并不能全面预测炸药的安全性。Browning^[8-9]和Scammon等^[10]根据恒定热流量假设, 结合点火时间、宏观压力和塑性剪切应变率, 得到炸药撞击点火的幂次函数, 时间项使数值模拟数据更有效印证了试验结果。Gruau等^[11-12]在Browning判据的基础上, 采用类混凝土本构模型, 对Steven试验进行数值模拟, 由于本构模型不能完整地表现出炸药的力学特性, 导致数值计算结果与实测数据有一定误差, 但是基本能够满足炸药安全性预测的要求。Ma等^[13]使用有限元分析软件, 以应力和应变率乘积的积分形式表征了PBX炸药的撞击点火阈值。以上研究从压力、应变率和加载时间等方面预测炸药的安全性, 但均基于有限元方法(FEM)对炸药的撞击响应特性进行分析, 而FEM在涉及大变形及材料发生塑性流动的情况下计算精度欠佳。针对这一不足, 研究者发展了多种无网格算法, 包括光滑粒子(SPH)法、物质点法(MPM)和离散元法(DEM)等^[14]。SPH法以SPH粒子替代有限元网格, 可以从根本上避免FEM因网格畸变造成的计算不收敛问题。但是, SPH法需耗费大量计算机资源, 计算效率远低于FEM, 在工程运算中受到较大的限制。

本研究结合FEM与SPH法的优点, 采用有限元法与光滑粒子(FEM-SPH)耦合算法对粘性炸药的撞击响应特性进行分析, 对比FEM的计算结果, 验证FEM-SPH耦合算法在模拟粘性炸药受撞击问题中的优越性。

传统FEM在处理大变形问题时有两种方法。第一种为采用侵蚀算法, 即当网格变形较大时删除畸变网格单元, 如图 1所示, 加载过程中3和7节点产生较大位移, 网格变形量达到失效值, 网格将被删除, 但是网格的删除, 不能反映该处的实际行为, 对加载的响应特性影响较大^[15]。

图 1(Fig. 1)

图 1 畸变网格删除过程示意图[16] Fig.1 Sketch of the deleted process for distortion mesh

第二种为采用共用网格节点随机分离技术^[16], 即对不同网格中相同位置的节点建立节点组元, 设置各个节点分离的失效应变值, 当网格发生较大变形时(达到设置的失效值), 网格间的节点发生分离, 网格分离过程如图 2所示, 但是这样处理同样避免不了单个有限元网格因畸变而导致计算不收敛的问题。

图 2(Fig. 2)

图 2 节点分离过程示意图 Fig.2 Sketch of nodes separation process

在FEM-SPH耦合算法中, SPH粒子与有限元单元面是采用主从节点耦合的算法, 即将SPH粒子看成是从节点, 将有限元节点作为主节点, 以点面算法来处理接触力。对于每一条位于耦合界面的单元边, 根据与它关联的粒子的应力确定一个等效的面力, 并认为其均匀作用在单元边上, 粒子对单元节点作用力示意图如图 3所示。

图 3(Fig. 3)

图 3 粒子对单元节点作用力示意图 Fig.3 Sketch of the forces on element nodes from particles

为了印证FEM-SPH耦合算法能够更准确地模拟炸药类颗粒材料的撞击响应特性, 采用不同算法分别对跌落试验及Steven试验进行数值分析。

跌落试验装置详细参数见文献[5], 对模型进行适当简化, 简化后的FEM计算模型及FEM-SPH耦合计算模型如图 4所示。

图 4(Fig. 4)

图 4 跌落试验计算模型 Fig.4 Calculation model of the drop test

Steven试验装置详细参数见文献[17], 对模型进行适当简化, 考虑到炸药及聚四氟乙烯(PTFE)环在撞击过程中形变均较大, 对FEM-SPH耦合算法模型中的炸药及聚四氟乙烯环分别填充SPH粒子, 得到FEM模型及FEM-SPH耦合模型如图 5所示。

图 5(Fig. 5)

图 5 Steven试验计算模型 Fig.5 Calculation model of the Steven test

壳体、撞针及弹丸为Q235钢, 采用Johnson-cook模型, 相关材料参数见表 1^[18]。状态方程参数为C=4569 m·s^-1, S=1.49, γ=2.17^[19]。对炸药和聚四氟乙烯环采用随动塑性模型, 硬化系数选为1.0(材料采用各向同性硬化模型), 具体参数如表 2所示。

表 1(Tab. 1)

表 1 钢外壳、撞针及弹丸材料参数 Tab. 1 Material parameters for steel shell, firing pin and projectile ρ/g·cm-3 G/GPa v A/MPa B/MPa n c m 7.85 78.9 0.3 245 900 0.94 0.0391 0.1515 Note: ρ is the density; G is the shear modulus; v is the Poisson′s ratio; A, B, n, c and m are J-C constants.

表 1 钢外壳、撞针及弹丸材料参数 Tab.1 Material parameters for steel shell, firing pin and projectile

表 2(Tab. 2)

表 2 炸药和PTFE材料参数 Tab. 2 Material parameters for explosives and PTFE material ρ/g·cm-3 E/GPa v σ/MPa β PBX-2[20] 1.84 10.1 0.30 44.85 1.0 Comp.B[21] 1.65 4.1 0.38 60 1.0 PTFE[11] 2.15 0.4 0.463 14 1.0 Note: ρ is the density; E is the elastic modulus; v is the Poisson′s ratio; σ is the yield stress; β is the hardening parameter.

表 2 炸药和PTFE材料参数 Tab.2 Material parameters for explosives and PTFE

在跌落试验有限元节点分离技术的模拟中, 由于炸药单元变形较大, 节点分离后形成较多单独的六面体单元, 这些分离的六面体单元受到持续的载荷作用, 相互之间发生碰撞产生较大的网格畸变, 计算不容易收敛, 较难得到有效撞击响应数据, 而且在网格数与粒子数相等的条件下, 节点分离法计算时长比FEM-SPH耦合算法要长, 所以下面重点对FEM和FEM-SPH耦合算法进行分析。

跌落试验模拟中, FEM、FEM-SPH耦合算法及实验得到的压力曲线和过载曲线分别如图 6和图 7所示。从图 6和图 7可以得出, 试验得到的最大压力出现在撞击后735 μs, 最大压力为0.162 GPa, 最大过载出现在撞击后775 μs, 最大过载为25400 m·s^-2; FEM-SPH耦合算法最大压力出现在760 μs左右, 最大压力为0.156 GPa, 最大压力误差为4%左右, 最大过载出现在撞击后745 μs, 最大过载为23400 m·s^-2, FEM-SPH耦合算法结果与实测值比较吻合; 而FEM得到的最大压力则出现在485 μs左右, 最大压力为0.136 GPa, 最大压力误差为16%左右, 最大过载出现在撞击后570 μs, 最大过载为25010 m·s^-2。FEM最大压力和最大过载出现时间明显比试验结果靠前, 曲线斜率较大, 造成这个误差主要是因为炸药由含能颗粒和少量粘结剂组成, 其在受到撞击的过程中, 内部颗粒发生相对滑动, 表现出一定的流动特性, 而有限元Lagrange网格不能表现出炸药单元间相对滑移的特性, 所以压力上升速率较快, 最大压力和过载持续时间较长, 加载总时长也相对较短。

图 6(Fig. 6)

图 6 试验与数值模拟压力-时间曲线的对比 Fig.6 Comparison of pressure-time curves between experimental and numerical simulation

图 7(Fig. 7)

图 7 试验与数值模拟加速度-时间曲线的对比 Fig.7 Comparison of acceleration-time curves between experimental and numerical simulation

应变对炸药的安全性预测非常重要, 两种算法在计算500 μs时的塑性应变云图如图 8所示。由图 8可以看出, 在炸药受撞击的中心, 两种算法的炸药应变云图较为相似, 但是由于FEM不能有效模拟颗粒材料内部的塑性流动特性, 使得炸药与钢壳拐角处的边界出现应力集中, 此处单元应变较大, 对炸药的安全性预测产生一定的干扰。

图 8(Fig. 8)

图 8 500 μs时FEM和FEM-SPH耦合算法炸药应变云图 Fig.8 Strain contours on explosive of FEM and FEM-SPH coupling algorithm at 500 μs

根据文献[17]的试验数据, 对184.2 m·s^-1和229.6 m·s^-1两种速度下的Steven试验分别进行模拟, FEM、FEM-SPH耦合算法及实验得到的压力数据如图 9所示。从图 9可以看出, 在较高的撞击速度下, FEM-SPH耦合算法压力曲线震荡较大, 这是由于SPH粒子受压后产生相对滑动使粒子受力不均匀, 但是压力曲线总体趋势与试验结果较为吻合。在较高速度的撞击下, FEM得到的压力上升速率较快, 压力保持平稳的时间较长, 与4.1中FEM计算的跌落试验压力趋势基本相似。炸药的安全性不仅与压力的瞬时值相关, 也与压力的时程相关, 对速度为184.2 m·s^-1撞击情况下的压力时程曲线进行积分, 由于在50 μs之后炸药可能发生爆炸(爆燃), 所以对0~50 μs之间的压力曲线进行积分, Steven试验、FEM-SPH耦合算法及FEM得到的压力积分值分别为10, 8.93 GPa·μs和7.11 GPa·μs, FEM-SPH耦合算法与试验的误差为10.7%, FEM为29%, FEM-SPH耦合算法结果更理想。

图 9(Fig. 9)

图 9 不同撞击速度下试验与数值模拟压力时程曲线的对比 Fig.9 Comparison of pressure histories curves between experiment and numerical simulation under different impact velocity

对于在229.6 m·s^-1速度下的Steven试验的模拟, 分别选取FEM和FEM-SPH耦合算法计算50 μs时炸药的应变云图, 如图 10所示, 两种算法得到的炸药应变云图较为相似。但FEM-SPH耦合算法得到的炸药最大应变出现在弹丸撞击的下方, 而FEM得到的最大应变出现在被“挤出”的撞击边缘部分, 这是由于采用连续Lagrange网格, 不能有效地表现出材料内部的塑性流动特性, 导致此处网格单元变形较大。

图 10(Fig. 10)

图 10 50 μs时FEM和FEM-SPH耦合算法的炸药应变云图 Fig.10 Strain contours on explosive of FEM and FEM-SPH coupling algorithm at 50 μs

对未爆试验回收的炸药样品进行分析可知, 炸药受撞击后发生较大塑性流动, FEM-SPH耦合算法模拟的炸药变形特点与实际情况较接近, 如图 11所示。

图 11(Fig. 11)

图 11 试验与FEM-SPH耦合算法模拟现象对比 Fig.11 Comparison between experimental and phenomenon of FEM-SPH coupling algorithm

(1) 在跌落试验的数值模拟中, FEM-SPH耦合算法得到的压力曲线及过载曲线与测量值较接近, 最大压力误差约4%;而FEM得到的压力及过载曲线上升速率均较快, 最大压力的误差约16%, FEM-SPH耦合算法结果更为理想。

(2) 在Steven试验的数值模拟中, 由于SPH粒子间会产生相对滑动, FEM-SPH耦合算法得到的压力曲线虽然震荡较大, 但总体趋势与试验结果吻合较好。FEM-SPH耦合算法及FEM的压力积分值与试验偏差分别为10.7%和29%。FEM-SPH耦合算法与试验未爆炸药的撞击变形特点符合较好。

(3) FEM-SPH方法能够在一定程度上描述颗粒类材料在加载过程中颗粒相互滑移及流动的特性, 同时FEM-SPH耦合算法提高了总体计算效率, 对于模拟材料产生较大塑性流动的问题有一定优势。采用FEM-SPH耦合算法能够更有效地对粘性炸药的安全性预测提供数据支持。