TNT空中爆炸冲击波的工程和数值计算

文章编号：1004-7182(2018)03-0098-05 DOI：10.7654/j.issn.1004-7182.20180319

TNT空中爆炸冲击波的工程和数值计算

辛春亮1，王俊林1，余道建1，李 通2，宋师军1

（1. 北京航天长征飞行器研究所，北京，100076；2. 中国运载火箭技术研究院，北京，100076）

摘要：针对TNT空气中爆炸冲击波压力，建立了空气自由场爆炸冲击波工程计算模型，并采用AUTODYN和LS-DYNA软件中的一维计算模型对TNT空气自由场爆炸冲击波压力衰减过程进行了数值模拟研究。为了捕捉峰值，数值计算模型中采用细密的网格和很小的人工粘性系数。数值模拟结果与工程计算结果大致吻合，但AUTODYN计算耗时远高于LS-DYNA。

关键词：冲击波；工程计算模型；数值模拟；TNT 中图分类号：TJ55；O383 文献标识码：A

Empirical Formula and Numerical Simulation of TNT Explosion Shock Wave

in Free Air

Xin Chun-liang1, Wang Jun-lin1, Yu Dao-jian1, Li Tong2, Song Shi-jun1

(1. Beijing Institute of Space Long March Vehicle, Beijing, 100076; 2. China Academy of Launch Vehicle Technology, Beijing, 100076)

Abstract: Empirical formula of blast waves in free air is presented in the paper. Calculation and numerical simulation of TNT

explosion in air are analyzed using these empirical formula, AUTODYN and LS-DYNA software respectively. Fine meshes and small hourglass coefficient are used in one-dimensional simulation model to capture peak pressure of shock wave. Simulated pressure profiles agree approximately with these calculated by empirical formula. The computational CPU time by AUTODYN is much longer than that by LS-DYNA.

Key words: Shock wave; Empirical formula; Numerical simulation; TNT

0 引 言

炸药在空气中爆炸后瞬间形成高温高压的爆炸产物，产物强烈压缩周围静止的空气，在空气中形成冲击波向四周传播，对结构造成破坏。许多学者对于TNT空中爆炸冲击波超压和传播规律进行了研究[1~5]，总结出了Brode、Henrych、Kingery-Bulmash、Kinney-Graham等超压经验和半经验公式，并编写了CONWEP、BEC、ATBLAST、BLAPAN、3D-BLAST、EBLAST、SHOCK、BEAM、BLASTX等工程计算程序，这些超压计算公式和工程计算程序是在总结了大量试验数据的基础上发展起来的，大都将炸药按爆热折算成等效TNT当量，不考虑负压区和冲击波的多次反射以及稀疏作用，并假设压力以指数规律衰减，因此主要适用于自由场环境下比例距离较大时超压的快速计算。数值模拟可以解决超压计算公式和工程计算程序难以解决的更为复杂的问题，这方面的研究成果也很多[6~9]，大都采用基于有限差分、有限体积法（如AUTODYN、Air3D、

收稿日期：2017-04-25；修回日期：2017-07-11

CTH、DYTRAN、DYSMAS、IFSAS、SHAMRC）以及有限元法软件（如LS-DYNA、EUROPLEXUS）中的显式积分算法，数值计算超压与试验测试结果的吻合程度也各不相同，大体来说，基于有限差分和有限体积法的软件计算准确度较高，有限元软件则差一些，其计算超压曲线的显著特点是：a）峰值之前有明显的压力爬升过程；b）超压峰值容易被抹平，峰值过后压力衰减过快；c）网格尺寸越小，计算结果越准确。

本文将空气冲击波正压区和负压区工程计算公式及其参数结合起来，形成TNT空气自由场爆炸冲击波工程计算模型，并采用有限差分软件AUTODYN和有限元软件LS-DYNA对无限空间TNT空气中爆炸冲击波传播过程进行数值模拟，最后将三者得出的计算结果进行对比分析。

1 空气冲击波工程计算模型

典型的TNT空气中爆炸冲击波曲线如图1所示。

作者简介：辛春亮（1973-），男，博士，研究员，主要研究方向为战斗部设计和数值模拟

第3期

辛春亮等 TNT空中爆炸冲击波的工程和数值计算

99

的Friedlander方程[2]较接近实际过程且又简单易于计算：

P=P0，t⎛t−ta

P=P0+Pi⎜1−

td⎝

⎞⎟e⎠

−b(t−ta)td

，ta≤t≤ta+td （2）

1984年，Kingery在对大量试验数据总结分析的基础上提出了Kingery-Bulmash空气冲击波参数计算模型[4,10,11]，其研究成果已被广泛采用，并被植入多种计

图1 自由场空气冲击波压力-时间曲线

算程序如CONWEP和LS-DYNA（即*LOAD_BLAST关键字）中。经过一系列的试验测试，1998年Kingery远场对Kingery-Bulmash超压模型作了一些修改[12,13]，

超压预测值比以前高出许多，此修正模型已在BEC工程计算程序中实现。

Kingery冲击波正压区参数计算公式：

FUNCTION=EXP(A+BlnZ+C(lnZ)+D(lnZ)+

（3） 456

E(lnZ)+F(lnZ)+G(lnZ))

2

3

Fig.1 Schematic Profile of Airblast Shock Wave in Free Air

由图1可知，在冲击波到达之前，该处的压力等于大气压力P0，冲击波在时间ta到达该处后，压力瞬间由大气压力突跃至最大值。压力最大值与P0的差值，通常称为入射超压峰值Pi。波阵面通过后压力迅速下降，经过时间td压力衰减到大气压力并继续下降，直至出现负超压峰值，在一定时间内又逐渐地回升到大气压力。负超压与P0的差值，通常称为负超压峰值Pn。

空气冲击波压力在正压段大致按指数规律衰减，有许多经验公式可以描述此压力衰减过程，其中修正

式中 FUNCTION代表冲击波到达时间ta、入射超压Pi及正压作用时间td等参数；Z为比例距离，Z=d/3W；

d为距离爆心的距离；W为装药质量。当比例距离 Z>40 m·kg

-1/3

时，对于大多数结构已无破坏能力。A，

B，C，D，E，F，G是系数，具体取值见表1。

表1 简化的Kingery空气冲击波参数

Tab.1 Shock Wave Parameters of Simplified Kingery Formula

冲击波到达时间ta/(ms⋅kg-1/3)

RANGE，Z (m⋅kg-1/3) 0.06~1.50 1.50~40

A B C D E F G -0.7604 1.8058 0.1257 -0.0437

-0.0310

-0.00669 0 -0.7137 1.5732 0.5561 -0.4213 0.1054 -0.00929 0 入射超压Pi/kPa

RANGE，Z (m⋅kg-1/3)

A B C D E F G

-0.3229 0.1117 0.0685 0 0 0.2~2.9 7.2106 -2.1069

2.9~23.8 7.5938 -3.0523 0.40977 0.0261 -0.01267 0 0 23.8~198.5 6.0536 -1.4066 0 正压作用时间td/(ms⋅kg-1/3)

RANGE，Z (m⋅kg-1/3)

0 0 0 0

A B C D E F G

-5.9667

-4.0815

-0.9149 0 0.2~1.02 0.5426 3.2299 -1.5931

1.02~2.80 0.5440 2.7082 -9.7354 14.3425 -9.7791 2.8535 0 2.80~40

-2.4608 7.1639 -5.6215 2.2711 -0.44994 0.03486 0

Martin Larcher[9]提出了下面的计对于衰减系数b，

算公式：

b=5.2777Z−1.1975 （4）

负压区冲击波参数虽然对刚性结构（钢筋混凝土）的设计不重要，但对柔性防护结构（一般指钢框架）尤其重要。Martin Larcher[9]采用双折线方程来计算负

100

导 弹 与 航 天 运 载 技 术 2018年

压区压力：

tn2Pn⎧

⎪P0−t(t−ta−td)，ta+tdn⎪t2Pn⎪

P=⎨P0−(ta+td+tn−t)，ta+td+ntn2⎪Ptttt，>++⎪0adn⎪⎩

只有一个二维四边形单元。而LS-DYNA可采用一维

梁单元球对称计算模型，如图2所示。在计算模型中采用多物质Euler算法和缺省的人工粘性系数1×10-6，时间步长缩放因子均为0.9。

式中 P0为大气压力；Pn为负超压；tn为负超压持续时间。

对于负压区参数（负超压及其持续时间），Martin Larcher对Krauthammer出下列计算公式。

负超压计算公式：

⎧0.355

10，Z≥3.5⎪

Pn=⎨Z （6）

4⎪10，Z<3.5⎩

[9]

[14]

提供的图表进行拟合，得

a）AUTODYN二维楔形网格轴对称计算模型

b）LS-DYNA一维梁单元球对称计算模型 图2 AUTODYN和LS-DYNA计算模型

负超压持续时间计算公式：

⎧0.0104W13，Z<0.3⎪

tn=⎨(0.003125log(Z)+0.01201)W13，0.3≤Z≤1.9（7） ⎪0.0139W13，Z>1.9⎩

Fig.2 AUTODYN and LS-DYNA Computational Models

计算模型中的炸药为1 kg TNT，计算空气域为 12 m，网格划分总数为6000，距离爆心1 m内的网格尺寸为1 mm，然后由小到大逐渐向外围渐变。 在上述计算模型中，空气密度为1.225 kg/m3，采用理想气体状态方程，γ=1.4。空气中初始压力为1个标准大气压（0.1 MPa）。

TNT炸药爆炸产物压力用JWL状态方程来描述。

⎛⎛ω⎞−R1Vω⎞−R2VωE

+B⎜1−+ （8） P=A⎜1−⎟e⎟e

⎜RV⎟⎜RV⎟V⎝⎠⎝⎠12

将上述空气冲击波正压区和负压区工程计算公式

结合起来，便可形成TNT空气自由场爆炸冲击波工程计算模型，该模型可以快速计算出不同TNT当量在不同距离处的冲击波压力曲线。

2 数值计算模型

由于炸药爆炸初期产生的冲击波是高频波，在数值计算模型中炸药及其附近区域需要划分细密网格才能反映出足够频宽的冲击波特性，否则计算出的压力峰值会被抹平。Martin Larcher[9]建议，对于1 kg TNT空气中爆炸数值模拟，炸药及其附近空气的网格尺寸在1 mm左右，如此细密的网格划分方式会导致网格数量极其庞大。无限空间TNT空中爆炸问题具有球对称性质，一维计算模型是最佳选择，可显著降低计算规模和计算时间。

AUTODYN软件中的二维楔形网格轴对称计算模型可用来等效一维计算模型，垂直于冲击波传播方向

式中 E为单位质量内能；V为比容；A，B，R1，R2，ω为常数。其中，方程式右端第1项在高压段起主要作用，第2项在中压段起主要作用，第3项代表低压段。在爆轰产物膨胀的后期，方程式前两项的作用可以忽略，为了加快求解速度，将炸药从JWL状态方程转换为更为简单的理想气体状态方程（绝热指数γ=ω+1）。TNT炸药爆炸产物JWL状态方程参数取值来自文献[15]。TNT爆炸产物JWL状态方程参数如表2所示。

表2 TNT爆炸产物JWL状态方程参数

Tab.2 JWL EOS Parameters of TNT Explosion Products

A/kPa

B/kPa

R1

R2

ω

Ρ/(kg·m-3)

D/(m·s-1)

E/(J·m-3)

PCJ/kPa

3.712×108 3.231×106 4.15 0.95 0.3 1630 6930 7.0×109 2.1×107

3 工程计算和数值计算结果的对比

图3是距离爆心不同距离处冲击波压力计算曲线。 图3中AUTODYN和LS-DYNA数值计算曲线上均存

多个峰值，第2个峰值是由TNT炸药的爆心汇聚反射

追赶造成的。由于空气密度和压强远小于爆轰产物的密度和压强，在冲击波形成的同时由界面向炸药中心反射

第3期

辛春亮等 TNT空中爆炸冲击波的工程和数值计算

101

回一个稀疏波，使爆轰产物发生膨胀，降低内部压力，由图3可知，AUTODYN和LS-DYNA数值计算此稀疏波在炸药中心汇聚后又向外传播一个压缩波，由于前导冲击波已经将空气绝热压缩，此压缩波的传播速度将大于前导冲击波，并逐渐向前追赶前导冲击波。如果测点离炸药不远，此二次压力波峰值足够高，就有可能将负压区强行打断，并再次衰减到负压。

a）距离爆心1 m

b）距离爆心2 m

c）距离爆心4 m

d）距离爆心6 m

图3 不同距离处冲击波压力-时间计算曲线

Fig.3 Comparition of Airblast Shock Wave Profiles of Different Ranges

曲线极为相似，走势与工程计算曲线基本一致。距爆心越远，3条曲线符合程度越高，且数值计算曲线提供的信息更加丰富。

数值计算曲线与工程计算曲线的区别主要在于： a）冲击波到达时间：工程计算曲线冲击波峰值最高，因此最早到达，其次是AUTODYN，LS-DYNA计算曲线冲击波到达时间最晚。由于数值计算模型中采用了细密网格和很小的人工粘性系数，计算压力时间曲线爬升段接近强间断。

b）正压峰值：工程计算值最高，其次是AUTODYN计算值，LS-DYNA计算值最低，且距离爆心越近，也就是比例距离越小，差别越明显。

c）正压作用时间：工程计算值略高于数值计算值。

d）正压区冲量：由b）和c）基本可以得出，工程计算值高于数值计算值，比例距离越小，差别越明显。

e）负压峰值：工程计算值与数值计算峰值基本相当。负压峰值到达时间提前。

f）负压峰值到达时间：AUTODYN计算曲线负压峰值到达时间最早，其次是LS-DYNA，工程计算曲线负压峰值最晚到达。

g）负压作用时间：数值计算峰值与工程计算值基本相当。

对于不同的比例距离，上述规律均相同。此外，由于工程计算模型简单，考虑因素少，耗时短，瞬间即可完成计算。LS-DYNA计算耗时4 min 52 s。AUTODYN软件本身计算速度比LS-DYNA慢，再加上二维四边形单元比一维梁单元计算耗费大，因此AUTODYN计算耗时最长，为34 min，是LS-DYNA的7倍。

4 结 论

a）本文将把Kingery、Bulmash、Martin Larcher以及Krauthammer等人在空气冲击波正压区和负压区参数计算研究成果结合起来，形成TNT空气自由场爆炸冲击波工程计算模型，通过该模型可以快速计算出TNT空气自由场爆炸冲击波压力衰减演化曲线。该模型的正压区参数计算公式由大量实验数据拟合而成，

准确可靠。由于过于简单地采用双折线模型将负压区持续时间均分，与实际会存在较大偏差。

b）AUTODYN数值计算曲线与工程计算曲线吻合程度最高，负压区包含的信息比工程计算曲线更为丰富，也更接近实际。

102

导 弹 与 航 天 运 载 技 术 2018年

[8] Alia A, Souli M. High explosive simulation using multi-material

formulations[J]. Applied Thermal Engineering, 2006, 26(10):1032-1042. [9] Martin L. Simulation of the effects of an air blast wave[R]. Ispra:

c）LS-DYNA计算曲线与工程计算曲线吻合程度稍差，但LS-DYNA计算时间远低于AUTODYN。

参 考 文 献

[1] Brode H L. Numerical solutions of spherical blast wave[J]. JAP, 1955,

26(6): 766-775.

[2] Baker W E. Explosions in air[M]. Austin: University of Texas Press, 1973. [3] Henrych J. The dynamics of explosion and its use[M]. Amsterdam:

Elsevier, 1979.

[4] Kingery, Charles N, Bulmash, Gerald. Airblast parameters from TNT

spherical air burst and hemispherical surface burst[R]. Maryland: Defence Technical Information Center, Ballistic Research Laboratory, Aberdeen Proving Ground, 1984.

[5] Kinney G F, Graham K J. Explosive shocks in air (2nd Edition)[M]. New

York: Springer-Verlag, 1985.

[6] Borgers J B. Vantomme J. Towards a parametric model of a planar

blastwave created with detonating cord[C]. Ralston: 19th military aspects of blast and shock, 2006.

[7] Clutter J K. Stahl M. Hydrocode simulations of air and water shocks for

facility vulnerability assessments[J]. Journal of Hazardous Materials, 2004(106): 9-24.

JRC41337, 2007.

[10] Army Engineer Waterways Experiment Station. Fundamentals of

protective design for conventional weapons[R]. Vicksburg: US Department of the Army, WES/IR/SL-1, 1987.

[11] Conwep Software, Hyde D W. Conventional weapon effects[R].

Albuquerque: US Army Engineering Waterways Experimental Station, 919167, 1992.

[12] TB 700-2, NAVSEAINST 8020.8 B. DoD ammunition and explosives[R].

Washington: Hazards Classification Procedures, DOD-4145.26-M-REV, 1998.

[13] Swisdak M M. Simplified kingery airblast calculations[C]. Florida:

Proceedings of the 26th DoD Explosives Safety Seminar, 1994.

[14] Krauthammer T, Altenberg A. Negative phase blast effects on glass

panels[J]. International Journal of Impact Engineering, 2000, 24(1): 1-18. [15] Dobratz B M, Crawford P C. LLNL explosives handbook - properties of

chemical explosives and explosive stimulants[M]. California: Lawrence Livermore National Laboratory, UCRL-52997, 1985.

（上接第97页）

参 考 文 献

[1] 顾乃威, 王丽伟, 等. 地面设备伪装隐身评估方法研究[J]. 导弹与航天

运载技术, 2016(6): 86.

Gu Naiwei, Wang Liwei, etc al. Study on camouflage and stealthy capability evaluation method of ground support equipments[J]. Missiles And Space Vehicles, 2016(6): 86.

[2] 贾鑫, 叶伟, 吴彦宏, 王宏艳. 合成孔径雷达对抗技术[M]. 北京: 国防

工业出版社, 2004.

Jia Xin, Ye Wei, Wu Yanhong, Wang Hongyan. Electronic countermeasure technology to synthetic aperture radar[M]. Beijing: National Defense Industry Press, 2004.

[3] 萨布林, 瓦切斯拉夫, 尼卡拉伊维奇. 侦察-打击一体化系统和对地观

测雷达系统[M]. 吴飞, 译. 北京: 国防工业出版社, 2005.

Sa Bulin, Wa Qiesilafu, Ni Kalayiweiqi. Reconnaissance and strike integration system and ground observation radar system[M]. Wu Fei, Translated. Beijing: National Defense Industry Press. 2005.

[4] Skolnik M I, 主编. 雷达手册[M]. 王军, 林强, 等, 译. 北京: 电子工业

出版社, 2007.

Skolnik M I. Radar handbook[M]. Wang Jun, Lin Qiang, et al. Beijing: Publishing House of Electronics Industry, 2007.

[5] 张考, 马东立. 军用飞机生存力与隐身设计[M]. 北京: 国防工业出版

社, 2002.

Zhang Kao, Ma Dongli. Military aircraft survivability and stealth design[M]. Beijing: National Defense Industry Press, 2002.

[6] 赵渊, 彭海军. 军事伪装与战争欺骗[M]. 北京: 化学工业出版社, 2013.

Zhao Yuan, Peng Haijun. Military disguise and war deception[M]. Beijing: Chemical Industry Press, 2013.

[7] 刘冶, 李竹影, 等. 组合型电磁隐身斗篷的超材料设计与仿真[J]. 功能

材料, 2013, 15(44): 2235-2238.

Liu Ye, Li Zhuying, et al. Design and emulation of combined-shaped electromagnetic stealthy cloak made of metamaterials[J]. Journal of Functional Materials, 2013, 15(44): 2235-2238.

[8] 杨长胜, 程海峰, 等. 智能隐身材料的研究现状[J]. 功能材料, 2005,

36(5): 643-647.

Yang Changsheng, Cheng Haifeng, et al. Present status of intelligent stealth material[J]. Journal of Functional Materials, 2005, 36(5): 643-647.