基于爆能等效原理大型模爆器燃气爆炸冲击加载的数值模拟＊

张秀华，张春巍，段忠东

（1.东北林业大学土木工程学院，黑龙江 哈尔滨150040；2.哈尔滨工业大学土木工程学院，黑龙江 哈尔滨150090）

近年来，随着各类爆炸事故的不断发生，开展工程结构抗爆分析、设计与防护问题研究逐渐成为土木工程及相关领域新的热点研究方向。欧进萍等［1］和张春巍［2］对大型核爆炸模拟器（模爆器）实验装置进行了系统改造，开发了采用乙炔－空气混合气体爆炸产生的爆炸冲击波对结构（构件）进行加载的加载方法。结合该实验系统的加载控制与实施，可采用数值模拟方法研究乙炔－空气混合气体爆炸产生的空气冲击波的传播规律。

气体爆炸是预混燃烧的一种。当爆炸波与周围介质发生相互作用，火焰加速非常快。爆炸包括预混火焰的初始爆燃以及随后更剧烈的爆轰，这取决于边界条件和超压水平。爆炸过程中，由火焰前阵面加速引起的初始超压导致的破坏通常比由随后的火焰引起的破坏严重得多。有关开敞空间气体爆炸的理论研究方法，包括TNT当量法、TNO多能模型、Kuhl模型以及数值模拟方法等。其中，TNT当量法属偏差较大的经验方法。而气体爆炸的数值模拟研究方法，因为能够模拟影响气体爆炸的众多因素，在预测可燃气体爆炸场的特性中应用越来越广泛［3］。本文中，运用LS－DYNA软件对化学炸药和乙炔－空气混合气体爆炸冲击波在空气中的传播进行数值模拟，并将结果与经典的经验公式进行比较，分析两种爆炸源产生爆炸作用的异同，提出气体爆炸的名义比例距离概念计算气体爆炸压力，验证采用乙炔－空气混合气体代替化学炸药作为爆炸源进行结构抗爆实验的可行性。

1 Euler控制方程

LS－DYNA程序具有Lagrange、Euler和ALE算法。在空间离散域描述方面，采用Euler算法时，有限元计算网格是固定不动的，网格节点就是空间点，而材料则在网格间流动，可以避免Lagrange算法经常遇到的网格大变形所导致的计算困难，特别适合分析流体问题。对于空气中的爆炸冲击波传播问题，由于爆轰产物和空气都属于流体，所以相对于Lagrange算法，采用Euler方法描述更合适和有效。质量守恒、动量守恒和能量守恒的Euler控制方程分别为［4］：

式中：ρ是密度，v是速度，σ是Cauchy应力张量，ε是应变率张量，b是体力，e是内能。

2 计算模型

2.1 材料模型和参数

采用高能炸药燃烧材料模型和状态方程描述高能TNT炸药爆轰产物的状态，其形式为［5］：

式中：p为静水压力（以压为正），V 是相对体积，E是炸药的体积初始内能，A1、B1、R1、R2、ω 均为JWL状态方程参数。TNT炸药主要材料参数分别为：密度ρ＝1 630kg／m3，爆速D＝6.930km／s，爆压pCJ＝21.0GPa，A1＝371.2GPa，B1＝3.231GPa，R1＝4.15，R2＝0.95，ω＝0.35，初始内能 E＝7.0GJ／m3。

乙炔－空气混合气体的爆炸是一种较典型的气相爆炸问题。为了便于对乙炔－空气爆炸后的流场进行数值模拟计算，必须对气相爆炸的爆炸参数进行求解。张秀华［6］给出了乙炔－空气混合气体中乙炔质量分数为7.75%时的爆炸参数计算方法，具体参数分别为：乙炔质量分数w＝7.75%，初始密度ρ0＝1.278kg／m3，比热比γ＝1.262，爆温TH＝3 560K，爆速D＝2.008km／s，爆压pH＝2.28MPa，瞬时爆压p＝1.14MPa，爆热Qv＝3 402kJ／kg，初始内能E0＝4 348kJ／m3。

蛋白质是细胞的产物，同时也是细胞中各种生命活动必不可少的物质，蛋白质分布在细胞的各个部位中。在构建模型时，学生应尽可能标记出各种高中阶段所学蛋白质及相应的生命活动，将复杂知识可视化、系统化。联系细胞的各个部位所进行的生命活动，思考会有哪些蛋白质的参与。如，细胞中的滑面型内质网上是否分布与脂质、糖合成的有关的酶、而粗面型内质网上是否分布与蛋白质加工有关的酶？细胞质基质中与细胞呼吸有关的酶（有氧呼吸第一阶段、无氧呼吸）有哪些？

在数值模拟中，空气和乙炔－空气混合气体均采用空物质材料模型和线性多项式状态方程描述［5］，具体参数分别为：空气初始密度ρ＝1.290kg／m3，C0＝－100kPa，C1＝C2＝C3＝0，C4＝C5＝0.4，C6＝0，E＝0.25MJ／m3，V0＝1.0；乙炔－空气混合气体ρ＝1.278kg／m3，C0＝C1＝C2＝C3＝0，C4＝C5＝0.262，C6＝0，E＝4.348MJ／m3，V0＝1.0。在标准状态下，乙炔密度ρ＝1.17kg／m3，初始内能E0＝58.3MJ／m3。

2.2 爆能等效原理与有限元模型

基于爆能等效原理，可求出乙炔－空气混合气体的等效TNT当量［7］，再进一步确定乙炔－空气混合气体的体积：

式中：W 为乙炔－空气混合气体的等效TNT当量，WEXP为爆炸物乙炔气体的质量，HTNT为TNT炸药的爆炸热量，HEXP为爆炸物乙炔气体爆炸的热量。

根据式（5），1m3混合气体（乙炔质量分数为7.75%）爆炸放出的热量相当于1.014kg TNT炸药爆炸放出的热量。取TNT炸药质量为1.63kg、TNT炸药计算域为边长0.1m的立方体，可以求出乙炔－空气混合气体的体积为Vh＝1.60m3、乙炔－空气混合气体计算域为边长1.15m的立方体。在进行数值模拟有限元建模时，假设空气、乙炔－空气混合气体和炸药为均匀连续介质，整个爆炸过程为绝热过程。综合考虑计算精度和计算成本，选取空气计算域为边长5m的立方体。

分别建立炸药和乙炔－空气混合气在空气中爆炸的有限元模型，为提高计算效率，基于对称性选取1／8模型进行计算。空气、炸药、乙炔－空气混合气体均采用实体单元Solid 164，采用Euler坐标描述方法。张秀华［6］对TNT炸药在自由空气中爆炸传播规律进行了数值计算，分析了网格尺寸对计算结果的影响，网格尺寸选取25mm可以满足精度要求。在对称面上施加对称边界条件，在空气计算域边界，除对称面均采用无反射边界条件。有限元模型如图1～2所示。

图1 有限元模型Fig.1 Finite element model

图2 目标测点Fig.2 Measuring target points

3 数值模拟分析

3.1 混合气体爆炸冲击波的传播特性

图3为两种爆炸源爆炸后空气冲击波传播的几个瞬时冲击波压力云图的比较。由图3（a）可见，乙炔－空气混合气体爆炸刚开始时，爆炸物并不是以球面波形式向外传播，但当超过爆炸物飞散区后，空气冲击波开始以球面波的形式向远处传播。由图3（b）可见，TNT炸药爆炸产物形成以起爆点为中心的球面波形式向外传播，正压区随着波的传播距离不断拉宽，波阵面压力随传播距离增加而迅速衰减。

由图3可以看出，以相同时间冲击波波阵面所能传播的距离衡量，在爆炸后期乙炔－空气混合气体爆炸产生的空气冲击波速度比TNT炸药爆炸产生的空气冲击波速度稍慢。

图3 不同时刻空气冲击波瞬时超压云图Fig.3 Transient over pressure contours of air blast at different time

图4 冲击波超压曲线Fig.4 Histories of air blast overpressure

图4为两种爆炸源距爆心径向距离分别为1.167和1.765m处（对应图2中点24 215、25 711）的冲击波超压曲线，冲击波超压峰值随着距爆心距离增加迅速衰减。在图4（a）中，以测点24 215为例，在t≈0.4ms时爆炸冲击波传播到测点处，在t＝0.8ms时达到峰值，峰值压力为0.381MPa，随后冲击波阵面上压力不断衰减并且出现了几次脉动小峰。而在图4（b）中，炸药爆炸空气冲击波超压达到峰值后则快速衰减，不会出现随后的脉动峰值。原因可能是由于炸药属于典型的点源爆炸，而乙炔气体爆轰是非点源爆炸，达到压力波峰值后，乙炔气体爆炸产物继续膨胀做功所致。

3.2 两种爆炸源计算结果比较

图5（a）为相同等效TNT当量的两种爆炸源、沿距爆心距离R＝1.167m（比例距离z＝0.99m／kg1／3）和R＝1.765m（比例距离z＝1.5m／kg1／3）时的冲击波超压曲线。由图可以看出，在z＝0.99m／kg1／3时混合气体爆炸的冲击波压力与化学炸药爆炸产生的爆炸冲击波超压误差较大，炸药爆炸在约0.477ms时到达峰值，峰值压力为1.019MPa，乙炔－空气混合气体爆炸在约0.797ms时到达峰值，峰值压力为0.381MPa，且达到峰值压力的时间比炸药爆炸稍有滞后。在给定空间位置上随着持续时间增加，两者峰值超压误差将变小，波形吻合逐渐趋好。随着距爆源距离的增加，两种爆炸源爆炸产生的冲击波超压峰值误差也将减小。

图5 两种爆炸源爆炸冲击波超压曲线比较Fig.5 Blast overpressure histories for two types of explosion source

3.3 乙炔－空气混合气体爆炸名义比例距离的确定

不同的爆源爆炸释放能量的方式不同，产生的压力波形和压力衰减的规律也不同，特别在离爆心较近的近区两者有很大差别，但是在远区差别不大。由于乙炔空气混合气体爆炸在爆炸初期并非为点源爆炸，随着波阵面不断向外传播，超出爆炸生成产物的极限区域后（远区），爆炸冲击波将脱离爆炸生产物开始独立存在，此时爆炸冲击波和爆炸生成产物无关，可以采用点源爆炸理论计算空气冲击波压力。对相同TNT当量的两种爆炸源，乙炔－空气混合气体初始体积比TNT炸药的初始体积大得多，因此在自由场条件下乙炔－空气混合气体爆炸的作用范围虽然能达到普通炸药爆炸达到的作用范围，但是其作用影响范围的实际比例距离比普通炸药小。乙炔－空气混合气体在爆炸时释放出的能量并没有完全转化为爆炸冲击波能量。因此对于乙炔－空气混合气体爆炸应该定义名义比例距离而不是与炸药爆炸（点源爆炸）采用相同的比例距离定义方法。

在离爆心距离较远（远区）时，基于某个比例距离下乙炔－空气混合气体爆炸产生的冲击波超压峰值与TNT炸药爆炸产生的冲击波超压峰值相等的原则，定义名义比例距离。

超压峰值p与爆炸源的TNT当量W、距爆炸中心的距离R有关。TNT炸药爆炸产生的爆炸冲击波超压峰值为：

类似地，定义乙炔－空气混合气体爆炸产生的爆炸冲击波超压峰值为：

式中：z0为乙炔－空气混合气体爆炸的名义比例距离。

根据爆炸冲击波超压峰值相等原则，即f（z）＝g（z0），可以建立z和z0的对应计算关系。选取两种爆源数值分析有限元模型对称面yOz上通过爆心对角线上的若干测点，根据爆炸冲击波超压峰值相等的原则，计算对应点源爆炸定义下的比例距离确定乙炔－空气混合气体爆炸计算的名义比例距离z0，在本文中选取若干个测点进行最小二乘拟合得到（见图6）：

式中：z为炸药爆炸定义的比例距离，z＝R／W1／3；W 为等效TNT当量，R为距爆心的实际距离；z0为乙炔－空气混合气体的名义比例距离。

由图5（b）可以看出，按上述名义比例距离计算乙炔－空气混合气体与TNT炸药爆炸冲击波超压曲线，超压峰值和波形较好吻合。

在离爆心距离较远时，z和z0为线性关系，根据比例距离的定义，从等效TNT当量的角度，也可以定义折算的乙炔气体实际质量，为简化计算定义z0＝R／W1／3ET，z＝R／W1／3，则z／z0＝W1／3ET／W1／3＝k，由此确定两种爆炸源在距爆源距离相等时产生爆炸冲击波超压峰值相等的TNT当量系数k。通过对测点数据进行最小二乘得到（见图7）：

则

式中：WET为乙炔－空气混合气体实际等效的TNT当量，W、R的单位分别为kg、m。

由式（9）可以计算，TNT当量系数k的取值范围为1.1～2.1，即为达到相同的爆炸作用参数需要增加乙炔质量，增加幅度与物理距离有关。

通过前面的分析，乙炔－空气混合气体自由空气中爆炸等效的TNT炸药的当量为：

式中：η为等效TNT炸药的修正系数，η＝1／k3。

图6 名义比例距离的拟合Fig.6 Fitting of nominal scale distance

图7 TNT当量系数的拟合Fig.7 Fitting of TNT equivalent coefficient

3.4 数值模拟与经验公式的比较

为了进一步验证数值模拟计算结果的可靠性，表1为沿爆心径向爆炸冲击波超压的数值模拟结果与前述经验公式计算结果的比较。表中，Brode、Baker和Josef Henrych经验公式分别见文献［7－9］。由表可以看出，化学炸药爆炸数值模拟结果与经验公式的计算结果在比例距离z＞1.0m／kg1／3的区域基本一致；而在爆点附近，各经验公式计算结果差别较大，LS－DYNA计算结果偏大。这与文献［10］结论一致。

由表1可以看出，除个别点外，采用名义比例距离拟合公式计算的超压峰值与炸药爆炸的超压峰值均较好吻合，表明名义比例距离公式（8）是合理的。在进行可燃气体爆炸分析时，利用公式（10）或（11）均可以进一步计算乙炔－空气混合气体实际等效TNT当量。

表1 TNT炸药和乙炔气体爆炸峰值压力Table 1 Peak overpressure of TNT dynamite and mixture explosive gas

分析结果表明，在比例距离z＞1.0m／kg1／3（远区）时，整个流场内冲击波峰值压力随距离的衰减规律和经验公式较好吻合，说明基于LS－DYNA数值模拟化学炸药和可燃气体爆炸冲击波在空气中的传播过程是可行的。随着距爆炸源距离的增加，两种爆源的计算结果误差逐渐趋小，两种爆炸源爆炸产生的峰值压力将逐渐接近。由于乙炔爆炸非点源爆炸，把乙炔－空气混合气体实际作用距离按名义比例距离进行修正计算，最大误差为11.6%，误差在合理的范围内，计算结果与化学炸药爆炸产生冲击波超压计算结果较一致。

4 结 论

（1）基于LS－DYNA软件模拟化学炸药和可燃气体爆炸冲击波在空气中的传播过程是可行的。提出了基于爆能等效原则的可燃气体爆炸冲击加载的数值模拟计算方法，通过控制乙炔－空气混合气体的等效TNT当量，可以控制空气冲击波超压和压力作用时间等。

（2）相同装药形状和相同等效TNT当量的两种爆炸源，随着距爆源距离的增加，两种爆炸压力波形间的误差逐渐减小，进一步分析结果表明需要加大乙炔气体等效质量提高冲击波传播速度和峰值压力；另外，采用名义比例距离公式修正后，乙炔－空气混合气体爆炸与炸药爆炸冲击波计算误差可以得到有效控制。当冲击波超压p＜0.5MPa时，采用乙炔－空气混合气体爆炸加载方式可以替代化学炸药用于结构（构件）爆炸冲击实验加载。

本文工作可为大型模爆器的结构（构件）可燃气体爆炸冲击实验加载控制以及冲击波超压计算提供参考。

［1］欧进萍，张春巍，张晓漪，等.大型模爆器综合试验系统研制、改造与应用［C］∥中国土木工程学会防护工程分会第十次学术年会.乌鲁木齐，2006.

［2］张春巍.结构抗爆抗冲击与振动控制的若干问题研究［R］.哈尔滨：哈尔滨工业大学，2007.

［3］李生娟，毕明树，章正军，等.气体爆炸研究现状及发展趋势［J］.化工装备技术，2002，23（6）：15－19.Li Sheng－juan，Bi Ming－shu，Zhang Zheng－jun，et al.Gas explosion present research situation and development trend［J］.Chemical Equipment Technology，2002，23（6）：15－19.

［4］Benson D J.Computational methods in Lagrangian and Eulerian hydrocodes［J］.Computer Methods in Applied Mechanics and Engineering，1992，99：235－394.

［5］LS－DYNA keyword user’s manual：Version 971［M］.Livermore California：Livermore Software Technology Corporation，2006.

［6］张秀华.气爆炸冲击作用下钢框架抗爆性能试验研究与数值模拟［D］.哈尔滨：哈尔滨工业大学，2011.

［7］Smith P D，Hetherington J G.Blast and ballistic loading structures［M］.Butterworth－Heinemann Ltd，1994：30－166.

［8］Mays G C，Smith P D.Blast effects on buildings［M］.London：Thomas Telford Publications，1995：24－50.

［9］Baker W E.Explosions in air［M］.Austin，TX：University of Texas Press，1973：7－15.

［10］王飞，朱立新，顾文彬，等.基于 ALE算法的空气冲击波绕流数值模拟研究［J］.工程爆破，2002，8（2）：13－16.Wang Fei，Zhu Li－xin，Gu Wen－bin，et al.Numerical simlation shock wave around－flow on basis of ALE algorithm［J］.Engineering Blasting，2002，8（2）：13－16.