陈鑫，高轩能――炸药近地爆炸的数值模拟及影响参数分析.doc

炸药近地爆炸的数值模拟及影响参数分析 陈鑫 1,2，高轩能 1 (1. 华侨大学土木工程学院，福建 厦门 361021；2. 福建农林大学交通与土木工程学院，福建 福州 350002) 摘 要：为了研究炸药近地爆炸冲击波的反射传播过程，应用 ANSYS/LS-DYNA 有限元程序建模，在验证模型及参数选取 正确可靠的基础上，研究了地面材料、炸药高度、空气域形状以及炸药当量等参数对马赫波超压的影响。结果表明：对于研 究非近距离的爆炸作用时，不同地面上的马赫波超压峰值相差不大，可简化为刚性地面；建模宜选择刚性壳体地面以及长方 体空气域；与经验公式相比，数值模拟的马赫波超压峰值误差随着炸药当量的增大而减小，随着炸药高度的减小而增大。 关键字：爆炸；马赫反射波；数值模拟；参数分析；TNT 炸药；地面刚度 中图分类号：O383 文献标识码：A Numerical Simulation and Analysis of Influence Parameters for Explosions near Ground CHEN Xin, GAO Xuan-Neng (1. College of Civil Engineering, Huaqiao University, Xiamen 361021, China; 2. College of Traffic and Civil Engineering, Fujian Agriculture and Forestry University, Fuzhou, 350002, China) The motivation of this paper is to study the shock wave reflection propagation process of explosions near ground, and a simulation model is established by employing finite element program ANSYS/LS-DYNA. Based on the proof of correctness and reliability of the model and parameters selection, the influences of different parameters, such as ground material, explosive height, air domain shape and explosive equivalent on the Mach wave overpressure are investigated. The numerical results show that: 1) Investigating the explosion of non-close range, the Mach wave overpressure peaks are similar on different grounds that can be simplified as rigid ground. 2) Rigid shell ground and cuboid air domain shall be chosen for establishing model. 3) Compared with empirical formula, the numerical simulation error of Mach wave overpressure peak decreases with the increase of explosive equivalent, and increases with the increase of explosive height. Abstract： Keywords: explosions; Mach reflection waves; numerical simulation; parameter analysis; TNT explosive; ground stiffness 近年来，针对各类建筑结构的恐怖爆炸袭击以及在工业与民用建筑中的意外爆炸事件时有发生，给社 会的和谐稳定和公民的生命财产安全造成了重大威胁。此类爆炸多为近地爆炸，爆炸冲击波经过地面反射 加强，比空爆更具破坏力。如何有效准确的模拟出近地爆炸中冲击波的反射传播过程，从而更精确的分析 近地爆炸对建筑结构以及人员的冲击作用，就显得尤为重要。许多学者[1-3]就近地爆炸对建筑结构产生的冲 击波超压，动力响应以及破坏模式进行了数值模拟，由于爆炸模拟的复杂性，在建模中通常将地面简化为 刚性体[3, 4]，未考虑不同材料地面反射冲击波能力的差异性。汪维[5]模拟了刚性地面上的爆炸，获得了不同 爆炸高度下的冲击波变化规律。成凤生[6]研究了刚性平面上方的传播反射规律，研究表明当比例距离大于 1.5 m/kg1/3 时，可以忽略炸药形状的影响。杨亚东[7]，Wang J [8]，李利莎[9]在数值建模中考虑实际的地面材 质，但其研究重点是炸药对土体或混凝土破坏成坑情况，文献[8]还对比了实测和模拟的土体上方 30cm 及 70cm 位置的超压峰值。刘伟[10]试验并模拟了 TNT 近地面爆炸的冲击波的传播，数值模拟结果与试验数据 吻合较好，建模中考虑了地面土壤材料，但未与刚性地面进行对比。李鑫[11]对比模拟了空中爆炸，刚性和 沙土地面上 TNT 爆炸对超压分布的影响进行了简单对比，证明刚性地面上起爆的超压最大。综上所述， 目前的文献中对将各种地面简化为刚性地面的适用性探讨不多，对各种近地爆影响因素的分析比较不够全 面。本文应用 ANSYS/LS- DYNA 有限元程序，模拟炸药近地爆炸建立 TNT 炸药爆炸的数值计算模型，研 究马赫波的传播特性，分析探讨不同参数对马赫波超压的影响。 -1- 1 刚性地面上炸药近地爆炸的数值模拟 1.1 计算模型的建立 计算模型模拟炸药在刚性地面上空爆炸，利用对称性建立 1/4 模型，以节约计算成本。1/4 的炸药尺寸 为 0.2m×0.2m×0.4m，炸药中心距离地高度 1.2m。空气域尺寸取 12m×12m×8m（长×宽×高） ，如图 1 所 示，地面以上 6m，地面以下 2m。炸药与空气选用 SOLID164 实体单元并采用多物质 ALE 网格。刚性地 面选用 SHELL163 壳体单元并采用拉格朗日网格，两种网格通过流固耦合算法相互作用，并在关键字中通 过*CONSTRAINED_LAGRANGE_IN_SOLID 约束来实现。在模型的 Oxz，Oyz 面上施加对称约束，其余 表面均采用透射边界以模拟无限空气域。模型中所有单元网格尺寸均取 0.2m。考虑到在后续不同地面的建 模中，混凝土和土体须采用实体单元模拟，因此建立了刚性实体单元地面的模型进行比对，在原模型中删 除刚性壳体单元，并将地面以下（-2~0m）的空气单元置换成刚性地面单元，并采用了共用节点算法考虑 材料间的相互作用。 空气 炸药 刚性壳体地面 图 1 计算模型 图 2 马赫波超压数值模拟与经验公式对比 Fig.2 Comparison of the Mach wave overpressure with numerical simulation and empirical formulas Fig.1 Computational model 1.2 材料参数取值 炸 药 的 爆 炸 采 用 高 能 炸 药 材 料 模 型 *MAT_HIGH_EXPLOSIVE_BURN 材 料 模 型 ， 并 用 JWL （Jones-Wilkins-Lee）状态方程（式 1）来描述，并用*INITIAL_DETONATION 语句控制炸药单点起爆， 起爆位置即为炸药中心，坐标（0，0，1.2m）。炸药的爆速 D 取 6718 m•s-1，爆压 PCJ 为 18.5 GPa，炸药 密度 ρ=1630 kg•m-3， p = A (1 - wE0 w w - RV - R V )e + B (1 )e + R1V R2V V 1 (1) 2 式中 p 为压力，A，B，R1，R2，ω 为 JWL 状态方程参数，其值由试验确定，V 为相对体积，E0 为初 始内能。JWL 状态方程参数取值如表 1 所示。 参数 取值 表 1. 炸药的材料参数[4] Tab 1. Material parameters of dynamite B/GPa R1 R2 ω 9.4 4.5 1.1 0.35 A/GPa 540.9 E0 / J•m-3 8×109 V0 1.0 表 2. 空气的材料参数[4] Tab 2. Material parameters of air 参数 取值 ρ0/ kg.m-3 1.290 C0 0 C1 0 C2 C3 0 0 C4 0.4 C5 0.4 C6 E/ J.m-3 V0 0 5 1.0 2.5´ 10 空气采用*MAT_NULL 材料模型，并用线性多项式方程*EOS_LINEAR_POLYNOMIAL 来描述，其表 达式如下： p = C 0 +C1m+C 2 m2 +C 3 m3 + (C 4 +C 5 m+C 6 m2 ) E -2- (2) 式中 μ = ρ / ρ0－1，ρ 为质量密度，ρ0 为参考质量密度，E 为单位体积的材料内能，C0，C1，C2，C3， C4，C5，C6 为实常数。参数按理想气体状态方程取值，如表 2 所示。表中 V0 表示初始相对体积。 刚性地面采用*MAT_RIGID 材料模型来控制，参数 CON1 取 7，表示约束 x，y 和 z 三个方向位移， 参数 CON2 取 7，表示约束 x，y 和 z 三个方向转角。 1.3 数值结果分析 本文在地表上距离炸药中心投影点 2~11m 的区间内，每隔 1m 选取一个测点，读取超压数值。这些测 点都位于马赫反射区，为了验证模型计算的准确性，将测点的超压值与经验公式计算值进行对比。马赫反 射区内空气冲击波超压的计算公式[12]为： DPm =DPG (1+cosj ) (3) 式中，ΔPm 为马赫波阵面的超压，φ 为马赫反射区内测点对应的入射波的入射角。ΔPG 为相应的地爆炸超 压，本文取如下 4 种经验公式分别计算（单位：0.1MPa， ）： 1.06 4.3 14 我国爆破安全规程 GB6722-2003[13]： DPG = + 2+ 3 R R H £ 0.35,1£ R £10~15 3 W R (4) 式中 R =R / 3 W 为比例距离（下同），R 表示炸药中心到测点的距离，单位：m，W 为 TNT 炸药当量，单 位：kg（公式(6)和(7) 中用 ton）。 Henrych[14]的空爆公式（对于刚性地面爆炸，计算时取 2 倍 TNT 炸药当量代入）： 2 3 4 DPf =14.0717/ R +5.5397/ R - 0.3572/ R +0.00625/ R 0.05 £ R £ 0.3 2 0.3£ R £1.0 3 DPf =6.1938/ R - 0.3262/ R +2.1324/ R 2 3 DPf =0.662/ R +4.05/ R +3.288/ R 1.0 £ R £10 Newmark[15]的常规炸药地表爆炸的超压峰值公式： DPG = 93 R 32 R 11.568 978.23 曹树鼎[16]通过实验拟合的地爆公式： DPG = + 2 R (5) 6784 + 3 R 0.8 £ R £ 3.0 (6) (7) 由经验公式求得地爆超压后再代入公式（2） ，计算马赫反射波的超压，将数值模拟结果与之相对比， 如图 2 所示，从图中可以看出，马赫波超压峰值衰减的整体趋势均与经验公式保持一致。近爆炸中心处， 与 Henrych 公式最为接近，当 R =0.495 m/kg1/3 时，刚性壳体地面和刚性实体地面与 Henrych 公式之间的误 差分别为 9.1%和 32.4%。而当 R 增大到 0.887 m/kg1/3 后，数值模拟结果与 Newmark 公式更为接近，且随 着 R 的增大，误差开始不断减小。 R =2.35 m/kg1/3 时，刚性壳体和刚性实体地面与 Newmark 公式的误差分 别减小为 3.64%和 1.75%。在 R <1.9 m/kg1/3 的范围内，刚性壳体地面的马赫波超压大于刚性实体地面，与 经验公式更为接近，而当 R >1.9 m/kg1/3 后，刚性实体地面上的超压值开始略大于刚性壳体地面，并更为接 近经验公式。通过对比提取的马赫波的压力时程曲线，发现刚性壳体地面上的马赫波到达时间要早于刚性 实体地面，表明两种地面单元及算法之间存在差异，刚性实体地面上反射波的产生时间相对滞后，与入射 波在合成为马赫波时间亦滞后，并造成了两种地面马赫波超压峰值的不同。在比例较小的区域内，由于爆 炸及与地面相互作用的复杂性，爆炸模拟的超压峰值与经验公式间存在一定的差异，有待进一步研究改进。 2 不同参数变化对爆炸的影响 2.1 不同地面材料的影响 计算模型模拟炸药在不同材料地面上空的爆炸，建模方式与 1.1 节模型基本相同，删除刚性壳体单元， 并分别将地面以下（-2~0m）的网格分别置换成砂质土，粘性土，混凝土三种材料，均采用多物质 ALE 网 格，同时按共用节点算法建模。此外，考虑到实际建筑结构的地面多为土壤地基上浇筑一定厚度的混凝土 层，因此以沙土为地基，在地表位置建立了一层 30cm 厚的混凝土层，此处混凝土采用拉格朗日网格，并 与其他材料按流固耦合算法相互作用。 -3- 炸药和空气材料的状态方程与参数取值均与 1.2 节相同，不再赘述。砂质土和粘性土均采用*MAT_ SOIL_AND_FOAM 材料模型，此材料模型没有应变硬化，屈服极限仅与压力 p 有关，即 12 2 s y =é ë3(a0 +a1 p +a2 p )ù û (8) 参数取值详见表格 4（kg-m-s 单位制） 。表格中 ρ0 为土体的密度（kg/m3） ，G 表示剪切弹性模量，BULK 表示卸载时的体积模量，a0、a1、a2 为式(8)中的塑性屈服常数，PC 为拉应力系数，VCR 是体积沙漏选项， REF 则为是否使用参考几何初始化压力的选项。 表 3. 砂质土和粘性土的材料参数 Tab 3. Material parameters of sandy soil and cohesive soil 参数 ρ0/ kg.m-3 G/Pa BULK/Pa a0 a1 a2 PC VCR REF 砂质土[17] 1850 1800 1.60E7 6.385E7 2.50E8 3.00E10 3.3E11 3.4E9 0 7.033E4 0 0.3 0 -6.90E3 0 0 0 0 T e0 emin 4e6 0.001 0.01 K1 K2 K3 85E9 -171E9 208E9 粘性土[8] 表 4. 混凝土的材料参数[18] Tab 4. Material parameters of concrete 参数 ρ0/ kg.m-3 G/Pa A 取值 2400 14.86E9 0.79 参数 取值 Smax PC/Pa μC 7.0 16E6 0.001 B 1.60 Plock 0.8E9 c 0.007 N 0.61 fC’ 48E6 μlock 0.1 D1 D2 1 0.04 混凝土材料采用*MAT_JOHNSON_HOLMQUIST_CONCRETE 材料模型，即 HJC 模型，是综 合考虑了高压效应、高应变率和损伤的混凝土或岩石类材料模型，其强度按无量纲等效应力表示为[18]： *N * s * =é ëA(1- D) +BP ùé ûë1- cln(e&)ù û (9) 式中 s * 为实际等效应力与准静态单轴抗压强度 fC’之比， P* 为无量纲压力， e&为无量纲应变率，A 为 无量纲粘聚强度，B 为无量纲硬化系数，c 为应变率系数，N 为应变硬化参数，参数取值列于表格 5 中（kg-m-s 单位制）。表格中 T 为无量纲最大拉伸静水压力， e0 为参考应变率， emin 为断裂前最小塑形应变，Smax 为无 量纲最大强度，PC 为压碎压力，μC 为压碎体积应变，Plock 为压实压力，μlock 为压实体积应变，D1、D2 为损 伤常数，K1、K2、K3 为损伤常数。 (b) (a) 图 3 在炸高 1.2m 时的马赫波超压-比例距离曲线 Fig.3 Mach wave overpressure-scaled distance curve when the height of burst equal to 1.2m 同 1.3 节相同，在地表上距离炸药中心投影点分别为 2~11m 的区间内，每隔 1m 选取一个测点，读取 马赫波超压数值。将各种材料的地面与刚性实体地面以及刚性壳体地面的马赫波超压峰值绘图（如图 3） 对比。由图可得，除去采用刚性壳体的超压峰值在 R 较小时有一定差异外，其他采用各种材料的实体单元 地面的超压峰值很接近，特别是到了在 R 较大的 2.35 m/kg1/3 处，各实体地面之间的最大误差仅为 1.99%， 即便在 R =0.495 m/kg1/3 处（地表上距离炸药中心投影点 2m 处） ，砂质土、粘性土、混凝土、土上覆 30cm 混凝土以及刚性实体地面的超压峰值（MPa）分别为 4.730、4.840、4.87、4.87 以及 4.99，误差不大，其中 -4- 两种混凝土地面上的超压峰值与刚性实体地面上的超压很接近（误差 2.4%），说明在本节模型的工况下， 混凝土地面可以认为刚性地面。粘性土地面的误差为 3.0%，砂质土地面的误差最大（5.2%） ，砂质土和粘 性土在炸药作用下的破坏并形成爆炸坑，爆炸能量被消耗，故峰值较小。因此，对于模拟土质地面上的爆 炸，若简化为刚性地面建模，需根据不同的情况，针对性的进行修正。 2.2 不同炸点高度的影响 由前一节的分析可以看到，当炸药高度为 1.2m 时，不同实体单元地面的马赫波超压相差不大，故尝 试将炸药高度进行改变，计算比较马赫波超压峰值的差异。建模选用 2.1 节模型中的砂质土、粘性土、混 凝土、土体上覆混凝土以及刚性实体五种地面情况来比对，保持模型单元、网格和材料参数不变。 数值模拟表明，当炸高较小（0.4m）时，随着炸药成坑体积增大，耗能增加，砂质土、粘性土地面的 冲击波压力场与刚性地面有明显区别，如图 4 所示，在一定时间内，砂质土地面上压力峰值集中在爆炸坑 中，而刚性地面却可以很好的反射冲击波。分别取炸高 0.4m 和 1.6m 对比不同地面的马赫波超压（如图 5 和图 6） ，炸高 0.4m 时，不同材料地面的产生超压差异性增大，在 R =0.433 m/kg1/3 处，刚性实体地面与砂 质土地面间的误差为 11.6%，到了 R =2.34 m/kg1/3 处，这一误差降为 4%。炸高 1.6m 时的超压规律与炸高 1.2m 时相似，并且各种地面的差异进一步缩小，各种地面的超压峰值基本相同，在 R =0.544 m/kg1/3 处， 刚性实体地面与砂质土地面的误差仅为 1.2%。将各种炸高下的超压与经验公式比较和后发现，随着炸高减 小，马赫波超压误差增大，刚性壳体地面上的爆炸马赫波超压与经验公式最为吻合。 (a) 刚性壳体地面 (b) 砂质土地面 图 4 冲击波压力场 Fig.4 Pressure field of shock wave 图 5 炸高 0.4m 时的马赫波超压-比例距离曲线 Fig.5 Mach wave overpressure-scaled distance curve when the height of burst equal to 0.4m 图 6 在炸高 1.6m 时的马赫波超压-比例距离曲线 Fig.6 Mach wave overpressure-scaled distance curve when the height of burst equal to 1.6m 2.3 不同空气域形状的影响 炸药冲击波的传播过程呈放射状，故分别建立球体，圆柱体和长方体的空气域，对应在划分网格时， 球体和圆柱体的网格在曲边上由内向外逐渐变大，（详见图 7），长方体的网格保持为均匀的立方体。建模 方式和参数取值保持不变，地面材料取砂质土、刚性实体和刚性壳体三种情况来比对。对比结果如图 8 所 示，比较可得，对于相同的地面材质，圆柱体和长方体的空气域下数值模拟的马赫波超压很接近，球体空 -5- 气域下的马赫波超压略小。以刚性壳体地面为例， R =2.35 m/kg1/3 处，球体，圆柱体和长方体的空气域的 马赫波超压分别为 1.41 MPa，1.43 MPa 和 1.43MPa。砂质土、刚性实体地面亦可以得到相同的规律。对于 R 较小的区域，采用长方体的空气域下模拟的马赫波超压最接近经验公式。 图 7 球体网格划分 Fig.7 Sphere meshing 图 8 不同空气域的马赫波超压-比例距离曲线 Fig.8 Mach wave overpressure-scaled distance curve with different air domain shape 2.4 不同炸药当量的影响 为了研究不同炸药当量对地爆产生的影响，选用了三种尺寸的立方体炸药进行模拟，炸药边长分别为： 0.2m、0.4m 和 0.6m（即 TNT 当量 W 分别为 13.04kg、104.32kg 和 352.08kg） ，为了减小网格差异对爆炸模 拟产生的影响，在水平投影面上，炸药周围 2m×2m 的区域内采用渐变网格，保证与炸药水平距离 2m 以 外的网格尺寸均为 0.2m，建模方式和参数取值保持不变。在地表上选取 R =1.30m/kg1/3 附近的测点（水平 距离分别为 3m，6m，9m）以及水平距离皆为 10m 的测点并分别将模拟所得马赫波超压与 Henrych 经验公 式进行对比，详见表 4。虽然各个测点的 R 略有差异，但因为随着 R 的增大，误差在不断减小，因而从表 中可得，随着炸药当量的增加，马赫波超压的误差在减小。 同时对比相同 TNT 时，刚性实体和刚性壳体地面上的马赫波超压误差，可以看出总体上，刚性壳体的 误差更小，更适合模拟刚性地面。 表 5. 不同炸药当量的马赫波超压对比 Tab 5. Comparison of shock waves overpressure with different explosive equivalent 炸药当量 kg 13.04 104.32 352.08 13.04 104.32 352.08 水平距离 m m/kg1/3 3.0 6.0 9.0 10.0 10.0 10.0 1.372 1.300 1.286 4.278 2.139 1.426 R 刚性壳体地面 马赫波超压/ 105Pa 与公式误差/ % 3.29 63.4 4.61 42.8 5.65 34.6 0.47 33.3 1.68 39.1 4.42 34.3 刚性实体地面 马赫波超压/ 105Pa 与公式误差/ % 2.28 74.7 3.43 61.5 4.54 47.5 0.46 34.8 1.74 37.0 3.68 45.3 3 结论 通过本文的模拟计算以及影响参数分析表明： （1） 、基于 LS-DYNA 有限元程序实现炸药近地爆炸的数值模拟计算是可行的，模拟结果与经验公式 的误差随着比例距离的增大而减小。 （2） 、除爆炸中心区域外，不同材料地面与刚性地面的差异性较小，混凝土地面的误差最小。综合不 同工况的地爆模拟结果，刚性壳体地面上的马赫波超压峰值误差较小，更适合模拟刚性地面；空气域形状 为长方体时，网格最均匀，马赫波超压误差最小。 （3） 、炸药高度对马赫波超压峰值的影响大，随着炸高减小，误差增大；随着炸药当量的增大，马赫 波超压误差减小。 -6- 参考文献 [1]丁阳, 汪明, 李忠献. 爆炸荷载作用下钢框架结构连续倒塌分析[J]. 建筑结构学报, 2012, 32(2):78-84. [2]高轩能, 刘颖, 杨维英. 马场坪收费站汽车爆炸破坏的数值模拟分析[J]. 振动与冲击, 2012, 31(21):184-189. [3]段雷琳, 高轩能, 江媛. 内爆炸下球面钢网壳结构的冲击波超压计算分析[J]. 华侨大学学报(自然科学版), 2013, 34 (5):557562. [4]高轩能, 王书鹏, 江媛. 爆炸荷载下大空间结构的冲击波压力场分布及泄爆措施研究[J]. 工程力学, 2010, 27(4):226-233. [5]汪维, 张舵, 卢芳云等. 大当量 TNT 装药爆轰的远场数值模拟及超压预测[J]. 弹箭与制导学报, 2010, 30(1):127-130. [6]成凤生, 宋浦, 顾晓辉等. TNT 装药爆炸波在刚性平面上方传播反射的数值研究[J]. 爆破器材, 2011, 40 (4):1-4. [7]杨亚东, 李向东, 王辉等. RDX 基含铝炸药和 TNT 浅层土壤中爆炸开坑的数值模拟与试验[J]. 火炸药学报, 2013, 36(2):24- 29. [8]Wang J. Simulation of landmine explosion using LS-DYNA3D software: Benchmark Work of Simulation of Explosion in Soil and Air[R]. Defence Science &Technology Organisation, 2001.10. [9]李利莎, 谢清粮, 郑全平等. 基于 Lagrange、ALE 和 SPH 算法的接触爆炸模拟计算[J]. 爆破, 2011, 28(1):18-22+27. [10]刘伟, 郑毅, 秦飞. 近地面 TNT 爆炸的试验研究和数值模拟[J]. 爆破, 2012, 29(1):5-9+26. [11]李鑫, 吴桂英, 贾昊凯. 挡墙对爆炸冲击波传播影响的数值模拟[J]. 工程力学, 2012, 29(2): 245-250. [12]李翼祺, 马素贞. 爆炸力学[M]. 北京:科学出版社, 1992: 284. [13]中华人民共和国国家标准. GB6722-2003 爆破安全规程[S]. 北京:中国标准出版社, 2003:43. [14]J. Henrych. 爆炸动力学及其应用[M]. 熊建国译. 北京:科学出版社, 1987:127. [15]Newmark N M, Hansen R J. Design of blast-resistant structures in shock and Vibration Handbook [M]. NewYork:McGraw-Hill, 1961. [16]]曹树鼎. 化爆空气冲击波参数实验数据的经验拟合[J]. 爆炸与冲击, 1984, 4(1):82-85. [17]尚晓江. ANSYS/LS-DYNA 动力分析方法与工程实例. 第 2 版. 北京:中国水利水电出版社, 2008:299. [18]白金泽. LS-DYNA3D 理论基础与实例分析. 北京:科学出版社, 2005:102. -7-