空中强爆炸冲击波地面反射规律数值模拟研究*

姚成宝，王宏亮，浦锡锋，寿列枫，王智环

（1. 北京大学数学科学学院，北京 100871；2. 西北核技术研究院，陕西 西安 710024）

炸药在距地面有限高度内爆炸时，冲击波将会在地面、建筑等障碍物的表面发生反射和绕射现象，传播规律变得更加复杂。相比于化学爆炸，强爆炸产生的初始压力可达上万吉帕，初始温度高达数十万摄氏度[1-2]，空气在冲击波作用下将发生剧烈的压缩和温升，热力学性质发生显著变化，需要建立高温高压下真实气体的状态方程。此外，由于强爆炸冲击波的传播距离远，其在传播过程中会受到空气不均匀性的影响，空气的初始压强和密度越低，冲击波的衰减越明显[2-3]。空中强爆炸问题具有的上述特点，对数值计算的精度和效率提出了很高的要求。

空中强爆炸问题属于可压缩多介质大变形问题的范畴，一般会存在多种介质的相互作用。由于不同介质的状态方程和初始条件通常存在较大的差异，会给数值计算特别是介质界面附近的数值模拟带来很大的困难。目前针对空中爆炸冲击波问题的主流数值方法是基于Euler 坐标系下的多介质流体数值方法，其包含两个关键步骤：界面追踪和介质间相互作用的处理。其中，常见的界面追踪方法包括VOF (volume-of-fluid)方法[4]、Level Set 方法[5]和Front Tracking 方法[6]等。当确定介质界面的位置后，需进一步准确描述不同介质间的相互作用，常见的处理方法包括虚拟流体类方法(ghost fluid method)[7-8]和切割单元法(cut cell method)[9]。

在爆炸冲击波地面反射规律研究方面，Crowl[10]给出了TNT 等化学炸药在不同爆高下的地面反射冲击波载荷分布规律，乔登江[1]和Glasstone 等[2]依据实验结果分别给出了标准大气状态下1 kt TNT 当量的空中强爆炸冲击波在大尺度范围内的等超压曲线。目前公开发表的文献大多针对TNT 等化学爆炸产生的冲击波在局部范围内的传播[11-15]，较少见到能够完整给出空中强爆炸冲击波在地面反射过程的大尺度范围数值模拟。

针对空中强爆炸问题的强对流性，本文在前期工作[16]的基础上进一步改进基于欧拉坐标系的多介质流体数值方法，并考虑高温、高压下的真实气体状态方程和空气随高度不均匀分布的影响[17]，结合网格自适应技术，计算1 kt TNT 当量空中强爆炸产生的冲击波在不同爆高下的地面反射完整过程，得到地面上距爆心投影点5 km 范围内的冲击波载荷分布。通过和实测结果进行比对，验证计算结果的正确性，并在此基础上研究地面反射冲击波载荷随爆高的变化规律。

1 物理模型和数值方法

1.1 物理模型

根据空中强爆炸问题的对称性，采用二维柱对称计算模型，控制方程组为：

式中： ρ 为密度，u、v分别为r,z方向的速度，E为体积总能量，p为压力。

由于强爆炸伴随着强烈的放热、电离等效应，爆炸产物迅速膨胀至气体状态，且周围空气在受到冲击波、热辐射等因素的综合作用下发生剧烈的加热和压缩，物理状态极其复杂。本文中采用基于Saha 电离平衡理论给出的、以数值表形式给出的状态方程[1]来描述强爆炸产物的热力学行为，并采用真实气体状态方程来描述空气：

式中： γ 为绝热指数，e为比内能。Symbalisty 等[18]给出了8 组参考密度下(10-6～10 kg/m3)的温度与内能、温度与绝热指数之间关系的实测结果，如图1 所示。通过对参考曲线进行双对数插值来获得各热力学变量之间的关系：首先，根据已知的比内能e，对图1(a)进行对数插值，得到相应的温度T；然后，对图1(b)进行插值，获得相应的绝热指数 γ -1 ；最后，将 γ -1 代入式(3)进行求解，最终完成整个状态方程的完整计算。

图 1 不同参考密度下真实气体状态方程的热力学参数关系Fig. 1 Plots of equations of state for real gas at different reference densities

实际空气的初始压力和密度随着高度的增加逐渐降低，且满足下列关系[1,17]：

式 中：Tref=288.16 K，pref=1.01325×105Pa，ρref=1.225 kg/m3分 别 为 标 准 大 气 温 度、压 强 和 密 度;µ=28.966 g/mol，g=9.8 kg/m2，R=8.31 m2·s2/K；z为垂直高度，b为拟合系数。

1.2 数值方法

采用基于欧拉坐标系、具有互不相溶特征的二维柱对称多介质流体计算模型来进行数值离散。在任意时刻，整个计算区域 Ω 可以分成2 个子区域(t) 和(t) ，且满足：

式中： Γ (t) 为2 种互不相溶流体之间的物质界面。

在前期工作[16]的基础上，采用Level Set 方法捕捉爆炸产物与空气的物质界面，根据Level Set 函数分片线性的特征清晰重构出具有分片线性的物质界面，并通过在界面的法向上求解多介质Riemann 问题的精确解来构造爆炸产物与空气之间的数值通量。整个计算流程如下。

1.2.1 物质界面的追踪和重构

Level Set 方法[5,19-20]把随时间运动的物质界面 Γ (t) 定义为距离函数 φ (x,t) 的零等值面。利用特征线方法求解Level Set 函数的演化方程：

式中：u˜ 为由物质界面的运动速度延拓得到的速度场[21]。

根据特征线 方 法 的思想，只 要 知道tn+1时 刻 网 格顶点Xi上的流体在tn时 刻的位置x(tn+1) ，就可 以 得 到新的Level Set 函数在Xi上的值 φ (x,tn+1) 。对每个顶点Xi，利用特征线方法可得：

再利用代数平均，可得新的Level Set 函数：

利用显式正系数格式[21]求解重新初始化方程：

以保持Level Set 函数的符号距离性质，其中， φ0为重新初始化前φ (x,t) 的值，ε 为一很小的正数（例如10-6）。显示正系数格式的求解过程[21]比较繁锁，本文不再说细展开。

当得到tn时刻的Level Set 函数 φn后，根据其分片线性的特征重构出界面单元内分片线性流形的物质界面，如图2 所示，图中nKi,n为物质界面的单位法向。假设单元Ki,n与单元Kj,n具有共边Sij,n，且 包含物质界面ΓKi,n。ΓKi,n将Ki,n和Sij,n依次分为2 个部分Ki±,n和Δtn，此时该单元包含2 种流体，且每种流体的守恒量分别定义为：

图 2 多介质流体计算模型示意图Fig. 2 Schematic diagram for the multi-media fluid calculation model

1.2.2 单元边界数值通量

单元边界的数值通量是指单元边界上同种流体之间的数值通量，且满足：

本文中采用local Lax-Friedrich 数值通量函数：

式中： λ =max{λKi,n,λKj,n} ， λKi,n、 λKj,n分别为和的最大信号速度。

1.2.3 物质界面数值通量

物质界面数值通量是物质界面两侧不同流体之间的数值通量，且满足：

式中： | ΓKi,n| 为 物质界面的长度；和分别为物质界面上的压力和法向速度，且通过在物质界面法向上求解局部一维多介质Riemann 问题的精确解来获得，具体求解过程见文献[22]。

1.2.4 守恒量的更新

当计算得到Ki,n的单元边界数值通量和物质界面数值通量后，可对该单元中每种流体的守恒量进行更新：

2 空中强爆炸冲击波地面反射

2.1 空中强爆炸冲击波反射数值模拟

利用数值模拟程序对当量为1 kt TNT、爆高为100 m 的空中强爆炸冲击波传播过程进行数值模拟。其中，强爆炸产物采用等温等压球模型，取爆炸总当量的85%作为爆炸产物的力学初始能量[1]，产物质量设为50 kg，初始半径设为30 cm。空气采用真实气体状态方程，取水平面处的初始压力为101.3 kPa，初始密度为1.29 kg/m3，并按照式(4)和(5)对不同高度处空气的初始压力和密度进行设置。

计算区域在r=0 和z=0 设为固壁反射边界条件，其余边界设为无反射边界条件。数值计算中采用h-网格自适应技术[23]，根据相邻网格的压力梯度来进行网格的局部加密或稀疏，得到典型时刻的冲击波压力等值线以及相应的网格自适应图，如图3 所示。

由图3 可知，当空中强爆炸产生的冲击波传播到地面时，最早在爆心投影点发生正反射，然后以逐渐增大的入射角发生斜反射，产生双激波结构的规则反射。随着冲击波向外继续传播，入射角不断增大。当入射角增大到临界角附近时将发生马赫反射，此时反射波阵面和入射波阵面的交点离开地面，形成垂直于地面的马赫杆。随着冲击波的进一步传播，马赫杆迅速增长。图4 给出了数值计算得到的地面冲击波载荷分布(峰值超压和冲量)与实测结果[2]的对比情况，其中在距爆心投影点1 km 以内，数值计算得到的峰值超压与实测结果符合较好，而在1 km 以外则比实测结果略低，这是由于随着计算区域的尺度增大，网格自适应后的单元尺寸不够精细，使得数值计算对峰值超压的捕捉能力减弱导致。此外，数值计算得到的冲量比实测结果略大，但最大误差不超过50%。

图 3 典型时刻的冲击波压力等值线图和网格自适应图Fig. 3 Pressure contours and adaptive meshes at typical times

图 4 地面不同距离处的冲击波峰值超压和冲量Fig. 4 Peak overpressures and impulses at different radii

2.2 不同爆炸高度的冲击波反射规律

利用数值模拟程序进一步对当量为1 kt TNT，爆高H依次为20、50、100、200、400 和800 m 的空中强爆炸冲击波传播过程进行了数值模拟，得到了地面上不同距离处的冲击波峰值超压和冲量，如图5 所示。

图 5 不同爆高下的地面冲击波峰值超压和冲量Fig. 5 Peak overpressures and impulses at different heights of burst

图5 的计算结果表明，对1 kt TNT 当量的空中强爆炸，在距爆心投影点300 m 范围以内，地面上的峰值超压和冲量随着爆高的增大而逐渐减小。在距爆心投影点300 m 范围以外，随着距爆心投影点的距离的增大，不同爆高条件下地面反射冲击波载荷的差距逐渐变小。其中，爆高约为200 m 的空中强爆炸产生的地面峰值超压和冲量处于最大值，其他爆高条件下的冲击波载荷随着爆高的进一步增大或减小均呈现逐渐减弱的趋势，其原因主要是由于空中强爆炸地面反射而引发的复杂冲击波传播引起的。当空中强爆炸形成的冲击波达到地面时，随着入射角的增大，会依次在地面形成规则反射和马赫反射等复杂的波系结构。其中，在规则反射区，反射冲击波载荷由入射波强度和入射角决定。对于距爆心投影点相同距离的观察点，爆高越大，冲击波载荷越弱。

当入射角增大到临界角附近时，冲击波将会发生马赫反射。马赫反射的特点是，反射波阵面和入射波阵面的交点离开地面，形成垂直于地面的马赫杆。马赫反射对地面冲击波载荷分布起到了双重效果，一方面马赫反射的出现使得冲击波载荷的强度局部增大，另一方面由于马赫杆高度的不断增大，加快了冲击波波阵面能量的扩散速度，使得马赫反射区地面的冲击波载荷衰减速度加快。因此，空中强爆炸冲击波在地面反射后产生的大尺度空间载荷分布与爆炸当量以及爆炸高度密切相关，在实际应用中需要根据预定的目的，正确地选择爆炸方式，可以有效地扩大冲击波破坏区域。

3 结 论

本文中建立了基于欧拉坐标系的互不相溶、具有清晰锐利界面的可压缩多介质流体数值方法，并研制了适用于模拟空中强爆炸冲击波传播的数值计算程序。利用Level Set 方法捕捉爆炸产物与空气的运动界面，并清晰重构出分片线性流形的物质界面。采用有限体积方法求解每种物质的控制方程组，并通过在物质界面法向求解局部一维多介质Riemann 问题的精确解来计算物质界面的数值通量。数值计算中采用了网格自适应技术，在捕捉激波峰值的前提下，有效提高了计算效率。

利用计算程序首先对空中强爆炸冲击波遇到刚性地面时的反射过程进行了数值模拟，得到了与实测数据一致的数值结果，验证了数值模拟程序的可靠性。在此基础上进一步对不同爆高条件下强爆炸冲击波的反射过程进行了数值模拟，给出了不同爆高条件下地面反射冲击波峰值超压和冲量在不同反射区域的分布情况，得到了冲击波载荷随爆高的变化规律，并基于数值模拟结果给出了强爆炸冲击波对典型目标结构毁伤破坏的最佳爆炸高度，为空中强爆炸冲击波的杀伤破坏效果预测以及爆炸方式的选取提供了研究基础。