在建筑防爆设计、毁伤评估、火力筹划等应用场景中，往往要求在较短时间内进行大规模爆炸载荷计算，以评估意外爆炸、恐怖袭击等对建筑结构和人员的危害^[1-6]。计算流体动力学（computational fluid dynamics, CFD）方法因为计算密集性和需要高层次的专业知识较少使用。相较于CFD，理论方法计算简便且相对精确，在工程实践中得到了广泛应用^[7-9]。自由场爆炸波冲击波建模是所有爆炸载荷理论计算方法的基础。

现有理论方法可分为2类。第一类是UFC 3-340-02、ASCE/SEI 59-11等建筑抗爆炸防护标准，上述方法以Kingery-Bulmash (K-B)曲线为基础，依据建筑结构特点和与爆心位置，在K-B曲线中确定观测点的超压、到时、冲量等冲击波参数，并按照虚拟冲击波超压-时间关系（通常简化为三角波），确定观测点受到的爆炸冲击波载荷^[10-12]。该类工程模型的特点为目标载荷波形与实际情况相差较大，不考虑直接冲击波和多重反射波的耦合效应，其计算结果是保守的，超压计算误差可能达到100%以上^[12]。第二类是基于多冲击波叠加法则的计算方法，该类方法依据虚拟镜像法生成多重反射镜像爆源，观测点的载荷为直射波与多重反射波叠加耦合的结果，该类方法以表征任意时间、位置处冲击波三个基本流体动力学参数（压力、密度、粒子速度）的爆炸冲击波模型为基础，相较于第一类方法，该类方法计算结果的波形与实际波形接近，载荷计算结果更为精确，且理论上适用于任意复杂结构^[8-9,13-15]。另一方面，越来越多的研究成果表明，冲击波波形及次级冲击波对建筑结构、人员等目标的损伤也有重要影响^[5,16-17]，因此，第二类方法愈发受到工程人员们的青睐。

作为应用最为广泛、研究最为充分的高爆炸药，TNT常用作计算各类炸药爆炸参数的基准^[18-19]。当前，TNT自由场爆炸主冲击波压力随时间、空间的变化规律研究最为充分。Kingery等^[9]和Kinney等^[20]分别依据数值计算结果和大量实验数据，建立了爆炸冲击波超压峰值、到时、最大超压冲量和部分其他参数随比例距离的关系，上述关系得到了广泛应用。冲击波压力随时间变化关系常用修正Friedlander方程来表征^[21-23]，修正Friedlander方程中衰减系数可由冲击波超压峰值、到时、最大超压冲量、正压持续时间经牛顿割线法迭代计算^[24]。Xue等^[25]以数值模拟为基础，提出了预测不同时刻、距离的自由场爆炸冲击波超压的新的经验公式。对于密度、粒子速度建模问题，Dewey^[26-27]依据300多次化爆实验结果，建立了预测主冲击波正相超压、密度、粒子速度随时间、距离变化关系的Excel工具，但该工具不包含主冲击波负相及后续脉冲的计算。目前，尚无公开发表的、完整描述任意时刻、距离下TNT自由场爆炸冲击波压力、密度、空气速度的经验模型。Fl-Blast、BlastX等商业软件假设冲击波阵面后任意位置、任意时刻的三个基本流体动力学参数均满足Rankine-Hugoniot关系，进而计算流场密度、粒子速度。显然，以该假设得出的计算结果是不严谨的^[14]。因此，建立TNT自由场爆炸条件下，冲击波压力、密度、粒子速度随距离、时间变化的经验模型，进而支撑复杂结构爆炸载荷理论计算，具有一定的工程实践意义。

本文通过数值模拟的方法，建立计算TNT自由场爆炸冲击波压力、密度、粒子速度的经验模型，利用爆炸冲击波地面反射和爆炸冲击波建筑后绕射两个典型工况波后流场的理论计算，阐释提出模型的应用场景，并与实验、数值模拟数据作对比，验证提出模型和相关理论方法的准确性。

1.   TNT自由场爆炸冲击波建模

利用LS-DYNA软件对球形TNT装药自由场爆炸进行了数值模拟，获得了三个基本冲击波参数及其冲量、到时等参数的变化规律，建立了描述爆炸冲击波压力p、密度ρ、粒子速度u随半径r、时间t的变化经验模型。

1.1   数值模型

采用一维多材料任意拉格朗日-欧拉算法（multi-material arbitrary Lagrangian-Eulerian, MM-ALE）计算球形装药自由场爆炸冲击波参数。装药设置为TNT，起爆方式为中心起爆，质量W分别为1、10、100、1000、10000 kg。TNT装药密度为1630 kg/m³，因此装药半径对应的比例距离为R=0.053 m/kg^1/3（比例距离R=r/W^1/3, r为半径、W为TNT当量）。在比例距离为0.1～10 m/kg^1/3范围内均布150个测点，以记录冲击波参数随时间变化情况。计算模型的半径相当于比例距离16 m/kg^1/3，以减少透射边界对测点数据的影响，计算截止比例时间为45 ms/kg^1/3，以保证冲击波扫过所有测点。TNT采用JWL物态方程，其参数如表1所示。空气采用线性物态方程，初始密度为1.225 kg/m³，初始压力为101332 Pa，比热比γ为1.4。

表  1  TNT的JWL物态方程参数^[28]

Table  1.  JWL EOS parameters of TNT^[28]

密度/(kg·m−1/3)	E0/GPa	A/GPa	B/GPa	R1	R2	ω
1630	7.0	373.8	3.75	4.15	0.9	0.35

下载: 导出CSV 

| 显示表格

分别在40000、80000和160000个网格情况下进行数值模拟，以验证网格收敛性。图1～图3分别展示了1 kg TNT装药爆炸冲击波压力p、密度ρ、粒子速度u在比例距离1、4、7、10 m/kg^1/3下的时程曲线，可以看到三个网格数量下得到了几乎相同的计算结果，仅有细微差别，综合计算经济性及准确度需求，选择160000为计算网格数量，则装药半径上有527个网格，网格数量满足高精细数值模拟的要求^[29-30]。10、100、1000、10000 kg装药的网格收敛性分析结果类似，不赘述。

图  1  不同网格数量下的压力时程数值模拟结果

Figure  1.  Numerical result of pressure changing with different grid quantity

下载: 全尺寸图片 幻灯片

图  2  不同网格数量下的密度时程数值模拟结果

Figure  2.  Numerical result of density with different grid quantity

下载: 全尺寸图片 幻灯片

图  3  不同网格数量下的粒子速度时程数值模拟结果

Figure  3.  Numerical result of particle velocity with different grid quantity

下载: 全尺寸图片 幻灯片

图4分别为不同比例距离下冲击波正相超压峰值Δp⁺、到时t_a、正相冲量i⁺的数值模拟结果与K-B曲线^[9]、K-G曲线^[20]的对比，可以看到，数值计算结果与经验曲线基本吻合，近场处，正相超压峰值数值模拟结果略小于经验曲线，各测点误差均小于5%，冲量数值模拟结果局部误差较大，达到30%，为LS-DYNA MM-ALE算法固有缺陷导致。文献[7]对近场处多种经验方法、试验测试数据进行了不确定性分析，局部30%误差在可接受范围内。

图  4  数值模拟结果与经验曲线对比

Figure  4.  Comparison between simulated results and empirical curves

下载: 全尺寸图片 幻灯片

1.2   冲击波参数-比例距离关系

由霍普金森比例定律，正相超压峰值Δp⁺、负相超压峰值Δp⁻可表示为比例距离的函数，正相冲量i⁺，负相冲量i⁻（负相超压对时间的积分），超压峰值到时t_a、正相超压持续时间$t_{\text{d}}^ +$、负相超压持续时间$t_{\text{d}}^ -$等与时间的相关的参数，其比例数值i⁺/W^1/3、i⁻/W^1/3、t_a/W^1/3、$t_{\text{d}}^ +$/W^1/3、$t_{\text{d}}^ -$/W^1/3也可表示为比例距离的函数^[21]。

图5～图7上述参数与比例距离关系的数值模拟数据和拟合曲线（curve fitting method, CFM）。可以看出，正相超压峰值Δp⁺随着比例距离R的增大而减小。在0.1 m/kg^1/3<R<0.6 m/kg^1/3范围内，负相超压峰值Δp⁻近乎水平，原因为负相超压峰值时刻，该处空气几乎被冲击波排空，负相超压接近极值；在其余部分，负相超压峰值随着比例距离的增大而减小。正相冲量i⁺整体随比例距离的增大而减小，但在0.3 m/kg^1/3<R<0.7 m/kg^1/3范围内趋势相反。正相冲量为超压时程在正相的积分，其数值与超压、正压持续时间两个参数相关，如图5(f)所示，正相持续时间在0.3 m/kg^1/3<R<0.7 m/kg^1/3范围内急剧上升，此时正向持续时间对正相冲量的增强效应占主导；而在其余位置，超压随比例距离减小占主导，该现象也可与文献[10,18]结果吻合。负相冲量i⁻与比例距离的关系更为复杂，也与负相超压、负压持续时间两个参数对负相冲量的贡献有关。超压峰值到时t_a与比例距离正相关。超压正相持续时间$t_{\text{d}}^ +$、超压负相持续时间$t_{\text{d}}^ -$整体与比例距离正相关，两个参数在0.1 m/kg^1/3<R<0.8 m/kg^1/3范围内的剧烈变化与该处爆轰产物与主冲击波未完全分离有关，文献[9,31-32]计算认为爆轰产物与主冲击波分离的位置约为比例距离0.84 m/kg^1/3处，与上述结果吻合。

图  5  超压峰值与比例距离关系

Figure  5.  Relationship between peak pressure and scaled distance

下载: 全尺寸图片 幻灯片

图  6  冲量与比例距离关系

Figure  6.  Relationship between impulse and scaled distance

下载: 全尺寸图片 幻灯片

图  7  时间参数比例距离关系

Figure  7.  Relationship between time parameters and scaled distance

下载: 全尺寸图片 幻灯片

当计算截止时，冲击波仍在向外扩张，因此计算区域的空气密度并未回到初始值，定义计算截止时刻各位置空气密度为中止密度ρ_t。为后续建模方便，文中所述超密度均为相对中止密度的增量，而非终态密度（数值上等于初始密度），并定义超密度积分i_ρ（超密度沿时间的积分）、密度峰值到时t_aρ（密度峰值对应时刻）等参数。

图8～图10为中止密度ρ_t、正相超密度峰值Δρ⁺、负相超密度峰值Δρ⁻、比例正相超密度积分极值$i_\rho ^ +$/W^1/3、比例负相超密度积分极值$i_\rho ^ -$/W^1/3、比例超密度峰值到时$t_{{\text{a}}\rho }^{}$/W^1/3、比例正相超密度持续时间$t_{{\text{d}}\rho }^ +$/W^1/3、比例负相超密度持续时间$t_{{\text{d}}\rho }^ -$/W^1/3等参数与比例距离关系的拟合曲线及数值模拟数据。可以看出，超密度相关参数也服从霍普金森比例定律，中止密度ρ_t整体上随比例距离的增大而增大，并趋近于初始空气密度；在0.1 m/kg^1/3<R<0.8 m/kg^1/3范围内爆轰产物对中止密度的变化趋势扰动较大，在0.63 m/kg^1/3<R<0.69 m/kg^1/3范围内，中止密度出现向下“间断”，这是由于在该范围内粒子速度时程波形发生较大变化（由双峰变为单峰，参见文献[9,31-32]）。正相超密度峰值Δρ⁺随比例距离的增大而减小，如图5(b)所示，图中黑色曲线为依据Rankine-Hugoniot关系计算的冲击波阵面后超密度曲线，在0.1 m/kg^1/3<R<0.57 m/kg^1/3范围内，其明显小于数值计算结果。这是由于在该范围内，冲击波阵面后爆轰产物的超密度峰值大于冲击波阵面后的空气超密度峰值。在0.1 m/kg^1/3<R<1.3 m/kg^1/3范围内，中止密度为负相超密度峰值Δρ⁻、正相超密度积分i_ρ、负相超密度积分$i_{\rho}^-$、正相超密度持续时间$t_{{\text{d}}\rho }^ +$、负相超密度持续时间$t_{{\text{d}}\rho }^ -$等5个参数变化趋势的主导因素，当R≥1.3 m/kg^1/3，负相超密度峰值、正相超密度积分、负相超密度积分随比例距离增大而减小，正相和负相超密度持续时间随比例距离的增大而增大。明显的，超密度到时t_aρ随比例距离的增大而增大。

图  8  中止密度和超密度峰值与比例距离关系

Figure  8.  Relationship of cut suspend density and peak density with scaled distance

下载: 全尺寸图片 幻灯片

图  9  密度积分与比例距离的关系

Figure  9.  Relationship between density integral and scaled distance

下载: 全尺寸图片 幻灯片

图  10  时间参数与比例距离关系

Figure  10.  Relationship between time parameters and scaled distance

下载: 全尺寸图片 幻灯片

同样地，定义速度积分i_u（粒子速度沿时间的积分）、速度到时t_au（粒子速度峰值对应时刻）等参数。图11～图13为正相速度峰值Δu⁺、负相速度峰值Δu⁻、比例正相速度积分极值$i_u^ +$/W^1/3、比例负相速度积分极值$i_u^ -$/W^1/3、比例速度峰值到时t_au/W^1/3、比例正相速度持续时间$t_{{\text{d}}u}^ +$/W^1/3、比例负相速度持续时间$t_{{\text{d}}u}^ -$/W^1/3与比例距离关系的拟合曲线及数值模拟数据。Δu⁺随比例距离增大而减小，在0.1 m/kg^1/3<R<0.7 m/kg^1/3范围内，依据Rankine-Hugoniot关系计算的冲击波阵面后粒子速度小于正相速度峰值，这是由于在该范围内，冲击波阵面后爆轰产物粒子速度峰值大于冲击波阵面后空气速度。受爆轰产物、空气相互作用影响，0.1 m/kg^1/3<R<0.8 m/kg^1/3范围内，$i_u^ +$、$i_u^ -$、$t_{{\text{d}}u}^ +$和$t_{{\text{d}}u}^ -$波动较为剧烈。需要特别说明的是，上述参数整体上也服从霍普金森比例定律；但计数精度引起计算初始条件的微小变化，即导致负相速度积分$i_u^ -$数值计算结果的局部发散，表明在0.1 m/kg^1/3<R<0.4 m/kg^1/3范围内负相速度积分对计算初始条件极其敏感。该范围内计算结果与R>0.4 m/kg^1/3范围内有2个数量级的差距，因此采用一次多项式拟合即可。

图  11  粒子速度峰值与比例距离关系

Figure  11.  Relationship between peak particle velocity and scaled distance

下载: 全尺寸图片 幻灯片

图  12  速度积分与比例距离关系

Figure  12.  Relationship between velocity integral and scaled distance

下载: 全尺寸图片 幻灯片

图  13  时间参数比例距离关系

Figure  13.  Relationship between time parameters and scaled distance

下载: 全尺寸图片 幻灯片

依据数值模拟结果，采用最小二乘法曲线拟合获得1.2节所述爆炸冲击波参数与比例距离的关系，拟合过程中所有数据均转化为国际单位制。使用的函数关系有多项式、有理函数、指数函数、高斯函数，其通用形式分别为

式中：P_ara为爆炸冲击波参数或其缩尺值，m、n为函数阶数，a、b、c为拟合系数。

以决定系数度量拟合效果，使得拟合结果的决定系数值大于0.999，则拟合曲线能很好地解释计算数据的规律。决定系数最小处为负相速度冲量拟合结果，与图12(b)局部数据较离散对应。

1.3   冲击波参数-时间关系

以改进修正Friedlander方程描述压力、密度、粒子速度三个冲击波基本参数与时间的关系。改进修正Friedlander方程的基本形式为

式中：p₀为环境压力；e、f、G为待定参数，可利用牛顿割线法迭代求解^[24]，求解条件为

改进修正Friedlander方程（式(5)）不仅包含了主冲击波超压时程关系，也包含了次级冲击波超压时程关系，更适用基于LAMB多冲击波叠加法则的第二类冲击波载荷理论计算方法。

当比例距离较小时（R≤0.57 m/kg^1/3），空气冲击波阵面后密度小于爆轰产物密度峰值，如图6(b)所示。为准确描述冲击波后密度与时间的关系，式(5)改写为

即当t_a≤t<t_aρ时，依据冲击波阵面后空气超压，通过Rankine-Hugoniot关系近似估计密度，γ=1.4与数值模拟空气参数相同。为保持形式一致，t≥t_aρ时，以中止密度ρ_t为基准，构建修正Friedlander方程。当R>0.57 m/kg^1/3时，空气冲击波阵面密度大于爆轰产物密度峰值，即t_a=t_aρ，此时密度时程波形简化为Friedlander方程指数衰减率。

同样的，构建描述粒子速度与时间关系的改进修正Friedlander方程为

当R<0.7 m/kg^1/3时，t_a<t_au，即空气冲击波阵面后粒子速度小于爆轰产物粒子速度峰值；当R≥0.7 m/kg^1/3时，t_a=t_au。式(10)和式(11)中待定参数计算方法与式(5)类似。

式(5)、(10)和(11)与冲击波参数-比例距离经验方程联立，即可计算任意距离、任意时刻压力、密度、粒子速度等三个基本冲击波参数。图14为R=0.17 m/kg^1/3处冲击波基本参数时程曲线的理论及数值模拟结果对比，为方便观察，密度曲线、速度曲线分别向后移动0.125和0.25 ms。可以看到，爆轰产物尚未与主冲击波分离，近场压力、密度、速度时程曲线均非典型的指数衰减，可呈现双峰甚至三峰形态，以改进修正Friedlander方程计算的理论曲线与数值模拟结果在波形上相差较大。同时，冲击波参数时程曲线在不同比例距离处差异极大，很难建立统一的数学形式。在R<0.7 m/kg^1/3范围内，改进修正Friedlander方程参数求解过程保证了参数峰值及其冲量与数值模拟结果相等，对于超压-冲量类毁伤判据目标的载荷计算已经足够。

图  14  R=0.17 m/kg^1/3处冲击波基本参数理论及数值模拟时程

Figure  14.  Simulated and theoretical basic shock wave parameter histories at R=0.17 m/kg^1/3

下载: 全尺寸图片 幻灯片

图15为爆轰产物扩张极限半径外压力、密度、速度时程曲线的理论及数值模拟结果对比，可以看到，不同比例距离下3个基本冲击波参数的时程曲线均遵循指数衰减规律，改进修正Friedlander方程的计算结果与数值模拟结果基本吻合。另外，改进修正Friedlander方程算得的次级冲击波时程曲线为平滑过渡，而数值模拟结果为阶跃上升，但其峰值基本一致。

图  15  冲击波基本参数理论及数值模拟时程曲线

Figure  15.  Simulated and theoretical basic parameter history of shock waves

下载: 全尺寸图片 幻灯片

2.   应用及验证

利用文献[33]TNT空爆地面反射和文献[34]TNT爆炸冲击波建筑后绕射两个工况冲击波参数的理论计算，阐释提出模型的典型应用场景，并与数值模拟、试验测试结果做对比验证。

2.1   空爆地面反射

2.1.1   实验布置

图16为空爆地面反射实验布置，9.28 kg球状TNT裸装药置于距地面高3 m处，以爆心在地面的投影为零点，4个超压测点T1、T2、T3、T4的坐标分别为(4.8 m, 3 m)、(5.8 m, 3 m)、(4.8 m, 0.2 m)和(5.8 m, 0.2 m)。采用LS-DYNA软件对该场景下爆炸冲击波传播过程进行二维精细数值模拟，TNT、空气材料参数与1.1节相同，空气域20 m×20 m，网格数量1600万，地面设置为刚性壁面。为提高计算精度，采用一维到二维冲击波载荷映射技术，一维计算至冲击波接触刚性壁面前截止。

图  16  空爆地面反射实验布置

Figure  16.  Experimental arrangement of ground reflection of air explosion

下载: 全尺寸图片 幻灯片

2.1.2   理论计算方法

假设空爆地面反射冲击波为爆源与虚拟爆源冲击波的叠加，虚拟爆源位置为反射平面另一侧，到反射面的距离为爆高H_C。按冲击波传播过程将刚性壁面反射冲击波流场划分为三个区域，入射波控制区域A、反射波控制区域B、马赫波控制区域C。正规反射时，区域B的边界为以虚拟爆源为圆心，以入射波为半径的圆弧；马赫反射时，区域C的边界为马赫波阵面及三波点迹线，区域B的边界为三波点迹线和以虚拟爆源为圆心，虚拟爆源到三波点连线为半径的圆弧^[9,35]，如图17所示。

图  17  冲击波结构示意

Figure  17.  Schematic of shock wave structure

下载: 全尺寸图片 幻灯片

三波点轨迹采用文献[36]提出的经验公式计算，该公式适用比爆高为0.5～3.5 m/kg^1/3和比例距离为0～15 m/kg^1/3的范围，计算表达式为

式中：G₀为马赫反射起始点位置，H_T为三波点高度。

对区域A，其冲击波参数计算方法与自由场爆炸相同。正规反射时，对区域B，其冲击波参数为入射波与虚拟入射波的叠加，应用LAMB多冲击波叠加法则计算叠加后冲击波参数

式中：i为爆源编号，n为爆源数量，u为粒子速度的矢量形式。

马赫反射时，反射波半径大于入射波，需依据反射波半径与入射波半径之比对观测点与虚拟爆源距离进行等比例缩放，再计算虚拟爆源冲击波参数，随后应用LAMB法则。区域C冲击波参数计算与区域B类似。

2.1.3   压力云图对比

图18为理论与数值计算得到的压力云图对比，左侧为理论计算结果，右侧为精细数值模拟结果，两云图的色表相同。t=5.3 ms时，入射波在刚性壁面上发生正规反射，可以看到，理论计算结果与数值模拟结果的压力分布基本相同；理论方法假定反射波为圆弧状，而数值模拟反射波受爆轰产物影响，呈喇叭状，与理论结果略有不同。t=7.9 ms时，冲击波阵面已扫过测点T3，理论计算结果为正规反射，而数值模拟结果已出现马赫反射；这是由于理论法严格按照式(12)～(13)划分冲击波结构，而数值模拟采用MM-ALE方法计算冲击波结构，两者方法不同，得出局部相左结果。t=10.5 ms时，冲击波阵面扫过所有测点，理论、数值模拟计算的冲击波结构均进入反射阶段。t=31.6 ms时，冲击波阵面继续扩张，两者均出现3重次级反射波；相较于理论计算结果，数值模拟结果次级冲击波与负相界面极为清晰，表明该处压力值出现明显“间断”，而理论计算结果该界面处界面较模糊，存在色彩过渡，这是改进修正Friedlander方程次级冲击波时程曲线为平滑过渡导致的，与图15结果对应。

图  18  理论与数值计算得到的压力云图

Figure  18.  Theoretical and simulated results of pressure contour

下载: 全尺寸图片 幻灯片

2.1.4   测点数据对比

图19为各测点处实验、数值模拟和理论计算超压、冲量时程曲线对比。可以看到，理论计算与数值模拟得到的测点T1、T2冲击波超压、冲量时程曲线基本吻合，其超压峰值误差分别为2.5%、3.4%，冲量峰值误差分别为4.1%、5.3%。理论计算与实测超压峰值误差分别为23.3%、26.8%，冲量峰值误差分别为19.3%、33.6%。相较于实验数据，理论计算与数值模拟结果超压峰值更低且正压持续时间更长，推测与试验场海拔、试验当天温湿度、爆心距测试精度和要求加工精度有关，数值建模精度贡献的误差小于5%。理论计算与数值模拟得到的测点T3、T4冲击波超压、冲量时程曲线差距较大，理论计算结果的主冲击波为双峰形态，表明三波点位于测点T3、T4下方，而数值模拟结果的主冲击波超压时程为单峰形态，表明三波点位于测点T3、T4上方。实验测试结果主冲击波超压时程为双峰形态，印证了理论计算结果的准确性；数值模拟结果的误差可能与建模过程中将地面视为固壁，未完全还原试验条件有关。文献[33]认为激波绕射产生的稀疏波和波流导致测试反射波小于入射波，实际结果应反射波大于入射波，与理论计算结果吻合。

图  19  不同测点超压、冲量时程

Figure  19.  Overpressure and impulse histories of different gauges

下载: 全尺寸图片 幻灯片

该算例数值模拟采用Xeon 8358 CPU，32线程并行，耗时3 h 23 min；理论计算采用相同硬件配置，单线程计算，耗时11 min，表明该方法在计算时间上具有较明显的优势。

2.2   建筑外爆绕射

2.2.1   实验布置

图20为爆炸冲击波建筑后绕射实验布置俯视图，0.2 kg立方体TNT裸装药置于钢混地面建筑背侧，距背墙0.5 m，距地面1.35 m。建筑尺寸4 m×4 m×2.7 m，距建筑左墙0.35 m，距地面1.35 m间隔排列了5个测点，其中T6测点未启用。以建筑左下拐角为0点，超压测点T5、T7、T8、T9的坐标分别为（−0.35 m, −0.5 m, 1.35 m）、（−0.35 m, 0.5 m, 1.35 m）、（−0.35 m, 1 m, 1.35 m）和（−0.35 m, 1.5 m, 1.35 m）。实验共进行了3次，每次试验均采集了T5测点数据，而T7、T8、T9测点每次实验只启用其中一个。采用LS-DYNA软件对该场景下爆炸冲击波传播过程进行3维精细数值模拟，TNT、空气材料参数与1.1节相同，计算区域8 m×8 m×2.7 m，网格数量6400万，地面、建筑均设置为刚性壁面。为提高计算精度，采用一维到3维冲击波载荷映射技术，一维计算至冲击波接触刚性壁面前截止。

图  20  爆炸冲击波建筑后绕射实验布置图

Figure  20.  Experimental arrangement of explosion shock wave diffraction behind the building corner

下载: 全尺寸图片 幻灯片

2.2.2   理论计算方法

图21为冲击波绕射理论计算空间划分俯视图，虚拟爆源分别位于建筑背墙内部和地面另一侧（图中未展示）。建筑背墙、墙内虚拟爆源与建筑背墙左墙交线l形成的平面、地面围成区域D；爆源与l形成的平面、墙内虚拟爆源与l形成的平面、地面围成区域E；爆源与l形成的平面、建筑左墙、地面围成区域F。区域D冲击波由爆源入射波、墙面反射波和地面反射波叠加计算；区域E冲击波由爆源入射波、地面反射波叠加计算；区域F冲击波由爆源绕射波、地面反射波的绕射波叠加计算。由于各测点不在马赫波影响范围内，且没有复杂条件下马赫波结构计算先验知识，计算过程中不考虑马赫反射。

图  21  冲击波绕射理论计算空间划分示意图

Figure  21.  Schematic diagram of space division of shock wave diffraction theoretical calculation

下载: 全尺寸图片 幻灯片

以T9为例简述绕射波参数理论计算方法。通过射线追踪法获得爆源射线绕过墙角连接T9的拐点，爆源拐点距离为l₁、T9拐点距离为l₂，拐点向T9的射线与爆源向拐点的射线夹角为θ，则爆源绕射到T9实际距离为l₁+l₂，其等效距离R_e=(l₁+l₂)(1+sin θ)^{[9, 37]}。T9处爆源冲击波参数等效为半径R_e处自由场冲击波参数。

2.2.3   压力云图对比

图22为z=1.35 m平面理论与数值计算结果的压力云图对比，左侧为理论计算结果，右侧为精细数值模拟结果，两云图的色表相同。当t=1.6 ms时，爆源入射波阵面超越测点T5，压力云图上可以明显看到入射波与墙面反射波结构，理论计算结果与数值模拟结果基本相同。数值模拟结果出现马赫波结构，而理论计算过程未考虑马赫波结构，导致该处压力小于数值模拟结果。t=2.4 ms时，爆源绕射波接近T7测点，墙面反射波扫过T5测点，理论计算结果与数值模拟结果基本相同，冲击波绕射后由较明显的“滞后”现象；t=5.5 ms时，爆源绕射波接近T9测点，地面反射波进入z=1.35 m平面，在数值模拟压力云图中可以看出，地面反射波阵面附近出现类似马赫波的增强效应；t=8.7 ms时，地面反射波逐渐追上墙面反射波，理论计算结果与数值模拟结果中均可在D区域观察到明显的三波叠加结构，F区域观察到入射波、地面反射波绕射后的叠加结构。理论计算各区域交界处可看到明显的界限，而数值模拟结果为平滑过渡，更接近真实情况。数值模拟云图右侧、下侧可观察到明显的反射波，由透射边界程序漏洞导致，在计算截止前不影响测点处数据。

图  22  z=1.35 m平面理论与数值计算结果的压力分布云图

Figure  22.  Theoretical and simulated results of pressure contours on the z=1.35 m plane

下载: 全尺寸图片 幻灯片

2.2.4   测点数据对比

图23为各测点处实验、数值模拟和理论计算超压、冲量时程曲线对比。测点T5处理论计算结果、数值模拟结果与三次试验结果超压时程波形基本一致，均为主冲击波双峰，随后有较明显次级冲击波结构。三个波峰结构分别对应爆源入射波、墙面反射波和地面反射波，与文献[34]结论吻合。理论计算的超压峰值、正相冲量峰值与数值模拟结果误差分别为8.6%、0.8%，与三次试验结果平均值误差分别为19.8%、24.6%，具体误差数值如表2所示。流场发展本身是非线性极强的物理过程，鉴于理论方法所采用的诸多假定，理论结果与真实情况还存在一定误差^[35]。T9的冲量误差较大，原因为T9点实测曲线在主冲击波处产生双峰形态，属偶然误差^[34]。T7、T8、T9的理论计算结果、数值模拟结果与三次试验结果超压时程波形基本吻合，在7 ms内均出现两个波峰，分别对应入射波、地面反射波的绕射波，与文献[34]结论吻合。

图  23  不同测点超压、冲量时程曲线

Figure  23.  Overpressure and impulse histories of different gauges

下载: 全尺寸图片 幻灯片

表  2  各测点的超压峰值（Δ$p$⁺）和冲量极值（$i$⁺）理论结果与数值模拟、试验结果对比

Table  2.  Comparison of peak overpressure (Δp⁺) and maximum impulse (i⁺) between theoretical values and experimental results, simulation results of each gauges

测点	δ1(Δp+)/%	δ2(Δp+)/%	δ3(Δp+)/%	δs(Δp+)/%	δ1(i+)/%	δ1(i+)/%	δ1(i+)/%	δ1(i+)/%
T1	23.5	19.5	19.3	8.6	2.4	39.6	31.7	0.8
T3	30.0	−	−	8.6	21.8	−	−	2.1
T4	−	23.3	−	6.9	−	23.1	−	4.8
T5	−	−	0.4	7.8	−	−	38.4	6.0
注：δk为理论结果相对于第k次（k=1、2、3）试验结果的误差，δs为理论结果相对于模拟结果的误差。

下载: 导出CSV 

| 显示表格

同样的，该算例数值模拟采用Xeon 8358 CPU，32线程并行，耗时3小时34分；理论计算采用相同硬件配置，单线程计算，耗时9 min（计算域较反射算例小1/4），表明该方法在计算时间上具有较明显的优势。

3.   结　论

本文通过精细数值模拟数据拟合的方法，建立了球状TNT装药空气自由场爆炸冲击波压力、密度、粒子速度与时间、距离相关的经验模型，并利用空爆地面反射和空爆建筑后绕射两个典型工况，阐述了提出模型的场景，对提出模型和相关理论方法进行了验证，主要结论如下：

(1) 提出了改进修正Friedlander方程，用以建模冲击波压力、密度、粒子速度等参数随时间变化关系，该模型考虑了爆轰产物对冲击波参数时程曲线形状的影响，近场处，压力、密度、粒子速度时程模型可保证与数值模拟结果的峰值、冲量峰值一致，中远场处，压力、密度、粒子速度时程曲线与数值模拟结果基本吻合，在次级冲击波过渡处略有不同；

(2) 空爆地面反射和空爆建筑后绕射两个工况下，数值模拟压力云图与理论计算压力云图基本吻合，在相同硬件条件下，理论计算耗时仅为千万级网格数值模拟的5%左右，在计算速度上有明显的优越性；

(3) 各测点超压峰值、正相冲量峰值的理论计算与数值模拟结果误差均在10%以内，与实验测试结果平均误差在30%以内，后续可进一步研究多冲击波叠加传播规律、异型装药冲击波传播规律、真实地面反射规律等，提高计算结果准确性。