2024年4月5日发(作者：何凝珍)
2018
年第
3
期
导
弹
与
航
天
运
载
技
术
No.3 2018
总第
361
期
MISSILES AND SPACE VEHICLES
Sum No.361
文章编号：1004-7182(2018)03-0098-05 DOI：10.7654/.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、
收稿日期：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=P
0
，
t
<
t
a
（1）
⎛
t
−
t
a
P
=
P
0
+
P
i
⎜
1
−
t
d
⎝
⎞
⎟
e
⎠
−
b(t
−
t
a
)
t
d
，
t
a
≤
t
≤
t
a
+
t
d
（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
)
)
23
Fig.1 Schematic Profile of Airblast Shock Wave in Free Air
由图1可知，在冲击波到达之前，该处的压力等于
大气压力
P
0
，冲击波在时间
t
a
到达该处后，压力瞬间由
大气压力突跃至最大值。压力最大值与
P
0
的差值，通常
称为入射超压峰值
P
i
。波阵面通过后压力迅速下降，经
过时间
t
d
压力衰减到大气压力并继续下降，直至出现负
超压峰值，在一定时间内又逐渐地回升到大气压力。负
超压与
P
0
的差值，通常称为负超压峰值
P
n
。
空气冲击波压力在正压段大致按指数规律衰减，
有许多经验公式可以描述此压力衰减过程，其中修正
式中
FUNCTION
代表冲击波到达时间
t
a
、入射超压
P
i
及正压作用时间
t
d
等参数；
Z
为比例距离，
Z
=
d
/
3
W
；
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
冲击波到达时间
t
a
/(
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
入射超压
P
i
/
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
正压作用时间
t
d
/(
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.2777
Z
−
1.1975
（4）
负压区冲击波参数虽然对刚性结构（钢筋混凝土）
的设计不重要，但对柔性防护结构（一般指钢框架）
尤其重要。Martin Larcher
[9]
采用双折线方程来计算负
100
导
弹
与
航
天
运
载
技
术
2018
年
压区压力：
t
n
2
P
n
⎧
⎪
P
0
−
t
(
t
−
t
a
−
t
d
)
，
t
a
+
t
d
<
t
≤
t
a
+
t
d
+
2
n
⎪
t
2
P
n
⎪
P
=
⎨
P
0
−
(
t
a
+
t
d
+
t
n
−
t
)
，
t
a
+
t
d
+
n
<
t
≤
t
a
+
t
d
+
t
n
（5）
t
n
2
⎪
Ptttt
，
>++
⎪
0adn
⎪
⎩
只有一个二维四边形单元。而LS-DYNA可采用一维
梁单元球对称计算模型，如图2所示。在计算模型中
采用多物质Euler算法和缺省的人工粘性系数1×10
-6
，
时间步长缩放因子均为0.9。
式中
P
0
为大气压力；
P
n
为负超压；
t
n
为负超压持续
时间。
对于负压区参数（负超压及其持续时间），Martin
Larcher
对Krauthammer
出下列计算公式。
负超压计算公式：
⎧
0.35
5
10
，
Z
≥
3.5
⎪
P
n
=
⎨
Z
（6）
4
⎪
10
，
Z
<
3.5
⎩
[9][14]
提供的图表进行拟合，得
a
）
AUTODYN
二维楔形网格轴对称计算模型
b
）
LS-DYNA
一维梁单元球对称计算模型
图
2 AUTODYN
和
LS-DYNA
计算模型
负超压持续时间计算公式：
⎧
0.0104
W
13
，
Z
<
0.3
⎪
t
n
=
⎨
(
0.003125log(
Z
)
+
0.01201
)
W
13
，
0.3
≤
Z
≤
1.9
（
7
）
⎪
0.0139
W
13
，
Z
>
1.9
⎩
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状态方程来描述。
⎛⎛
ω
⎞
−
R
1
V
ω
⎞
−
R
2
V
ω
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
，
R
1
，
R
2
，
ω
为常数。其中，方程式右端第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×10
8
3.231×10
6
4.15 0.95 0.3 1630 6930 7.0×10
9
2.1×10
7
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 (2
nd
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.
2024年4月5日发(作者：何凝珍)
2018
年第
3
期
导
弹
与
航
天
运
载
技
术
No.3 2018
总第
361
期
MISSILES AND SPACE VEHICLES
Sum No.361
文章编号：1004-7182(2018)03-0098-05 DOI：10.7654/.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、
收稿日期：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=P
0
，
t
<
t
a
（1）
⎛
t
−
t
a
P
=
P
0
+
P
i
⎜
1
−
t
d
⎝
⎞
⎟
e
⎠
−
b(t
−
t
a
)
t
d
，
t
a
≤
t
≤
t
a
+
t
d
（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
)
)
23
Fig.1 Schematic Profile of Airblast Shock Wave in Free Air
由图1可知，在冲击波到达之前，该处的压力等于
大气压力
P
0
，冲击波在时间
t
a
到达该处后，压力瞬间由
大气压力突跃至最大值。压力最大值与
P
0
的差值，通常
称为入射超压峰值
P
i
。波阵面通过后压力迅速下降，经
过时间
t
d
压力衰减到大气压力并继续下降，直至出现负
超压峰值，在一定时间内又逐渐地回升到大气压力。负
超压与
P
0
的差值，通常称为负超压峰值
P
n
。
空气冲击波压力在正压段大致按指数规律衰减，
有许多经验公式可以描述此压力衰减过程，其中修正
式中
FUNCTION
代表冲击波到达时间
t
a
、入射超压
P
i
及正压作用时间
t
d
等参数；
Z
为比例距离，
Z
=
d
/
3
W
；
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
冲击波到达时间
t
a
/(
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
入射超压
P
i
/
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
正压作用时间
t
d
/(
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.2777
Z
−
1.1975
（4）
负压区冲击波参数虽然对刚性结构（钢筋混凝土）
的设计不重要，但对柔性防护结构（一般指钢框架）
尤其重要。Martin Larcher
[9]
采用双折线方程来计算负
100
导
弹
与
航
天
运
载
技
术
2018
年
压区压力：
t
n
2
P
n
⎧
⎪
P
0
−
t
(
t
−
t
a
−
t
d
)
，
t
a
+
t
d
<
t
≤
t
a
+
t
d
+
2
n
⎪
t
2
P
n
⎪
P
=
⎨
P
0
−
(
t
a
+
t
d
+
t
n
−
t
)
，
t
a
+
t
d
+
n
<
t
≤
t
a
+
t
d
+
t
n
（5）
t
n
2
⎪
Ptttt
，
>++
⎪
0adn
⎪
⎩
只有一个二维四边形单元。而LS-DYNA可采用一维
梁单元球对称计算模型，如图2所示。在计算模型中
采用多物质Euler算法和缺省的人工粘性系数1×10
-6
，
时间步长缩放因子均为0.9。
式中
P
0
为大气压力；
P
n
为负超压；
t
n
为负超压持续
时间。
对于负压区参数（负超压及其持续时间），Martin
Larcher
对Krauthammer
出下列计算公式。
负超压计算公式：
⎧
0.35
5
10
，
Z
≥
3.5
⎪
P
n
=
⎨
Z
（6）
4
⎪
10
，
Z
<
3.5
⎩
[9][14]
提供的图表进行拟合，得
a
）
AUTODYN
二维楔形网格轴对称计算模型
b
）
LS-DYNA
一维梁单元球对称计算模型
图
2 AUTODYN
和
LS-DYNA
计算模型
负超压持续时间计算公式：
⎧
0.0104
W
13
，
Z
<
0.3
⎪
t
n
=
⎨
(
0.003125log(
Z
)
+
0.01201
)
W
13
，
0.3
≤
Z
≤
1.9
（
7
）
⎪
0.0139
W
13
，
Z
>
1.9
⎩
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状态方程来描述。
⎛⎛
ω
⎞
−
R
1
V
ω
⎞
−
R
2
V
ω
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
，
R
1
，
R
2
，
ω
为常数。其中，方程式右端第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×10
8
3.231×10
6
4.15 0.95 0.3 1630 6930 7.0×10
9
2.1×10
7
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 (2
nd
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.