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

2018年第3期导弹与航天运载技术 No.3 2018 总第361期 MISSILES AND SPACE VEHICLES Sum No.361

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

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

文章编号：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-liang 1, Wang Jun-lin 1, Yu Dao-jian 1, Li Tong 2, Song Shi-jun 1

(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 、

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

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

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

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

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

第3期

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

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

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

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

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

0P P =，a t t < （1）

a d

()a

0i d 1e b t t t t t P P P t ?????=+?

???

，a a d t t t t ≤≤+ （2）

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

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

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

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

()()()()()2

456

FUNCTION EXP(ln ln ln ln ln ln )

A B Z C Z D Z E Z F Z G Z =++++

++（3）式中 FUNCTION 代表冲击波到达时间a t 、入射超压i P

及正压作用时间d t 等参数；Z 为比例距离，Z d =；

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

冲击波到达时间a t /(-1/3ms kg ?)

RANGE ，Z (-1/3m kg ?) A B C D E F G 0.06~1.50 -0.7604 1.8058 0.1257 -0.0437

-0.0310

-0.00669 0

1.50~40

-0.7137 1.5732 0.5561 -0.4213 0.1054 -0.00929 0

入射超压i P /kPa

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

A B C D E F G

0.2~2.9 7.2106 -2.1069

-0.3229 0.1117 0.0685

0 0

2.9~2

3.8 7.5938 -3.0523 0.40977 0.0261 -0.01267 0 0 23.8~198.5 6.0536 -1.4066 0

0 0

正压作用时间d t /(-1/3ms kg ?)

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

A B C D E F G

0.2~1.02 0.5426 3.2299 -1.5931

-5.9667

-4.0815

-0.9149 0

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

对于衰减系数b ，

Martin Larcher [9]提出了下面的计算公式：

1.19755.2777b Z ?= （4）

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

导弹与航天运载技术 2018年

100

压区压力：

()()n n 0a d a d a d n n n 0a d n a d a d n n 0a d n 22

22P t P t t t t t t t t t P t P P t t t t t t t t t t t P t t t t ?

???+<≤++???

=?++?++<≤++??>++???

，，，（5）

式中 0P 为大气压力；n P 为负超压；n t 为负超压持续时间。

对于负压区参数（负超压及其持续时间），Martin Larcher [9]

对Krauthammer [14]

提供的图表进行拟合，得

出下列计算公式。

负超压计算公式：

n 40.3510 3.510 3.5Z P Z Z ?≥?

=??

，，（6）

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

()1313n 130.01040.30.003125log()0.012010.3 1.90.0139 1.9W Z t Z W Z W Z ?

=+≤≤??>?

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

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

2 数值计算模型

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

称性质，一维计算模型是最佳选择，可显著降低计算规模和计算时间。

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

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

可采用一维

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

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

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

Fig.2 AUTODYN and LS-DYNA Computational Models

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

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

1212

1e 1e

R V R V E

P A B R V R V V ωωω??????=?+?+???????????? （8）式中 E 为单位质量内能；V 为比容；A ，B ，1R ，2R ，ω为常数。其中，方程式右端第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

R 1

R 2

Ρ/(kg·m -3)

D /(m·s -1)

E /(J·m -3)

P CJ /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 炸药的爆心汇聚反射

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

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

101

第3期

回一个稀疏波，使爆轰产物发生膨胀，降低内部压力，此稀疏波在炸药中心汇聚后又向外传播一个压缩波，由于前导冲击波已经将空气绝热压缩，此压缩波的传播速度将大于前导冲击波，并逐渐向前追赶前导冲击波。

如果测点离炸药不远，此二次压力波峰值足够高，就有可能将负压区强行打断，并再次衰减到负压。

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可知，AUTODYN 和LS-DYNA 数值计算曲线极为相似，走势与工程计算曲线基本一致。距爆心越远，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 数值计算曲线与工程计算曲线吻合程度最高，负压区包含的信息比工程计算曲线更为丰富，也更接近实际。

导弹与航天运载技术 2018年102

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. [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:

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.