Chinese Journal of Ship Research
Development of in-house high-resolution hydrocode for assessment of blast waves and its application
XU Weizheng1，2，WU Weiguo1，2 1 Key Laboratory of High Performance Ship Technology of Ministry of Education，Wuhan University of Technology， Wuhan 430063，China 2 School of Transportation，Wuhan University of Technology，Wuhan 430063，China Abstract：The propagation and evolution characteristics of blast waves in confined spaces are complicated due to the constraint of the surrounding walls, by which the enhanced reflected shock waves will cause more serious damage to the internal structures, facilities and personnel. In order to investigate the characteristics of explosions in confined spaces, an in-house high-resolution hydrocode was developed in this present work. The third-order WENO finite difference scheme (weighted essentially non-oscillation scheme) was implemented in the code to capture the shock waves generated by cylindrical explosives. The Sod shock tube problem, interacting blast wave problem and blast in air problem were simulated to validate the code. The validated code was then used to simulate the blast waves generated by condensed explosives in closed, vented and connected spaces. The propagation of blast waves and the characteristics of blast load were subsequently investigated. The developed code appears to accurately predict the process of explosions in confined spaces. This high-resolution hydrocode can be used to study the propagation paths of blast waves in complicated spaces and evaluate the internal blast load, which can provide reliable input for the design of explosion-resistant structures. Key words：explosive shock wave；confined space；numerical simulation；high-resolution hydrocode； WENO scheme；program development
0引言
随着现代反舰武器的迅速发展，各种高性能的半穿甲反舰导弹已成为水面舰船水线以上部分舷侧的主要威胁，其对舰攻击破坏特点是穿透舰船舷侧外板在舱室内发生爆炸。当爆炸发生在舱室内时，由于冲击波在传播过程中受到舱室壁面的约束限制，将产生持续时间短暂、峰值逐渐衰弱的冲击波以及持续时间较长的准静态超压［1-2］，并对舱室结构、内部设施及人员造成严重的损伤。近年来，国内外学者针对舱室、约束空间内爆炸载荷，在试验、数值模拟、高精度数值计算方法等方面开展了大量研究工作。Edri等［3］在设置有泄压口的长方体房间内开TNT展了 药柱爆炸实验，分析了装药质量对内爆Wu等［4］开展了封闭空间内炸载荷的影响规律。爆炸实验，研究了炸药形状和起爆位置对爆炸冲击波的影响规律。胡洋等［5］用压力传感器记录了长方体混凝土密闭空间内爆炸壁面上爆炸载荷的压力时间历程曲线，分析了壁面上爆炸载荷的分布规律。侯海量等［6］开展了舱室内爆炸冲击波载荷特性的实验研究。孔祥韶等［7］实验研究了角隅结构形式对舱内爆炸载荷的影响规律。侯海量等［8］ MSC.DYTRAN基于 软件，数值研究了舱室内爆炸MSC.DYTRAN ［9］冲击波的载荷特性。孔祥韶等 基于软件数值研究了舰船舱室内战斗部爆炸及爆炸毁AUTODYN伤效应。丁阳等［10］利用 商用程序中的Remap技术对室内爆炸进行模拟，研究了壁面爆LS-DYNA炸载荷的分布规律。樊壮卿等［11］基于软件，数值模拟了大体积、复杂结构的典型舱室内爆炸载荷传播特性。综上所述，当前对约束空间内部的爆炸载荷2种研究主要采用试验测试和商用软件数值模拟方法。鉴于爆炸试验昂贵且存在一定的风险，数值模拟成为研究爆炸波传播及爆炸载荷特性的主LS-DYNA，MSC.DYTRAN要手段。然而，包括 和AUTODYN在内的现有商用程序的求解器主要采用传统的二阶精度算法来模拟爆炸过程，如有限Euler-FCT Roe差分 和有限体积 方法。这些方法因计算精度较低，会严重抹平爆炸冲击波峰值压力。准确预报爆炸载荷是合理设计和评估抗爆结构的依据，而目前发展成熟的高精度计算方法尚未集成到商用程序中，故研究和开发高精度爆炸波数值计算程序具有重要的工程意义。对于属于高压力比、高密度比的炸药爆炸问题，其数值模拟对激波捕捉格式提出了更高的要 Liu 1994 WENO等［12］于求。为此， 年提出了 格式（Weighted essentially non-oscillation scheme），Shu WENO等［13-15 ］发展了该格式。 格式在有效性、通ENO量的光滑性和收敛解的稳定性方面均优于WENO格式［16］，所以将 格式用于模拟爆炸过程是一个比较好的选择。FORTRAN WENO本文将基于 平台，采用三阶有限差分格式，自主开发舱室内的爆炸波高精度数值计算程序，利用所开发的程序研究密闭空间、泄压空间及连通空间内部的爆炸波传播规律与爆炸载荷特性。该研究工作可为后续爆炸波传播、爆炸载荷以及抗爆结构设计提供基础。
1 程序开发 1.1 欧拉方程
基于瞬时爆轰假定，将炸药等效为高压、高密度气体。爆炸流场的控制方程采用三维可压缩欧拉方程［17］进行描述： ¶U + ¶E + ¶F + ¶G =0 （1） ¶t ¶x ¶y ¶z式中，U 为守恒通量；E ，F ，G 为数值通量，并
分别表示如下：
ρ ρu ρu ρu2 +p U= ρv ， E= ρuv ， ρw ρuw E u(E + p)
ρv ρw ρvu ρwu （2） F= 2 G= ρv +p ， ρwv 2 ρvw ρw +p v(E + p) v(E + p) E = ρe + 1 ρu2 + 1 ρv + 1 ρw2 （3） 2 2 2 2 （4） p= (γ - 1) ρe式中： ρ 为密度； u v w 分别为 x y z 方向上的速度分量；p 为流体压力；E 为单位体积流体的总能量；e为比内能；γ为气体的绝热指数，本文数1.4。值计算中取为
1.2 数值离散方法
Strang采用 维数分裂法，将欧拉方程分解为WENO x y z 方向进行求解。采用三阶 有限差分格式对欧拉方程的空间项进行数值离散，具体离散过程如下。单元面中心点 xi 处的数值通量 fi 为1/2 + 1/2
fi = 1( fi + fi 1) - 1 ω 0(D - D1) （5） + 1/2 2 + 2 0式中：fi ，fi 分别为点 xi ，xi 处的数值通量； + 1 + 1 ω 为非线性权重；D = fi - fi ；D1 = fi - fi 。0 0 - 1 + 1式（5）中的 ω 由下式求得： 0 1α0 ω 0 å αk
s =0 dk ；k=0，1 αk = (ε + βk) d = 1 ，d1 = 2 6 （） 0 3 3 =1.0×10-6；k式中：ε 为避免分母为零的小数，取 ε WEND为三阶 格式子模权的个数； dk 为三阶WENO格式的线性权值； αk 为转换函数，三阶WENO 2格式的 个子模板为{xi 1 xi} ，{xi xi 1} ； - + βk (k = 0 1) 为光滑因子，其表达式如下［18］： βk 2 = D0 7 （ ） 2 β1 = D1 TVD-RK法［19］对欧拉方程的时间项采用三阶进行数值离散，具体离散格式如下： (1) n n) U= U +DtL(U (2) n (1) + 1 DtL(U (1)) 8 = 3 + 1 U U U （） 4 4 4 n+ 1 n (2) + 2 DtL(U (2)) = 1 + 1 U U U 3 3 3 n (1) (2)式中：U 为n时刻的守恒通量；U ，U 为中间n+1变量；U n+ 1 为 时刻的守恒通量；Dt 为时间步长；L(×)表示运算算子。
2 算例验证 2.1 Sod激波管
该算例的初始条件如式（9）所示［20］，两端边界条件设置为流出边界，网格数为200个，计算结束0.18。图 1时间为 所示为计算结束时的无量纲压力曲线。
(1 0 1)T 0 x < 0.5
（9）
( ρ u p)T = (0.1250 0.1) T 0.5 x 1式中：ρ u p分别为无量纲密度、速度和压力。2.2 双爆轰波碰撞该算例的初始条件如式（10）所示［21］，两端边800界条件设置为反射边界，网格数为 个，计算结0.038。图 2束时间为 所示为计算双爆轰波碰撞结束时的无量纲压力曲线。(1 0 1000)T 0 x < 0.1 （10） ( ρ u p)T = (1 0 0.01)T 0.1 x 0.9 (1 0 100)T 0.9 x 1
2.3 空中爆炸
5 000 mm，装 100 mm。设计算域长度 药半径炸药瞬时爆轰等效后的高压、高密度气体参数如=1 630 kg/m3，P=3.057 9×109 Pa。经多次数值下：ρ 10 mm。图3试验，网格尺寸取 所示为选取典型位2，3 m AUTODYN置（即 ）的压力时间历程输出与商用程序输出结果的对比。1、图2、图3由图 的对比结果可知，本文所开发的数值计算程序具有较高的精度和一定的可靠性。
3 密闭空间内爆炸 3.1 初始条件设置
1 800 mm×800 mm×800 mm，设密闭空间尺寸为5个测点（No.1~No.5）对爆炸超压时间历壁面设置 4）。500 g程进行输出（图 柱形炸药放置在密闭空间中心位置，等效后的高压、高密度气体参数如99 mm 80 mm，ρ = 203.75 kg/m3，下：半径 ，高度=3.822×108 Pa。计算初始条件如图5（a P ）所示，考虑到计算时间及精度的要求，经多次数值试验， 14.4 90×40×40 5（b网格数取 万个（ ），如图 ）所示。壁面边界条件设置为刚性反射边界条件，这里不考虑冲击波与结构的耦合作用。
3.2 爆炸波传播过程
6图 所示为不同时刻输出时密闭空间内部的6爆炸初期压力分布云图。由图 可以分析出爆炸初期爆炸波的传播过程：高压气体首先进行三维柱对称自由膨胀，当爆炸波初次到达壁面时发生 6（a））；由于四周壁面约束，爆炸波正规则反射（图向密闭空间长度方向的端面传播并在两壁面交线
处形成局部压力汇聚现象，且局部汇聚压力峰值6（b））；当爆炸波到沿密闭空间长度方向传播（图达密闭空间端面处时，三壁面角隅附近区域形成6（c）），随后两端面反射波以近压力汇聚现象（图 6（d））。似平面波的方式向密闭空间中部传播（图
3.3 爆炸载荷特性
7图 所示为不同时刻输出时密闭空间内部的爆炸壁面测点超压及冲量（I）时间历程曲线。由7（a）可知，不同测点处爆炸前期爆炸波超压峰图值大小有所不同，而准静态超压峰值几乎相同（图7（a）中绿色粗实线），这说明密闭空间内部爆炸形7（b）成的准静态超压峰值在空间上是均匀的。由图 可知，密闭空间内部的爆炸壁面测点冲量时间历程曲线呈现通过原点的直线分布规律，且不同测点处的差异性很小。
4 泄压空间内爆炸 4.1 初始条件设置
3.1泄压空间的尺寸及测点布置与第 节的设置相同，仅在泄压空间长度方向的某一壁面上设8置了一个方形泄压口，如图 所示。方形泄压口320 mm，以下图中记为L320 mm。769 g的边长为柱形炸药放置在泄压空间中心位置，等效后的高100 mm，高度压、高密度气体具体参数如下：半径120 mm，ρ = 203.75 kg/m3，p = 3.822×108 Pa。计算9初始条件如图 所示，考虑到计算时间及精度的14.4万个（90×要求，经多次数值试验，网格数仍为40×40 ）。壁面边界条件设置为刚性反射边界条件，泄压口处边界条件设置为透射边界。
4.2 爆炸波传播过程
10图 所示为不同时刻输出时泄压空间内部10的爆炸初期压力分布云图。由图 可知，泄压空间与密闭空间内部爆炸波传播过程基本上是遵循相同的规律，最大差别是：由于高压气体从泄压口处泄出，使得泄压空间泄压口区域附近的压力降10（a）），但这种泄压过程对爆炸初期高强度低（图爆炸波的传播过程影响较小。 10图 泄压空间内爆炸初期压力分布云图Fig.10 Pressure distribution at early stage in venting space
4.3 爆炸载荷特性
11图 所示为泄压空间内部爆炸壁面上所有11（a）可知，测点超压及冲量时间历程曲线。由图不同测点处爆炸前期爆炸波超压峰值大小有所不11（a）中同，而准静态超压时间历程几乎相同（图绿色粗实线）。这说明泄压空间内爆炸形成的准静态超压时间历程在空间上是近似均匀的。准静态超压曲线遵循指数衰减规律，与文献［22］中基于大量实验数据的假定一致。
11（b由图 ）可知，泄压空间内爆炸壁面测点冲量时间历程曲线呈现通过原点的抛物线分布规1，3，4号律，且不同测点处存在一定的差异性。1（x=1 350）内，冲量时间历测点处于同一切平面2，5 2（x=程基本相同。 号测点处于同一切平面1 650）内，冲量时间历程基本相同。而由于切平1 2面 相较于切平面 更靠近泄压口，使得其冲量时间历程衰减速率相对快一些。
4.4 泄压口位置的影响
本小节主要探讨泄压口位置对泄压空间内爆3.1炸载荷的影响规律。泄压空间的尺寸与第 节3的设置相同，但设定了 种不同的泄压口位置，如12图 所示分别为：泄压空间x方向半长度平面偏450 mm 0 mm、偏 450 mm。图 13左 、中间 右 所示3为不同泄压口位置爆炸工况下典型测点 的超压及冲量时间历程曲线。13由图 可知，泄压口位置对泄压空间内爆炸准静态超压载荷几乎没有影响，而对冲量时间历程产生了一定的影响，缩小泄压口与炸药的相对距离能降低冲量的大小。
4.5 泄压口大小的影响
本小节主要探讨泄压口大小对泄压空间内爆3.1炸载荷的影响规律。泄压空间的尺寸与第 节3的设置相同，但设定了 种不同边长的方形泄压L80 mm，L160 mm L320 mm。图口，分别记为 和14 3所示为不同泄压口大小爆炸工况下典型测点的超压时间历程曲线。
14由图 可知：泄压口大小对泄压空间内爆炸准静态超压时间历程影响显著，泄压口边长越大， 14（a））；泄压口大小对泄准静态超压衰减越快（图压空间内爆炸冲量时间历程影响显著，增大泄压14（b））。通过分口大小能显著降低冲量大小（图析可知，这主要是由于随着泄压口边长的增大，在
相同爆炸波强度条件下，单位时间内从泄压口处
泄出的能量越多。
5 连通空间内爆炸 5.1 初始条件设置
2连通空间由左、右 个尺寸相同的约束空间和中间方形连接导管组成，连接导管截面边长为400 mm 4 点（No.1~No.4）对。壁面上设置 个测 爆15所示。769 g炸超压时间历程进行输出，如图 柱形炸药位于左约束空间内部，等效后的具体参数4.1 16（a）所示。在计同第 节。计算初始条件如图算过程中，采取正交规则网格，考虑计算时间及精
度的要求，经多次数值试验，单个网格尺寸为 20 mm×20 mm×20 mm，如图16（b）所示。壁面边界条件设置为刚性壁面反射边界。
5.2 爆炸波传播过程
17图 所示为连通空间内爆炸初期压力分布17云图。由图 可以分析出爆炸初期爆炸波的传播过程：高压气体首先进行三维柱对称自由膨胀，爆炸波初次到达壁面时发生正规则反射，部分高压气体从泄压导管处泄出到右约束空间，使得泄17（a））；由于四周压导管区域附近的压力降低（图
壁面约束，爆炸波在左约束空间角隅附近区域开17（b始形成压力汇聚现象（图 ））；爆炸波通过泄压导管传播到右约束空间内部，并在右约束空间17（c））；爆炸波到达右约内部沿长度方向传播（图束空间长度方向端面后进行反射，并在右约束空间角隅附近区域形成压力汇聚现象，左约束空间内部由于爆炸波的多次反射在角隅处形成压力汇17（d））。聚现象（图
5.3 爆炸载荷特性
18图 分别给出了左约束空间（爆炸空间）内1，2爆炸壁面测点 的超压及冲量时间历程曲线。18 3.3由图 可知，左约束空间内爆炸载荷特征与 小
节密闭空间内爆炸具有相同的特征：不同测点具有几乎相同的准静态超压峰值及冲量时间历程。
19图 分别给出了右约束空间（容爆空间）内3 4爆炸壁面测点 ，的超压及冲量时间历程曲线。19 3.3由图 可知，右约束空间内爆炸载荷特征与小节密闭空间内爆炸具有一个相同的特征：不同测点具有几乎相同的准静态超压峰值，而不同测点处的冲量时间历程却具有一定的差异性。
20图 所示为连通空间内爆炸壁面上所有测20（a）可知，壁面不同点超压时间历程曲线。由图测点位置处爆炸前期爆炸波超压峰值大小有所不20（a）中绿色同，而准静态超压峰值几乎相同（图粗实线）。这说明连通空间内爆炸形成的准静态超压峰值在空间上是均匀的，这一特征类似于密闭空间内爆炸。然而，不同测点的冲量时间历程存在较大的差异性，左约束空间内部测点的冲量20（b））。明显高于右约束空间内部测点的冲量（图分析发现，这主要是由于爆炸发生在左约束空间，较多的能量集中在左约束空间。
6结论
FORTRAN WENO本文基于 平台，采用了三阶有限差分格式，自主开发了约束空间内部的爆炸波高精度数值计算程序，并开展了密闭空间、泄压空间和连通空间内部的爆炸波数值计算，分析了约束空间内部的爆炸波传播路径及爆炸载荷特性，通过研究主要得到如下结论： 1 ）密闭空间内部的爆炸冲量时间历程曲线呈现直线分布规律，泄压空间内部的爆炸冲量时间历程曲线呈现抛物线分布规律。2 ）泄压口位置对泄压空间内部的爆炸准静态超压载荷几乎无影响，而对冲量时间历程则具有一定的影响，缩小泄压口与炸药的相对距离可降低冲量的大小；泄压口大小对泄压空间内部的爆炸载荷影响显著，增大泄压口大小可显著加快准静态压力衰减速率，降低冲量大小。3 ）测点位置对连通空间内部的爆炸准静态超压载荷几乎无影响，而对冲量时间历程的影响则较显著，炸药所在约束空间的内部冲量明显高于爆炸波泄入空间的内部冲量。
参考文献：
［1］ ESPARZA E D，BAKER W E，OLDHAM G A. Blast pressures inside and outside suppressive structures： ADA025504［R］. San Antonio， TX： Southwest Re⁃ search Institute，1975. ［2］ KEENAN W A，TANCRETO J E. Blast environment from fully and partially vented explosions in cubicles： ADA019026［R］.［S.l.］：Department of the Army Pica⁃ tinny Arsenal，1975. ［3］ EDRI I，SAVIR Z，FELDGUN V R，et al. On blast pressure analysis due to a partially confined explosion： I. Experimental studies［J］. International Journal of
74 Protective Structures，2011，2（1）：1-20. ［4］ WU C Q ，LUKASZEWICZ M，SCHEBELLA K，et al. Experimental and numerical investigation of confined explosion in a blast chamber［J］. Journal of Loss Pre⁃ vention in the Process Industries， 2013， 26（4）： 737-750. ［5］ 胡洋，朱建芳，朱锴.长方体单腔室空腔环境内爆炸
J］. 2016，36（3）：效应的实验研究［ 爆炸与冲击， 340-346. HUY ZHU J F ，ZHU K. Experimental study on ex⁃
， plosion effect in a closed single rectangular cavity［J］. Explosion and Shock Waves，2016，36（3）：340-346 （in Chinese）. ［6］ 侯海量，朱锡，李伟，等.舱内爆炸冲击载荷特性实验研究［J］. 船舶力学，2010，14（8）：901-907. HOU H L ，ZHU X， LIW ，et al. Experimental studies on characteristics of blast loading when exploded in⁃ side ship cabin［J］. Journal of Ship Mechanics，2010， 14（8）：901-907（in Chinese）. ［7］ 孔祥韶，吴卫国，李俊，等.角隅结构对舱内爆炸载
J］. 2012，53（3）：荷影响的实验研究［ 中国造船， 40-50. KONG X S， WU W G LIJ ，et al. Experimental re⁃ ， search of influence of corner structure on blast loading under inner explosion ［J］. Shipbuilding of China， 2012，53（3）：40-50（in Chinese）. ［8］ 侯海量，朱锡，梅志远.舱内爆炸载荷及舱室板架结构的失效模式分析［J］. 爆炸与冲击，2007，27（2）： 151-158. HOU H L ，ZHU X，MEI Z Y. Study on the blast load and failure mode of ship structure subject to internal explosion［J］. Explosion and Shock Waves，2007，27 （2）：151-158（in Chinese）. ［9］ 孔祥韶，吴卫国，李晓彬，等.舰船舱室内部爆炸的J］. 2009，4（4）：数值模拟研究［ 中国舰船研究， 7-11. KONG X S， WUWG LI X B ，et al. Numerical simu⁃
， lation of cabin structure under inner explosion［J］. Chi⁃ nese Journal of Ship Research，2009，4（4）：7-11（in Chinese）. ［10］ 丁阳，陈晔，师燕超.室内爆炸超压荷载简化模型［J］. 工程力学，2015，32（3）：119-125，133. DING Y，CHEN Y，SHI Y C. Simplified model of overpressure loading caused by internal blas［t J］. En⁃ gineering Mechanics，2015，32（3）：119-125，133 （in Chinese）. ［11］ 樊壮卿，王伟力，黄雪峰，等.典型舱室内爆炸仿真分析［J］. 工程爆破，2015，21（3）：13-17. FAN Z Q ，WANG W L，HUANG X F，et al. Simula⁃ tion analysis on typical cabin internal explosion［J］. 12 Engineering Blasting，2015，21（3）：13-17（in Chi⁃ nese）. ［12］ LIU X D ，OSHER S，CHAN T. Weighted essentially non-oscillatory schemes［J］. Journal of Computation⁃ al Physics，1994，115（1）：200-212. ［13］ JIANG G S，SHU C W. Efficient implementation of weighted ENO schemes［J］. Journal of Computational Physics，1996，126（1）：202-228. ［14］ SHU C W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws［M］//QUARTERONI A. Advanced numerical approximation of nonlinear hyperbolic equations. Berlin Heidelberg： Springer， 1998： 325-432. ［15］ SHU C W. High order weighted essentially nonoscilla⁃ tory schemes for convection dominated problems［J］. Siam Review，2009，51（1）：82-126. ［16］ HARTEN A，ENGQUIST B，OSHER S，et al. Uni⁃ formly high order accurate essentially non-oscillatory schemes，III［J］. Journal of Computational Physics， 1987，71（2）：231-303. ［17］曹乐. level set利用 方法捕捉气、水界面的三维数值研究［D］.合肥：中国科学技术大学，2009. CAO L. Three-dimensional computations on captur⁃ ing of gas-water interface by level set method［D］ .He⁃ fei：University of Science and Technology of China， 2009（in Chinese）. ［18］ 谢昌坦. R-M气—水可压缩流物质界面的 不稳定性研究［D］.哈尔滨：哈尔滨工程大学，2011. XIE C T. The study on R-M instability of the material interface in gas-water compressible flow［D］. Har⁃ bin：Harbin Engineering University，2011 （in Chi⁃ nese）. ［19］ SHU C W ，OSHER S. Efficient implementation of es⁃ sentially non-oscillatory shock-capturing schemes ［J］. Journal of Computational Physics， 1988， 77 （2）：439-471. ［20］ SOD G A. A survey of several finite difference meth⁃ ods for systems of nonlinear hyperbolic conservation laws［J］. Journal of Computational Physics，1978，27 （1）：1-31. ［21］ WOODWARD P，COLELLA P. The numerical simu⁃ lation of two-dimensional fluid flow with strong shocks ［J］. Journal of Computational Physics， 1984， 54 （1）：115-173. ［22］ ANDERSON C E Jr，BAKER W E，WAUTERS D K，et al. Quasi-static pressure，duration，and im⁃ pulse for explosions（e.g. HE）in structures［J］. Inter⁃ national Journal of Mechanical Sciences，1983，25 （6）：455-464.