爆炸冲击波是爆炸对周围目标产生破坏效应的重要因素之一，自20世纪40年代以来就受到广泛关注^[1]。爆炸冲击波在空气中传播，遇到地面或其他结构时会发生反射。为了准确评估破坏效应，需要确定冲击波反射后的流场分布。本文中拟建立一种理论方法，用于对任意当量、爆高条件下冲击波在固壁反射后的流场物理量进行快速计算。

常用的快速计算方法，大致分为2类。一类是基于数值或实验数据建立物理量的拟合关系式^[2-5]，该方法的精度取决于数据规模、自变量数目和关系式类型等，所得到的关系式通常具有一定的适用范围和场景。另一类方法是将固壁假想为对称面，在固壁另一侧建立与真实爆炸对称的虚拟爆炸，然后根据一维球爆炸自由场参数确定真实和虚拟爆炸流场，再将两者进行线性或非线性叠加，从而得到冲击波反射后的流场分布，这种方法称为镜像方法^[6-12]。与第一类方法相比，镜像方法更易于建立和使用，可以给出较清晰的物理图像，且不受限于特定的爆炸场景。

爆炸冲击波在固壁可能发生正规或非正规反射，对应的流场分布及演化是不相同的。在镜像方法中，目前已有的工作大都假定冲击波始终保持球面，即仅考虑了正规反射的情形，对于非正规反射的情形则没有涉及。另外，压力和动压均是评估冲击波毁伤的重要参量，其中流场压力可采用非线性叠加模型LAMB^[12]计算，而对于密度和速度则采用线性叠加计算，所得到的结果与真实值相差较大^[2]。

为了能够将镜像方法推广到任意反射类型时爆炸流场的计算，本文中建立正规和非正规反射阶段流场中冲击波传播的计算方法，并基于LAMB模型发展新的叠加模型LAMBR (LAMB revisied)。针对任意爆炸当量、爆高条件，首先确定任一时刻流场中的冲击波结构，然后利用镜像方法和LAMBR叠加模型，快速计算流场中压力、密度和速度等物理量分布。通过与数值模拟结果和美国结构抗偶然爆炸手册(UFC 3-340-02，后文简称为UFC)^[13]中的数据比对，验证该理论方法的合理性。

1.   理论方法

1.1   冲击波传播

当爆炸冲击波传播到达固壁时，首先发生正规反射。在t=t₀时刻，入射波阵面与固壁夹角超过某一临界值，正规反射转换为非正规反射。常用的转换准则，包括脱体准则、声速准则和力学平衡准则，考虑到爆炸冲击波在固壁反射属于非定常流动，本文中选取声速准则^[14]计算非正规反射的起始位置T₀。在t＞t₀时刻，三波点T逐渐向着远离固壁的方向运动，见图1。图1中O为爆炸中心，H为爆高，固壁与y=0重合，RR、IR分别表示正规反射和非正规反射。

图  1  爆炸冲击波反射示意图

Figure  1.  Schematic diagram of blast wave reflection

下载: 全尺寸图片 幻灯片

对于正规反射和非正规反射，均假定反射波是以点(0,–H)为中心的圆弧，根据不同时刻入射波传播距离，确定反射波形状。对于非正规反射，假定马赫杆是以爆心在固壁投影点为中心的圆弧，通过求解三波点T的位置，确定马赫杆形状。已有的计算点T位置的方法，可以分为基于几何相似^[15]和基于流体动力学^[16-19]两类。

1.1.1   基于几何相似的方法

图1中，三波点由点T₀开始，随时间发展逐渐远离壁面，在t₁时刻到达点T。此时，入射波阵面延长线、马赫杆分别与固壁交于点B、F。点A为直线OB与t₀时刻冲击波阵面的交点，直线AC与OT₀平行，CE与OB垂直，CD、EF与固壁垂直，则△BFE与△BDC是相似的。马赫杆与固壁交点F的坐标，可以表示为：

式中：R=R(t)为 t 时刻入射波半径，θ₀、θ_B 分别为入射波阵面在点T₀、B处的传播方向与水平方向的夹角(以逆时针方向为正)，g为与爆炸条件相关的常数。三波点T的坐标为：

易仰贤^[15]假定常数g仅与比例爆高$\bar H = H/{W^{1/3}}$有关(W为炸药质量)，并给出了核爆炸情况下的关系式g= g($\bar H$)，但是该关系式并不能直接应用于化爆情况。因此，本文中根据UFC^[13]中提供的三波点坐标数据，按照上述过程反推得到适用于化爆情况的g。定义绝对偏差：

式中：${\bar y_{T{\rm{,Yi}}}}{\rm{ = }}{y_{T{\rm{,Yi}}}}/{W^{1/3}}$为基于给定的g值和上述几何相似方法计算得到的三波点的比例纵坐标，${\bar y_{T{\rm{,UFC}}}}{\rm{ = }}{y_{T{\rm{,UFC}}}}/{W^{1/3}}$为UFC中提供的三波点的比例纵坐标，${\bar x_T} = {x_T}/{W^{1/3}}$为三波点的比例横坐标，根据UFC中数据确定其取值范围为${\bar x_{T{\rm{,1}}}} \text{≤} {\bar x_T} \text{≤}{\bar x_{T{\rm{,2}}}}$。最佳的g值，应满足“对于给定比例爆高$\bar H$，ε 取得最小值”。表1中列出了本文中按照这一准则所确定的g值，记为${g_{{\rm{UFC}}}}$；同时，表中还列出了易仰贤^[15]所给出的g值，记为${g_{{\rm{Yi}}}}$。对比${g_{{\rm{UFC}}}}$和${g_{{\rm{Yi}}}}$，发现两者具有明显的差别。

表  1  不同比例爆高条件下g 的取值

Table  1.  The value of g at different scaled heights of explosion

$\bar H$/(m·kg−1/3)	${\bar x_{{{T,1}}}}$/(m·kg−1/3)	${\bar x_{{{T,2}}}}$/(m·kg−1/3)	gUFC	gYi
0.396 7	0.793 4	1.895 8	0.837 7	0.361 1
0.595 1	0.793 4	2.444 0	0.774 2	0.299 5
0.793 4	0.793 4	3.424 8	0.638 8	0.254 0
0.991 7	1.586 8	4.020 8	0.537 2	0.218 9
注: UFC 中数据采用英制单位，本文使用时换算为国际单位制。

下载: 导出CSV 

| 显示表格

需要说明的是，由于UFC中没有提供非正规反射起始点T₀的坐标信息，因此本文中基于声速准则和UFC中爆炸自由场数据计算点T₀的位置。

1.1.2   基于流体动力学的方法

在任意时刻 t (t＞t₀)，如果已知三波点T的坐标 (${x_T}$,${y_T}$)和该时刻三波点迹线与水平方向的夹角χ，则有：

因此，给定三波点T起始位置T₀，按时间逐步推进，即求得三波点轨迹。

将三波点附近流动近似看作拟定常流动，可以采用von Neumann的三波理论^[16]求解χ。点T的运动速度D_T与χ之间满足：

式中：D_I为入射波传播速度。建立与点T固连的坐标系${\Re _T}$，并在坐标系${\Re _T}$中观察流场运动，则穿过入射波后流动方向偏转角度δ_I满足：

式中：ξ为波后与波前流场压力之比，Ma₀为波前流场马赫数。穿过反射波和马赫杆后流动方向偏转角度，也按式(6)计算。根据界面两侧连续条件，联立冲击波关系式，建立迭代求解过程，即可得到χ。

也可基于几何激波动力学(geometrical shock dynamics, GSD)^[20]理论求解χ^[17-19]。在三波点两侧，入射波与马赫杆的传播角度θ_I和θ_M满足：

式中：A_K = f (Ma_K)，下标K=I，M分别对应入射波和马赫杆。通过求解上式得到Ma_M，进而有：

在式(7)和(8)中，关系式A_K=f (Ma_K) 可参照Itoh等^[21]文中式(26)确定。

为了记述方便，将基于几何近似和表1中g_UFC值的方法记为Yi-UFC方法，将基于von Neumann三波理论的流体动力学方法记为VN方法，将基于GSD理论的流体动力学方法记为GSD方法。通过与数值结果和UFC数据比对，从中选取合理的方法。

1.2   叠加模型LAMBR

入射波后、反射波前的流场，与真实爆炸自由场无异。反射波后的流场，则利用镜像方法计算，见图2，将真实爆炸和虚拟爆炸自由场叠加，即可给出反射波后的流场分布。

图  2  镜像法示意图

Figure  2.  Diagram of the method of image for blast wave reflection

下载: 全尺寸图片 幻灯片

当冲击波在固壁发生正规反射时，见图2(a)，对于反射波后某一点P，记真实爆炸自由场在该点处的物理量为p₁、ρ₁、u₁，虚拟爆炸自由场在该点处的物理量为p₂、ρ₂、u₂，流体满足完全气体状态方程 p = (γ−1)ρe，γ为比热比。则点P处物理量按下述模型计算。

1.2.1   压力

按LAMB模型计算流场压力，有：

式中：p₀、ρ₀分别为入射波前未扰动流场的压力、密度。

对于冲击波在固壁发生正反射的情况，按照式（9）计算得到的压力与精确解完全一致。当波后流场密度较小时，可能得到$\bar \rho$＜0，进而给出数值极大、不符合物理的p。本文中建议选取$\bar \rho$=max{$\bar \rho$, min{ρ₁, ρ₂}}。

1.2.2   密度

跨越冲击波是一个非线性且熵增的过程，因此波后密度ρ大于线性叠加得到的密度$\bar \rho$，同时又小于流体经历等熵过程后的密度。因此，根据p与p₁的相对大小计算波后密度ρ：

1.2.3   速度

根据动量守恒计算速度预估值：

考虑到式(12)本质上是线性叠加，因此仅适用于压力较小时的情况，需要引入一定的修正。本文中利用激波关系式得到固壁处冲击波阵面后的速度，然后在保证该速度幅值不变的前提下，对式(12)给出的流场速度的幅值进行等比例修正，给出流场速度u。

将该模型记为LAMBR (LAMB revised),其与LAMB模型的主要区别在于密度和速度的计算。当冲击波在固壁发生非正规反射时，见图2(b)，仍然采用LAMBR模型计算点P处的物理量。与正规反射不同，此时反射波R超越入射波I传播，即反射波半径大于入射波半径。在确定点P处的爆炸自由场物理量时，应先根据当前时刻反射波半径与入射波半径之比，对点P到爆心距离进行等比例缩放，再根据修正后的爆心距确定自由场物理量。

2.   结果与讨论

利用本文中建立的理论方法，计算1 kg TNT爆炸冲击波在固壁反射后的流场，并与UFC数据和动力学软件AUTODYN计算结果进行对比，检验理论方法的合理性。初始时刻空气压力p₀=101.325 kPa，密度ρ₀=1.225 kg/m³，比热比γ=1.4，且处于静止状态。炸药TNT选用JWL状态方程描述:

式中: V=ρ₀/ρ为相对体积，e为质量内能，A、B、R₁、R₂和w为常数。TNT的JWL 状态方程参数^[22]：初始密度ρ₀=1 630 kg/m³, A=371.2 GPa, B=3.231 GPa, R₁=4.15, R₂=0.95, w=0.30, 爆速D=6.930 km/s, 初始体积能量E₀=7 GPa, CJ压力p_CJ=21 GPa。当炸药体积膨胀到约10倍初始体积时，为了减小计算误差，AUTODYN会改用理想气体状态方程来描述炸药产物的力学行为，取比热比γ_TNT=1.3。

2.1   爆炸自由场

在AUTODYN中建立一维轴对称楔形模型（r_min=1.00×10⁻³ m, r_max=6.00 m），利用Euler-Multi-Material 求解器，计算1 kg TNT在空气中的爆炸自由场，炸药半径r=5.27×10⁻² m。采用2种网格划分方案，网格尺寸dr分别取为1.00×10⁻⁴ m和5.00×10⁻⁵ m。图3中给出了爆炸自由场冲击波峰值超压Δp_s、走时t_a随爆心距 r变化的计算结果，并列出了UFC中的数据。利用不同网格方案得到的结果基本重合，验证了数值计算的网格收敛性。在爆心距 0.3 m ≤ r ≤ 3.0 m范围内，计算结果与 UFC 中数据基本吻合，验证了数值计算的合理性。在后文的理论分析和数值模拟中，均选取dr = 5.00×10⁻⁵ m的计算结果作为爆炸自由场数据。

图  3  1 kg TNT 爆炸自由场冲击波参数

Figure  3.  Free field blast wave parameters for 1 kg TNT explosion

下载: 全尺寸图片 幻灯片

2.2   爆炸冲击波在固壁反射

选取4种不同的爆高条件，见表1，利用理论方法对爆炸流场进行计算。同时，在AUTODYN中建立二维轴对称模型，利用Euler-Multi-Material 求解器进行数值计算。计算区域为0≤x≤4H, 0≤y≤3H，以y轴为对称轴，在边界x=4H, y=3H处施加 outflow 边界条件，在边界y=0处施加固壁边界条件，爆心位于(0, H)处。采用2种网格划分方案，网格尺寸分别取dx=dy=H/500 和 dx=dy=H/667。初始时刻流场，由t^* 时刻爆炸自由场映射得到，要求冲击波传播距离R ＜H。4种爆高条件对应的t^*，分别为 9.60×10⁻⁵、1.88×10⁻⁴、3.16×10⁻⁴和 4.80×10⁻⁴ s 。

2.2.1   三波点轨迹

图4给出了不同爆高条件下三波点轨迹的理论和数值结果，并列出了UFC中数据。其中，dx=dy=H/500, dx=dy=H/667分别为相应网格划分方案下的数值结果，通过读取不同时刻压力云图中三波点位置得到。

图  4  不同爆高情况下三波点轨迹

Figure  4.  Triple point trajectories at different heights of explosion

下载: 全尺寸图片 幻灯片

对比不同网格方案得到的结果，两者基本重合，验证了数值计算的网格收敛性。对比数值结果和UFC数据，发现在不同爆高H条件下，数值结果与UFC数据在y_T＜0.4～0.5 m范围内存在一定程度的偏差，可能是由于数值模型对爆炸产物或物理过程的描述与真实情况之间还存在差异，也可能是由于UFC数据本身的误差。随着三波点远离固壁或是爆高增加，两者趋于吻合，表明数值结果是合理的。

将理论结果与UFC数据、数值结果进行对比。首先，对比非正规反射起始点位置，见表2，理论与数值计算结果基本一致，表明选取声速准则作为正规与非正规反射的转换准则是合理的。而UFC中缺少该数据，因此未参与比较。接下来，对比三波点轨迹。VN与GSD同属流体动力学方法，仅三波点附近流场求解方法不同，因此当H=0.396 7 m和0.595 1 m时，两者的结果基本一致，但与UFC和数值结果差异较大；当H=0.793 4 m和0.991 7 m时，两者的结果在靠近非正规反射起点附近基本一致，随着三波点远离壁面，VN结果略大于GSD。与VN和GSD相比较，对于不同爆高条件，Yi-UFC方法的结果与UFC数据和数值结果的吻合程度均优于VN和GSD方法，表明Yi-UFC方法能够更好地描述三波点的运动。因此，选取Yi-UFC方法用于流场冲击波结构的理论计算。

表  2  非正规反射起始点坐标${x_{{T_0}}}$

Table  2.  The ${x_{{T_0}}}$ value of the starting point of IR

H /m	${x_{{T_0}}}$/m
	数值计算	理论分析
0.396 7	0.333 2	0.328 4
0.595 1	0.473 4	0.488 4
0.793 4	0.620 2	0.646 4
0.991 7	0.796 0	0.805 9

下载: 导出CSV 

| 显示表格

2.2.2   不同时刻流场分布

以爆高H=0.793 4 m为例，图5、6分别给出了不同时刻流场压力云图以及沿直线y=0.03 m的物理量分布。其中，数值计算选用dx=dy=H/667的网格划分方案。

图  5  H=0.793 4 m时，不同时刻流场压力云图(黑线为物质界面)

Figure  5.  Pressure contours at various time instants (black line stands for material interface) for H=0.793 4 m

下载: 全尺寸图片 幻灯片

在t=4.76×10⁻⁴ s 时刻，入射波在固壁发生正规反射。理论方法给出的流场压力云图与数值结果较接近，见图5(a)，能够反映流场中的波系结构。在真实情况中入射波后流场是非均匀的，反射波的形状势必不能保持圆弧；同时，反射波会在物质界面处发生折射，形成新的反射波和透射波，物质界面的运动速度也逐渐减小。理论方法近似假定反射波为圆弧，且物质界面按照原一维球对称传播规律运动，这与真实情况有差异，因此理论结果中反射波和物质界面的形状与数值结果不完全一致。直线y=0.03 m依次与反射波、入射波相交，沿该直线的物理量分布出现2个峰值，见图6(a)，基于同样的原因，理论结果中反射波位置与数值结果相差约为0.03 m。除此之外，压力、速度的理论与数值结果基本重合。

图  6  H=0.793 4 m时，不同时刻流场物理量沿y=0.03 m的分布

Figure  6.  Parameters along y=0.03 m at various time instants for H=0.793 4 m

下载: 全尺寸图片 幻灯片

在t=9.96×10⁻⁴ s时刻，入射波在固壁发生非正规反射。对比流场压力云图，理论与数值结果的差异主要集中在对称轴附近区域，而在其他区域则较接近，见图5(b)。对比沿直线 y=0.03 m的物理量分布，理论与数值结果中的马赫杆位置基本重合，表明理论方法中将马赫杆近似为圆弧是合理的，与真实情况较符合，这一点在后续时刻的比较中也可以得到说明。Hu等^[16]曾指出，当爆高小于某一临界值时，入射波在固壁的非正规反射会依次经历双马赫反射、过渡马赫反射和单马赫反射等3个阶段。根据数值结果判断，此时流场中的反射为过渡马赫反射，沿直线 y=0.03 m的物理量分布出现了2个峰值，见图6(b)。基于反射波的圆弧假定，理论方法仅能构建单马赫反射类型，所以在过渡马赫反射的影响区域内，理论结果中并未出现第2个峰值。

在t=1.92×10⁻³ s时刻，非正规反射类型转换为单马赫反射，此时流场波系结构与理论方法中的假定是一致的，因此流场压力分布的理论与数值结果较吻合，见图5(c)，沿直线y=0.03 m物理量分布的理论和数值结果也趋于吻合，见图6(c)。此时，在爆心附近出现次级冲击波并向外传播，形成高压区域。受到反射波的影响，数值结果中该高压区域的位置相较于理论结果更靠近爆心上方。考虑到流场发展本身具有极强的非线性，而本文中所用的叠加方法并不能完全精确地求解流场中的非线性作用，这就导致在非线性作用显著的区域，诸如靠近爆心和对称轴附近的区域，理论与数值结果存在一定的偏差。

在t=3.60×10⁻³ s时刻，流场中的次级冲击波向外传播，并在固壁发生反射，同时次级冲击波后流场压力逐渐减小。理论方法给出的流场压力云图能够反映这一特征，但是爆心附近的压力小于数值结果。对比沿直线y=0.03 m物理量分布，见图6(d)，在次级冲击波反射点之外的区域，即x＞1.47 m时，理论结果与数值结果基本吻合；而在x＜1.47 m的区域中，理论结果则小于数值结果，次级冲击波的位置相差约0.46 m。

通过以上对比可以发现，尽管理论方法对于流场发展过程进行了一定简化，但无论流场发展是处于正规反射阶段还是非正规反射阶段，理论方法所得到的结果均合理地反映了流场中波系发展特征，尤其是在靠近冲击波的区域，不同时刻的物理量分布规律与数值结果基本吻合，物理量峰值与数值结果的偏差大都小于10%，见表3，从而验证了理论方法的合理性。表3中同时列出了基于LAMBR和LAMB模型得到的理论结果，两者所给出的压力峰值相同，但对于密度和速度峰值，LAMBR模型的结果明显优于LAMB模型。

表  3  沿直线y=0.03 m的物理量峰值

Table  3.  Peak values of parameters along the line y=0.03 m

t/s	ps/(105 Pa)				ρs/(kg·m−3)						us/(m·s−1)
	NS	TA	ε/%		NS	TA1	ε1/%	TA2	ε2/%		NS	TA1	ε1/%	TA2	ε2/%
4.76×10−4	40.82	36.56	−10.44		7.81	7.86	0.64	6.54	−16.26		728.93	688.84	−5.50	513.35	−25.45
9.96×10−4	9.73	10.70	9.97		4.42	5.70	28.96	5.71	29.19		732.33	766.03	4.60	482.88	−34.06
1.92×10−3	4.46	4.41	−1.12		3.19	3.37	5.64	3.76	17.87		418.37	413.92	−1.06	342.21	−18.20
3.60×10−3	2.51	2.57	2.39		2.29	2.36	3.06	2.59	13.10		239.09	245.65	2.74	229.38	−4.07
注：NS 为数值结果，TA 为理论结果，ε=(TA-NS)/NS 为偏差。下标 1、2 分别表示基于 LAMBR 和 LAMB 模型的理论结果。

下载: 导出CSV 

| 显示表格

在运行时间上，以截止时刻t=3.60×10⁻³ s为例，采用相同的计算区域和网格划分方案，数值计算耗时约70 h，理论方法仅需约15 s，远远快于数值计算。

3.   结　论

爆炸冲击波遇到固壁发生反射，使得波后流场改变。本文中基于镜像法，在固壁另一侧建立与真实爆炸对称的虚拟爆炸，将爆炸冲击波在固壁的反射等效为真实与虚拟爆炸流场的相互作用，建立了适用于爆炸冲击波反射后流场的快速计算方法。

该方法首先将反射波、马赫杆分别简化为以虚拟爆源和地面爆心投影点为中心的圆弧，选取声速准则作为正规与非正规反射类型的转换准则，利用基于几何近似的Yi-UFC方法，计算不同时刻流场中冲击波结构。然后，根据冲击波结构，将流场划分为不同的区域。对于位于反射波前、入射波后的区域，直接根据爆炸自由场数据确定流场分布。对于反射波后的流场，则通过本文新发展的流场叠加模型LAMBR进行计算。通过与数值结果和已有数据进行对比，理论方法给出的压力、密度和速度等物理量分布规律、峰值等与数值结果基本吻合，能够识别流场中的次级冲击波传播等物理过程，验证了理论方法的合理性。与LAMB模型相比，LAMBR模型得到的密度和速度等物理量的峰值与数值结果更为吻合。同时，该理论方法所需要的时间，较数值计算大幅度减少。

需要注意的是，流场发展本身是非线性极强的物理过程，鉴于理论方法所采用的诸多假定，在流场非线性作用显著的区域中，理论结果与真实情况还存在一定差异。