不同真空度下空中爆炸近场特性的数值模拟研究

空中爆炸是爆炸力学研究中的一个基本现象，深入分析爆炸空气冲击波的传播规律，特别是凝聚炸药空中爆炸时的近场特性，对于研究武器试验、爆破工程、炸药工业等领域的安全防护问题有重要意义，如露天爆破、爆炸焊接时冲击波等有害效应的控制，爆炸容器等结构的抗爆防护设计等。

众所周知，空中爆炸时，爆轰产物高速膨胀强烈压缩周围的空气形成空气冲击波，而表征空中爆炸近场特性的特征参量包括峰值超压Δpm、比冲量i、正压作用时间t+、冲击波到达时间t′等(如图1所示)，为计算这些参量，学者们通过大量试验提出了各自的计算公式[1-8]，这些经验公式具有不同的适用范围，相互之间也存在一定的差异。

另一方面，当外界空气环境(如压力、温度、湿度、风力等)改变时，上述特征参量的规律将会随之发生变化[9]，即满足萨克斯比例定律[10]，Dewey等[11]通过调节高度模拟器中的环境压力和温度，测定了多种尺寸裸彭托利特炸药的爆炸波参数，得出的结论证实了萨克斯定律，该定律同样由Ericsson等[12-13]的实验所证实；Veldman等[14]研究了环境压力对反射冲击波超压和冲量的影响；Silnikov等[15]对比分析了高空中和常压环境下爆炸冲击波的参数变化及对飞行器的破坏影响。另一方面，Jack等[16]在炸药近旁和模拟高空情况时的试验指出，上述条件下爆炸波的性质将会发生改变，萨克斯定律此时不能用于特征参量的换算，如何准确而又完整地描述爆炸近场特征参量的变化规律是本文将要解决的主要问题之一。

图1 典型的空中爆炸p-t曲线
Fig.1 A typical pressure-time curve for air explosion

AUTODYN有限元程序可很好地用于模拟爆炸冲击波的传播过程[17-20]，本文将基于该平台建立空中爆炸的球对称楔形模型，在与标准状态经验公式对比的基础上，引入特征比例距离分析不同真空度下空中爆炸近场的特征参量变化规律，结论可为高空爆炸时目标的毁伤效应以及爆破工程的冲击波防护提供参考。

1 计算模型与状态方程

1.1 计算模型

本文采用AUTODYN程序中的球对称楔形模型来模拟自由场中球形装药爆炸冲击波的传播，建立如图2所示的计算模型。模型采用多物质Euler算法，选用与实际对应的ANFO炸药，炸药半径与网格尺寸之比λ=100，根据文献[18]和[19]的分析结果，所用计算模型和网格尺寸足以保证计算精度，从炸药边界开始设置一系列监测点(编号1#～96#)用于特征参量的结果输出。

图2 自由场空中爆炸楔形计算模型
Fig.2 Wedge model of free-space air explosion

1.2 参数设置

1.2.1 炸药

采用AUTODYN标准材料模型库里的ANFO炸药，其状态方程为经典的JWL状态方程，具体为

(1)

式中：P为爆压，V为爆轰产物相对比容(爆轰产物体积与炸药初始体积之比)，E0为炸药单位体积的初始内能，A、B、R1、R2、ω为炸药的特征参数(常数)，本文ANFO炸药各参数取值如表1所示。

1.2.2 空气

空气的状态方程近似采用理想气体状态方程，可描述为

表1 ANFO炸药参数
Tab.1 Parameters of ANFO

(2)

式中：P为空气压力，γ为理想气体绝热指数，取1.4；ρ为空气密度；ρ0为标准状态下的空气密度，取1.225×10-3 g/cm3；E0为空气的初始比内能，取2.068×105J/m3。本文通过改变空气的初始密度近似实现不同的初始压力，即实现不同真空度的空气环境。

2 爆炸相似律与经验公式

为检验模拟结果的可靠性，通过对比经验公式计算的结果分析数值计算的精度。根据萨克斯比例定律可知，包含超压Δpm、比冲量i、正压时间t+和某些空气参数的无量纲群

特定地的是(Rp01/3/E1/3)的单一函数。其中，R为爆心距离，p0为外界空气初压，c0是环境空气的声速，E是炸药的总能量(通常用装药质量W代替)。当炸药种类、装药密度和空气状态不变时，即为一般形式的爆炸相似定律，也称霍普金森比例定律，即

是比例距离的单一函数。

基于以上规律，学者们通过大量试验各自得到了计算冲击波峰值超压Δpm、正压作用时间t+和比冲量i的经验公式，但这些经验公式的适用范围和计算结果互有差异，以爆炸波峰值超压为例，图3给出了不同冲击波峰压经验公式计算得到的结果与数值模拟结果之间的关系(R0为装药半径)。

图3 不同冲击波峰压经验公式与模拟结果的关系

Fig.3 Comparison of simulated peak pressure with calculated results by empirical formulas

从图3中可以看出：

① 各经验公式在R>10～12R0后的计算结果趋于一致，因此经验计算不适用于近场，其适用范围一般为R>10～12R0；

② 数值计算结果与J. Henrych公式最为吻合，在R>6R0时，误差在10%以内，在R>12R0时，误差小于5%；

③ 在装药近旁(R<1.6R0)，峰压与距离之间存在其他对应关系。

图4和图5则分别给出了利用不同经验公式计算得到的正压作用时间和比冲量。

图4 不同经验公式的正压作用时间计算结果
Fig.4 Calculated results of positive pressure duration by empirical formulas

可以看出，不同经验公式关于正压作用时间和比冲量的计算结果存在很大差别，无法进行参照对比，因此这两个参量的模拟结果不适合与经验公式作比较，文献[21]中通过数值计算得到了如图6和图7所示的不同密度PETN炸药的比冲量和正压作用时间变化曲线，实线1代表密度1.6 g/cm3的PETN，虚线2密度为0.4 g/cm3，此处Re=(E/pM)1/3，tM=(E/pM)1/3/(pM/ρM)1/2，iM=pMtM。图中同时给出了本文常压下ANFO炸药的模拟结果，虽然与文献选用的炸药种类不同，但两者在变化规律上基本一致。

图5 不同经验公式的比冲量计算结果
Fig.5 Calculated results of specific impulse by empirical formulas

图6 比冲量的数值计算结果对比
Fig.6 Comparison of simulated specific impulse curves

图7 正压作用时间的数值计算结果对比
Fig.7 Comparison of simulated positive pressure duration curves

通过以上分析和对比可知，本文采用的数值方法是可行的，利用AUTODYN有限元程序计算空中爆炸的结果具有较高置信度。

3 不同真空度下的模拟结果与分析

3.1 真空度对峰值超压的影响

为研究外界空气压力(即真空度)对凝聚炸药爆炸峰压的影响，将外界空气域的初始压力依次从1 atm，0.5 atm，0.2 atm，0.1 atm改变至完全真空，考察同一装药下20R0范围内的峰值超压变化，以比例距离为横坐标，峰压取对数坐标，绘制曲线群如图8所示。

图8 不同真空度下凝聚炸药爆炸的峰值超压曲线

Fig.8 Peak overpressure curves of condensed explosive at different vacuum degrees

通过观察图8各曲线发现，爆炸峰压的变化呈现明显的区间性，装药近旁和远处具有不同的变化规律，峰压首先沿p0=0 atm(即完全真空)曲线变化，在某一距离后，爆炸峰压的衰减即为空气冲击波的衰减，这是因为装药近区的特征参量受爆轰产物的影响。图9显示了1标准大气压下距爆炸中心不同距离处压力随时间变化的曲线群，在冲击波传播的初始阶段(4#Gauss点)，压力出现了两个峰值，分别对应爆轰产物-空气界面前端的冲击波峰值和界面后端的产物压力峰值，在最初阶段，空气冲击波峰压小于爆轰产物峰压(图8中临界点前的初始段)，随着爆轰产物的急剧膨胀，冲击波峰压迅速超过产物峰压，因此图8中各真空度峰压曲线的大部分体现的是空气冲击波峰压的衰减过程。这样，通过对爆轰产物自由膨胀曲线(p0=0 atm)和空气冲击波峰压衰减曲线分别拟合，就可以得到爆炸场全范围内峰值超压的表达式。

图9 不同距离处压力-时间变化曲线
Fig.9 Pressure-time curve at different distances

从图9中还可以看出，爆轰产物对压力的影响逐渐减小，图中56#监测点距离爆炸中心12R0，该处的压力衰减基本趋于平滑，此处也正是爆轰产物极限膨胀体积的位置，而大部分经验公式的适用范围也正是在极限膨胀体积之后。

为分析外界初始空气状态对爆炸特征参量的影响，我们定义特征比例距离

(3)

式中:fd=(ρ0/ρM)1/3=[(p0/pM)/(T0/TM)]1/3，称为初始状态比例系数，下标0表示初始空气状态，M表示参考状态(pM=0.101 325 MPa，TM=273 K)；R为爆心距离，m；W为炸药TNT当量，kg。

图10给出了以特征比例距离Z表示的不同空气环境的峰压变化规律，计算结果与相关实验结果[22]符合良好。需要说明的是，由于本文旨在研究初始压力的影响，因此计算结果中空气的初始温度都取定为288 K。

图10 以特征比例距离表示的冲击波峰压曲线

Fig.10 Peak overpressure of blast wave with the scaled distance Z as variables

考察图10中本文计算结果可以发现，在该坐标系中，对于同一炸药，不同初始空气压力的冲击波峰压衰减趋于同一变化规律，拟合后的方程形式为

(4)

式中:Δpm为峰值超压，MPa；p0为空气域初始压力，MPa。对于本文所选炸药，α=5.565，β=1.078。

分析图8中p0=0 atm的曲线，爆轰产物自由扩散时峰压的变化规律可拟合为

(5)

式中:Δpm为峰值超压，为比例距离，m/kg1/3。本文中，λ=-3.902，ξ=5.970。

于是，我们就可以得到不同初始空气压力下全范围峰压的分段函数

(6)

其中，α、β、λ和ξ由炸药的性质所决定，对于AUTODYN程序提供的ANFO炸药，上述系数分别为5.565、1.078、3.902和5.970。Rs为不同初始压力对应的临界距离，它与真空度的关系如图11所示。

图11 不同真空度的峰压临界距离
Fig.11 Critical distance vs vacuum degree

拟合后得

(7)

式中:R0为装药半径，为标准大气压，0.101 325 MPa。

3.2 比冲量的变化规律

图12给出了1 atm、0.5 atm、0.2 atm和0.1 atm下同一装药爆炸时比冲量的变化规律，在以Z为横坐标，比例冲量为纵坐标的坐标系中，不同真空度下的冲量变化都是非单调的，起初各比冲量相继达到极小值，随后在Z>0.4后便具有相同的变化轨迹，在Z=0.7时统一升高到一个极大值，最后比冲量又不断减小。比冲量的这种变化规律可为抗爆安全防护提供参考，可综合防护效果和工程造价将结构尺寸设计在比冲量极小值位置。

图12 不同真空度下的比冲量变化
Fig.12 Specific impulse curves of various vacuum degrees

为具体分析比冲量极值点的位置，我们以1标准大气压时的情况为例，给出了八个相继时刻下20R0范围内的相对压力分布(如图13所示)。

图13 1 atm下不同时刻的压力分布曲线
Fig.13 Pressure distributions at different times

首先，我们对t=0.085 ms时刻(分布4)进行单独考察。图14给出了该时刻下的压力分布曲线及云图，图中产物-空气界面前是空气冲击波阵面，而在爆轰产物内部，从界面反射回的稀疏波又从爆心产生了又一次反射，稀疏波尾部的压力达到最小值，而主冲击波后的压力还没能正常衰减到稀疏波尾部的状态，因此在稀疏波尾部就形成了如图所示的间断面，该间断面最初是在运用特征线法进行爆炸数值计算中发现的[23]，一些文献中阐述此间断面时称其为“二次冲击”(secondary shock)[24-25]。

图14 不同真空度下的比冲量变化
Fig.14 Specific impulse curves of various vacuum degrees

继续观察图13，间断面在初始阶段由于爆轰产物的推动作用由内向外移动，直到分布7(间断面到达约10R0处)之后才开始向爆心运动。图12中，曲线波谷对应图13中的分布4，该时刻稀疏波尾部的压力等于初压p0；而曲线的波峰对应图13中的分布7，该时刻则是间断面传播方向的转折点，即由向外移动开始向爆心运动。其他真空度下的比冲量规律类似，这种非单调性变化都反映了爆轰产物中稀疏波尾部间断面形成和传播的波动特性。

3.3 正压作用时间

不同真空度下正压作用时间的变化规律如图15所示，可以发现，正压作用时间随比例距离Z呈分段变化规律，各折点位置同样与稀疏波尾部的间断面相关，其中A点对应图13中的分布4，即此时稀疏波尾部压力等于p0；在此之前，真空度小的曲线始终位于真空度大的曲线上方，此后，各曲线重叠在一起，即可用特征比例距离Z为变量的同一方程来描述不同真空度下的正压作用时间变化规律。折点B则对应图13中的分布7，即此时间断面开始向爆心运动；折点C对应爆轰产物膨胀的极限体积处，而在BC之间的范围内，正压作用时间与特征比例距离呈线性关系。

图15 不同真空度下正压作用时间的变化规律
Fig.15 Positive pressure duration curves of various vacuum degrees

3.4 产物界面和波阵面的运动

通过程序给定的监测点的超压数据可以得到峰压的时程曲线，即为空气冲击波波阵面的运动轨迹，如图16所示；而通过追踪爆轰产物界面的运动轨迹，则可以得到产物界面运动的时程曲线，如图17所示，两者经特征化后的曲线也相应地给出，图中tM含义与前述相同(能量E用质量W代替)。

图16 空气冲击波波阵面时程曲线
Fig.16 Time-history curves of shock wave front

分析图16和图17可知，在常规坐标系中，真空度越大，冲击波阵面和产物界面运动速度越快，且在一定距离后，冲击波速度各自趋于一恒定值，产物界面速度则趋于0，即爆轰产物膨胀至极限体积，而完全真空状态时，爆轰产物膨胀速度为一恒定值(U∞≈6 600 m/s)。不同真空度下产物界面和冲击波阵面时程曲线在相应的特征坐标系中都各自归结为一条运动轨迹，即在以Z表示的系统中，尽管外界的真空环境不同，但产物界面和冲击波阵面的运动规律都是相同的。

图17 爆轰产物界面时程曲线
Fig.17 Time-history curves of detonation product interface

不同真空度下爆轰产物的极限膨胀体积与初始空气压力的关系如图18所示，可以推算，它的垂直渐近线为Y轴，水平渐近线为X轴；而极限体积与临界距离之间呈线性关系(图19)，即

(8)

图18 极限膨胀体积与初始压力的关系
Fig.18 Ultimate expansion volume vs. initial pressure

图19 极限膨胀体积与临界距离的关系
Fig.19 Ultimate expansion volume as a function of critical distance

4 结 论

本文基于有限元程序AUTODYN-2D，建立了空中爆炸的一维球对称楔形模型，通过改变空气域的初始参数，得到了不同真空度下爆炸近场特征参量的变化规律，主要结论如下：

(1)常压下的模拟结果与峰压经验公式符合良好，尤其与J. Henrych最为吻合，在R>12R0时，误差小于5%，因此数值计算具有较高的置信度。

(2)不同真空度下爆炸峰值超压的衰减规律呈区间性，其中爆轰产物峰压对应完全真空时的衰减曲线，可用变量来表达，而不同真空度下的空气冲击波峰压则可在以特征比例距离Z表示的坐标系中用同一方程描述，这样，用以变量和Z表示的分段函数可表述全范围内爆炸峰压的变化规律，该方法同样适用于其他凝聚炸药，通过当量换算可确定对应系数。

(3)不同真空度下比冲量i和正压作用时间t+的变化具有非单调性，它们的极值点和折点与稀疏波尾部间断面密切相关，当稀疏波尾部压力等于p0时，i和t+达到极小值，而当间断面的传播方向由外向爆心转变时，i达到极大值，t+则出现折点。不同真空度在稀疏波尾部压力等于p0之后具有相同的比冲量和正压作用时间；此外，不同真空度下的爆轰产物界面和冲击波阵面时程曲线各自重合，极限体积与临界距离之间呈线性关系。