从一般的火炮到集各种先进技术于一身的导弹，大多数现代武器的战斗部速度均大于音速。战斗部在高速运动中爆炸所带来的各种特殊效应与战斗部静爆时相比有着特殊之处^[1]。因此，研究战斗部动爆条件下的爆炸冲击波特性，对于战斗部的威力评价、武器毁伤效能评估以及动爆试验测试方案设计均有重要意义^[2]。

冲击波静爆超压计算公式大多是建立在爆炸相似律的基础上，根据实测的多组数据拟合出来的经验公式，进行修正得到^[3]。比较著名的经验公式有Henrych公式、萨多夫斯基公式、Brode公式等^[4]。在动爆超压计算模型方面，现有的动爆理论中对于运动装药的处理，一般采取等效药量的方法。从能量相似原理出发，可以将运动TNT装药导致的威力增加与静止TNT装药量的增加相等效。根据等效药量公式结合上述静爆超压经验公式，可得到相应的动爆超压计算公式^[5]；杜红棉等^[1]得到了运动装药在空中爆炸时形成的冲击波正方向和反方向的超压计算公式。蒋海燕等^[6]利用AUTODYN软件对装药动爆冲击波场进行仿真，分析了动爆冲击波场的分布规律，并对数据进行回归分析，建立了与动爆试验结果符合较好的工程计算模型；聂源等^[7]采用高精度显示欧拉流体力学软件SPEED，模拟了球形装药在空气介质中的爆炸过程，得到了相似的结果。

本文中，通过AUTODYN软件进行建模计算，研究装药速度对爆炸冲击波流场演化的影响。分析研究动爆冲击波的冲击波超压、正压作用时间、比冲量的特性。并且从物理原理出发，利用模拟数据进行拟合分析，研究动爆冲击波对角度的依赖关系。最后建立修正因子δ的表达式，据此出发可以得到动爆冲击波的超压计算公式，与已有的结果相比应具备更好的准确度。

1.   数值模拟

1.1   建立数值计算模型

为准确分析爆炸瞬时装药运动速度对冲击波威力场的影响规律，选用TNT球形裸装药进行计算。采用中心起爆，装药质量W=0.170 1 kg，装药密度ρ=1.63 g/cm³。采用二维轴对称模型，计算空气域尺寸为2 400 mm×1 200 mm，网格尺寸为2 mm×2 mm。采用多物质欧拉算法，材料直接从AUTODYN材料库中选取，空气采用理想气体状态方程，炸药采用JWL状态方程。在空气域的x轴中间填充质量为0.170 1 kg的球形TNT炸药。边界条件设定为压力流出边界条件，以模拟无限空气域。如图1所示，将TNT炸药中心定义为坐标原点，各监测点成环形围绕着坐标原点排列于空气域中，距爆心300 mm起间隔200 mm放置4个监测点，每列监测点之间夹角为30°，共7列28个监测点。定义方位角θ为爆心连线与装药运动方向的夹角，即爆心连线与x轴正方向的夹角。

图  1  参数设置

Figure  1.  Parameter setting

下载: 全尺寸图片 幻灯片

为能较好地得到贴合实际的动爆冲击波特性，分别计算v为0、272、340、680、1 020和1 700 m/s等6种速度条件下的冲击波场，分别对应0Ma、0.8Ma、1.0Ma、2.0Ma、3.0Ma和5.0Ma的情况。

1.2   数值模拟结果

分析装药动爆冲击波流场演化图（见图2），可以发现：不同运动速度的装药在无限空中爆炸形成的冲击波均以球面波形式扩展；当v=0 m/s（静爆）时，压力分布呈现规则的球形，以x=0为对称轴呈对称结构，超压峰值约为0.37 MPa，因此在静爆条件下球形裸装药的爆炸冲击波场是规则的球对称结构；当v=272，340，680，1 020 m/s（动爆）时，冲击波压力分布随着装药运动速度的增大而严重畸变，与速度方向夹角越小的地方压力峰值越大，相同运动速度时的冲击波超压峰值出现在装药运动速度方向上。对比爆轰产物分布状态：当v=0 m/s时，TNT材料分布以x=0为对称轴呈左右对称状态，形成蝴蝶状；随着速度的增加，爆轰产物分布逐渐变形，形状变得不规则；当v=1 700 m/s时，超压峰值达到了0.96 MPa，是v=0 m/s（静爆）时超压峰值的2.6倍，此时的爆轰产物分布则不呈现蝴蝶状，形状不规则。

图  2  当t=0.4 ms时不同速度下的冲击波压力云图

Figure  2.  Pressure contour of shock waves under different velocities at t=0.4 ms

下载: 全尺寸图片 幻灯片

为了验证数值模拟的有效性，将速度v=0 m/s、方位角θ=0°的数据与Henrych公式计算结果^[8]进行比较（见表1）。数值计算的结果与经验公式相比较，误差都小于5%，这验证了数值计算结果的有效性和准确性。

表  1  超压数据与公式计算结果比较

Table  1.  Comparison of peak overpressure results with theoretical values

Δp/MPa		η/%
模拟	Henrych公式
2.49	2.38	4.41
0.97	0.94	3.78
0.47	0.47	1.07
0.28	0.27	2.22

下载: 导出CSV 

| 显示表格

本文中，主要研究起毁伤破坏作用的冲击波超压峰值Δp、冲击波正压作用时间τ、冲击波比冲量I等3个威力参数。从图3可以发现，当比例距离$\overline R$变大时，不同装药速度的超压越来越小，冲击波超压随$\overline R$增大衰减得较快。从超压的角分布来看：当θ=90°时，不同装药运动速度的动爆冲击波超压曲线汇聚于一点，数值大小与静爆时相近；当θ＞90°时，超压曲线随装药速度大小由下至上排列，可以得到装药速度越大冲击波超压越小的结论；当θ＜90°时，超压曲线随装药速度大小由上至下排列，可以得到装药速度越大冲击波超压越大的结论。因此可以得到结论：当90°＜θ＜180°、比例距离一定时，装药速度越大冲击波超压峰值越小；当0°＜θ＜90°、比例距离一定时，装药速度越大冲击波超压峰值越大。这个结果与冲击波流场演化图中观察到的结果一致。

图  3  不同比例距离时超压峰值随方位角的变化

Figure  3.  Peak overpressure Δp changes with θ under different $\overline R$

下载: 全尺寸图片 幻灯片

为了描述对称分布的监测点测得的冲击波超压之间的大小关系，引入一个增大系数r，定义为以中轴线为对称轴对称分布、相同比例距离处的监测点测得的超压峰值之间的比。

增大系数描述了装药运动速度导致θ＜90°方向与θ＞90°方向的冲击波超压峰值变化程度。由图4可以看出，当速度值小于2.0Ma时，运动导致的超压峰值增幅小于2倍，然而当速度大于2.0Ma时，增幅愈加明显。从2.0Ma开始，速度增幅分别为50%和150%，然而增大系数的增幅达到78%和418%，超压增幅程度远大于速度增幅程度。后两个对称方向的增大系数均小于θ=0°方向与θ=180°方向，说明装药速度一定时，超压峰值沿θ=0°至θ=180°呈现一个由大到小的分布状态。

图  4  不同方向的增大系数

Figure  4.  Increasing ratios r in different directions

下载: 全尺寸图片 幻灯片

如图5所示，比冲量也明显呈现出随方向角分布的特点：相同比例距离下，比冲量随角度的增大而减小；大于90°，装药速度越大比冲量越大；小于90°，装药速度越小比冲量越小。

图  5  不同比例距离时比冲量随方位角的变化

Figure  5.  Specific impluse I changes with θ under different $\overline R$

下载: 全尺寸图片 幻灯片

总体上，有装药速度越大、正压作用时间越大的趋势。由图6可见：当θ＜90°时，有装药速度越小、正压作用时间越大的趋势；当θ＞90°时，有装药速度越大、正压作用时间越大的趋势。此外，当比例距离较小时，正压作用时间的变化规律不是很明显，可以知道距离爆炸中心近区的物理参数变化情况比较复杂，变化的规律性不强。

图  6  不同比例距离时正压时间随方位角的变化

Figure  6.  Position pressure time τ changes with θ under different $\overline R$

下载: 全尺寸图片 幻灯片

1.3   建立计算模型

对动爆冲击波场进行简单分析，刚形成的冲击波可以近似为以爆炸中心为球心的球形冲击波。球面上的冲击波速度可以通过静爆冲击波速度${{{v}}_{\rm s}}$与装药运动速度${{v}}_0$的矢量相加，近似求出该处的动爆冲击波速度：${{{v}}_{\rm s}} = {{{v}}_{{\rm s} }}+{{{v}}_0}$。显然，随着方位角增大（0°＜θ＜180°），冲击波速度逐渐减小，因此冲击波超压大小变化与角度相关。这与前文的仿真数据分析得到的结论一致。

在描述动爆超压时可以引入修正因子δ^[7]，定义δ为：

修正因子δ是与装药运动速度v、方位角θ和比例距离$\overline R$等3个因子相关的量，即：

式中：f (v)为在特定的方位角θ^*和特定的比例距离$\overline R ^*$处(θ^*和$\overline R ^*$均为任取)、不同装药速度v时的δ；f ($\overline R$)为在特定的方位角θ^*、取一系列装药速度v时、不同比例距离$\overline R ^*$处的δ与$\overline R$ = $\overline R ^*$时的δ之比；f (θ)为取一系列装药速度v和比例距离$\overline R$时、不同方位角θ处的δ与θ=θ*时的δ之比。

取θ^*=0°、$\overline R ^*$ = 0.903 m/kg^1/3，得到不同装药运动速度下的修正因子δ，可拟合得到该平均值随装药运动速度v的关系曲线（见图7），即函数f(v)为：

图  7  拟合曲线f (v)

Figure  7.  Fitting curve f (v)

下载: 全尺寸图片 幻灯片

式中：装药速度v的单位为m/s，c₀为声速，取340 m/s。

取θ^* = 0°，得到一系列不同装药运动速度v、不同比例距离$\overline R$处的修正因子δ。在各种不同装药运动速度下，对不同比例距离$\overline R$时的修正因子δ取平均值，可拟合得到该平均值随比例距离$\overline R$的关系曲线（见图8），即函数f($\overline R$)为：

图  8  拟合曲线f ($\overline R$)

Figure  8.  Fitting curve f ($\overline R$)

下载: 全尺寸图片 幻灯片

式中：比例距离$\overline R$的单位为m/kg^1/3。

取一系列不同装药运动速度v和比例距离$\overline R$时，画出不同比例距离$\overline R$处的修正因子δ与方位角θ之间的关系曲线。在各种不同装药运动速度和比例距离下，对不同方位角θ时的修正因子δ取平均值，可拟合得到该平均值随方位角θ的关系曲线（见图9），即函数f (θ)为：

图  9  拟合曲线f (θ)

Figure  9.  Fitting curves f (θ)

下载: 全尺寸图片 幻灯片

式中：方位角θ的单位为rad。

得到修正因子δ的表达式为：

式中：修正因子的适用范围为$0.5\;{\rm{m}}/{\rm{k}}{{\rm{g}}^{1/3}} {\text{≤}} \overline R {\text{≤}} 1.5\;{\rm{m}}/{\rm{k}}{{\rm{g}}^{1/3}}$。将δ与各种经验公式组合便能得到不同公式的动爆冲击波计算模型。

1.4   计算模型验证

为了验证更高的可靠度，仍选用TNT球形裸装药进行计算，TNT装药质量W=0.170 1 kg，装药密度ρ=1.63 g/cm³。计算得到装药运动速度v为100、500、1 500 m/s时的动爆超压场的压力分布。分别选取比例距离$\overline R$为0.72、1.81、1.26、0.54 m/kg^1/3，方位角为θ为15°、73°、164°处的10组超压数据，见表2。

表  2  仿真计算结果对比

Table  2.  Comparison of simulation results

v/(m·s−1)	$\overline R$/(m·kg−1/3)	θ/(°)	Δp/kPa			η/%
			公式	静爆	模拟
100	0.72	15	1 609.03	1 516.35	1 663.80	3.29
100	1.81	15	214.91	219.84	236.42	9.09
100	1.26	73	394.31	480.98	493.13	20.04
100	0.54	164	1 155.98	2 552.41	2 322.79	50.23
500	0.72	15	2 324.81	1 516.35	2 367.87	1.82
500	1.81	15	310.52	219.84	317.83	2.30
500	1.26	73	569.71	480.98	533.51	6.79
500	0.54	164	1 670.22	2 552.41	1 579.85	5.72
1 500	0.72	15	5 551.67	1 516.35	5 249.82	5.75
1 500	1.81	15	741.52	219.84	690.90	7.33

下载: 导出CSV 

| 显示表格

采用运动装药爆炸冲击波的测试数据^[9]，对公式进行验证，见表3。实验所采用的是质量为0.17 kg的B炸药，其与TNT当量换算系数取1.3。

表  3  实验结果对比

Table  3.  Comparison of experimental results

v/(m·s−1)	$\overline R$/(m·kg−1/3)	θ/(°)	Δp/kPa		η/%
			实验	公式
534.31	1.37	15	627.42	638.32	1.71
534.31	1.37	45	565.37	572.21	1.20
534.31	1.37	105	361.97	375.33	3.56
587.65	1.19	15	751.53	926.03	18.84
587.65	1.19	45	599.84	830.13	27.74
587.65	1.19	105	379.21	544.50	30.36

下载: 导出CSV 

| 显示表格

静爆冲击波由Henrych公式计算^[8]：

式中：Δp_m的单位为MPa，${\overline R }$的单位为m/kg^1/3。

通过误差分析，得到比例距离、角度、装药速度不全相同的10组数据的公式验证，其中8组的误差小于10%，模型计算结果基本达到预期效果。在与实验数据的检验中，误差的平均值小于14%，考虑到当时实验器材精度状况以及实验误差，这个结果符合得比较好。

2.   结　论

通过使用AUTODYN软件进行仿真计算，分析动爆冲击波场的变化，研究了冲击波超压峰值、冲击波正压作用时间、冲击波比冲量3个重要参量，得出了相应的变化规律。主要有：装药速度对小于90°的方向的冲击波超压、比冲量起增强的作用，对90°到180°的方向的冲击波超压、比冲量起减弱的作用；而装药速度对小于90°的方向的冲击波正压作用时间起减小作用，对于大于90°的方向的正压作用时间起增大作用；并且基于仿真结果进行分析，得到了动爆冲击波的一个工程计算模型。将本文得出的计算模型与软件仿真算出的不同装药速度下的结果作对比。读取不同比例距离、角度处测得的冲击波峰值超压，并将数据与根据模型计算得出的动爆超压峰值以及与实验测试数据对比分析其误差，结果表明计算模型能得到较好的结果，因此本次研究得出的计算模型具备一定的应用价值。