基于Mie-Grüneisen状态方程的水下爆炸数值模拟

吴宗铎，严 谨，蒋 颉，董美余，黄 技

(1. 广东海洋大学 海洋工程学院，广东 湛江 524088；2. 中国船舶研究设计中心，上海 201108；3. 上海振华重工集团股份有限公司，上海 200125)

基于Mie-Grüneisen状态方程的水下爆炸数值模拟

吴宗铎1，严 谨1，蒋 颉2，董美余3，黄 技1

(1. 广东海洋大学 海洋工程学院，广东 湛江 524088；2. 中国船舶研究设计中心，上海 201108；3. 上海振华重工集团股份有限公司，上海 200125)

由于Mie-Grüneisen状态方程形式比较复杂，给界面的处理带来很大难度。本文通过对Euler方程做分离变量和引入质量分数，完成了Mie-Grüneisen状态方程下的多介质可压缩流动的数值模拟，使计算过程得到简化，并通过算例验证了该方法的可靠性。在多项式形式的Mie-Grüneisen状态方程下，利用该方法模拟了球形炸药在水中爆炸后，气体和水相互作用的近场情况。在计算模型中引入的气体质量分数，很好反映了流场中不同区域内气体、水及水气并存3种不同的状态。

水下爆炸；数值计算；Mie-Grüneisen状态方程；多介质

0 引 言

由于水下爆炸在军事和国防领域有着极大的重要性，因而一直受到大家的关注。多数时候，炸药都以球形装药的形式出现。对于球形炸药来说，得到其冲击波压力并不容易。早期对冲击波压力的研究以实验为主[1 – 2]。其后，虽然出现了许多数值模拟，但水下爆炸的实验仍未中断[3 – 4]。但考虑保密性的要求，很多结果未能发表。

然而，由于实验条件成本大，实施难，水下爆炸已经越来越多的被数值模拟所代替。但球形炸药的三维特性使得它在普通的数值模拟中很难得到广泛应用。其中很大一个原因就是普通的数值模拟中设置三维计算网格将消耗更多的计算机资源。因此，关于水下爆炸的研究工作多基于商业软件，如Autodyn[5 – 6]，LS-DYNA[7]，Abaqus[8]，MSC-Dytran[9]等。而目前的水下爆炸数值模拟仍然多停留在二维平面上，如张阿曼的无网格SPH方法[10]、Daramizadeh[11]的五方程方法。除此以外，还有一个原因也对球形炸药的数值模拟起着限制作用，那就是涉及到复杂状态方程下的Riemann问题求解。目前对于水下爆炸问题，采用的状态方程多为形式复杂Mie-Grüneisen状态方程[5, 7, 10]，给冲击波的模拟造成很大的难度。

本文将求解坐标系由传统的二维网格体系转到球坐标体系下。由于药球爆炸后冲击波传播的球面对称形式，冲击波的参数仅和传播距离有关[12]。这样，可将冲击波物理参数表示成球坐标系半径的函数，无需建立三维的直角坐标网格来求解，为计算带来很大的方便。对于Riemann问题的处理，本文在Mie-Grüneisen状态方程下对能量守恒Euler方程做变量分离，通过近似处理把参考压力、参考单位内能等复杂的物理量分离出来，使Euler方程组简单易于计算。利用这种算法，对水下爆炸的近场情况进行数值计算。在激波的传播过程中，整个流场被视为一个理想气体和水的混合体，当中的气体和水用质量分数加以区分。由于药球爆炸后，水和水蒸气共存的状态大量存在于流场周围，该方法能很好地反映这一状态，同时也能较好地计算出水下爆炸周围流场的情况。

1 控制方程

水下爆炸问题中，无论是炸药还是介质，其物理参数密度ρ、压力p和单位体积内能ρe都可写成Mie-Grüneisen状态方程的形式[12]，其表达式如下：

其中，Г为Mie-Grüneisen系数，可以表示成

因此，Г可以写成密度ρ的函数。pref和eref为参考点的压力和单位体积内能，与密度ρ有关。假设流场无粘，那么根据质量、动量和能量守恒关系，可以得出欧拉方程：

式中：ρ，p，u，E分别为密度，压力，速度和单位体积能量，并且对于单位体积能量E，可以写成内能和动能的总和，于是有

在炸药和水的接触间断面上，大多数物理参数为间断状态但流体的压力和速度仍然保持连续，间断的物理参数导数可以视为无穷大，因此u和p的导数相对其他导数来说可以忽略不计，那么式（3）中质量守恒和能量守恒可以近似表示成：

将方程（1）代入式（3）中的能量守恒方程，有

由于Г，pref和eref都可以写成只和密度ρ有关的数，再忽略p的变化，因此，可以对式（5）作变量分离：

这样Mie-Gruneisen状态方程下的m种介质的控制方程可以表示成包含ρ，p，u，E和Г，pref，eref以及质量分数αi的如下的欧拉关系式。方程（3）中的变量U可以写成：

通量F则表示成：

2 计算方法

本文采用MUSCL-TVD格式的有限体积法进行计算。假设问题为一维Euler方程，其有限体积形式有：

3 计算实例

半径为R0裸药球在无限水域中爆炸，假设流体无粘且药球质量分布均匀，t=0时刻，炸药爆炸并形成高压气体，气体密度ρ0等于炸药球的密度，即1 630 kg/m3，爆压pH=2.1×1010Pa，气团半径等于药球半径[1]。假设TNT炸药爆炸后形成的气体为理想气体，状态方程如下：

式中：μ=ρ/ρ0–1，其系数如表1所示。

这样，转化成如（1）所示的Mie-Grüneisen的形式以后，相关参数为：

表 1 炸药和水的状态方程参数Tab. 1 The EOS parameters for explosive and water

图1为压力峰值与文献[2]的对比结果。从图中可以看出，2条压力峰值曲线与大致保持一致，且在距爆点约12R0处吻合得比较好。当冲击波过了12R0以后，冲击波的压力衰减速度降低，曲线将有所缓和，与本文的计算存在少量差距。而在距爆点很近的地方（6R0以内），误差会有所增加。由于这部分区域的爆炸冲击波的作用很强，压力变化非常剧烈，有些实验结果甚至未能在6R0以内给出合理的经验公式[1]。

图2为距爆心不同距离处的压力峰值曲线，横坐标是测点与爆心的距离R和初始半径R0的比值。本文计算结果与文献[2]的经验公式相比，近场压力值总体有些偏大。但基本反映了压力在冲击波前后的变化情况。在冲击波传播到某点时，该点压力瞬间上升到峰值；随后冲击波继续前进，此时由于冲击波作用减弱该点压力则迅速下降。随后，向内收缩的稀疏波聚集到爆心处又重新扩散开来，造成压力的脉动变化，形成第2次峰值。

图3为不同时刻压力的分布图，从图中可以看出，随着时间的推移，压力曲线逐渐变得平稳，并在最大峰值过后一定距离出现了压力脉动，形成二次脉动压力。图4给出了不同位置的二次压力峰值，可以看出其数值比较小。值得注意的是，由于本文的计算是基于能量守恒的假设，未考虑能量的损耗，因此本文的计算结果比实际情况要大，由脉动压力形成的第2次的压力峰值也比较明显。实际情况下，第2次压力峰值很难在曲线上清楚显现出来[1]。

图5理想气体的质量分数α和质量ρα在不同时刻的变化情况，从图6中可以看出，冲击波传播开来以后，中央的气团逐渐由单一的气体状态变为气体和液态水并存的一种状态。这情况与实际的气水模型存在一定差别，实际的气水模型中并未有爆轰气体溶解在水中，因此本计算模型仍然存在着一些数值上的不准确，但并不影响冲击波压力的计算。

4 结 语

利用本文的分离变量方法能很好地解决Mie-Grüneisen状态方程下的多介质可压缩流动的激波问题，并且能有效地简化计算过程。该方法在处理水下爆炸近场问题的时候能取得很好的效果。按本文的计算方法处理得到的计算结果和文献[1]对比，取得了不错的效果，并对计算第2次的脉动压力有着很好的帮助。

炸药在水中爆炸以后，水的状态包括液态水、水蒸气和液气混合态。再加上爆炸产生的气体、液态水、气体和气水混合体共存于流场中。利用质量分数可以区分开流场中的各种成分。炸药爆炸后，早期形成高密度气团，随后气团扩散很快，在气团扩散过程中，气体的质量分数和密度都急剧下降，最后形成一个中心大部分区域稀薄，两边小部分区域稠密的状态。

由于本文仅考虑了动量和能量守恒的理想情况，而实际情况比较复杂，能量损耗不可避免，因此，本文所计算的压力值比实际数值大，冲击波的传播速度也比实际要快。

[ 1 ]COLE R. H. Underwater explosions[M]. New York: Dover Publications, INC. 1965.

[ 2 ]ZAMYSHLYAYEV B V, YAKOVLEV Y S. Dynamic loads in underwater explosion, AD-757183 [R]. Naval Intelligence Support Center, 1978.

[ 3 ]杨振, 沈晓乐. 爆破战斗部水中兵器爆炸威力评定方法研究[J]. 爆破, 2015, 32(2): 51–53. YANG Zhen, SHEN Xiao-le. Research on evaluation method of underwater blast brisance of weapon's explosive[J]. Blasting, 2015, 32(2): 51–53.

[ 4 ]鲁忠宝, 胡宏伟, 刘锐, 等. 典型装药水下爆炸的殉爆规律研究[J]. 鱼雷技术, 2014, 22(3): 230–235. LU Zhong-bao, HU Hong-wei, LIU Riu, et al. Research on law of sympathetic detonation of typical charge[J]. Torpedo Technology, 2014, 22(3): 230–235.

[ 5 ]张社荣, 李宏璧, 王高辉, 等. 水下爆炸冲击波数值模拟的网格尺寸确定方法[J]. 振动与冲击, 2015, 34(8): 93–100. ZHANG She-rong, LI Hong-bi, WANG Gao-hui, et al. A method to determine mesh size in numerical simulation of shock wave of underwater explosion[J]. Journal of Vibration and Shock, 2015, 34(8): 93–100.

[ 6 ]姜涛, 王桂芹, 詹发民, 等. 基于AUTODYN的潜艇典型舱段水中爆炸冲击损伤研究[J]. 爆破器材, 2015, 44(6): 61–64. JIANG Tao, WANG Gui-qin, ZHAN Fa-min, et al. Impact damage analysis of typical submarine compartment subjected to underwater blasting based on AUTODYN[J]. Explosive Materials, 2015, 44(6): 61–64.

[ 7 ]梅群, 侯中华, 朱俊锋, 等. 水下爆炸冲击波压力时程的数值模拟[J]. 河南科技大学学报(自然科学版), 2010, 31(4): 57–59 MEI Qun, HOU Zhong-Hua, ZHU Jun-Feng, et al. Numerical simulation of underwater explosion shock wave[J]. Journal of Henan University of Science and Technology: Natural Science, 2010, 31(4): 57–59.

[ 8 ]贾则, 权琳, 张姝红, 等. 舰艇结构在水中兵器静态爆炸作用下的冲击响应[J]. 舰船科学技术, 2015, 37(3): 55–58. JIA Ze, QUAN Lin, ZHANG Shu-hong, et al. The research in shock response of warship under static explosion by underwater weapon[J]. Ship Science and Technology, 2015, 37(3): 55–58.

[ 9 ]袁建红, 朱锡, 张振华. 水下爆炸载荷数值模拟方法[J]. 舰船科学技术, 2011, 33(9): 18–23. YUAN Jian-hong, ZHU Xi, ZHANG Zhen-hua. Numerical simulaton method study of underwater explosion load[J]. Ship Science and Technology, 2011, 33(9): 18–23.

[10]ZHANG A-man, YANG Wen-shan, YAO Xiong-liang. Numerical simulation of underwater contact explosion[J]. Applied Ocean Research, 2012, 34(1): 10–20.

[11]A. DARAMIZADEH a, M. R. ANSARI. Numerical simulation of underwater explosion near air–water free surface using a five-equation reduced model[J]. Ocean Engineering, 2015, 110(12): 25–35.

[12]张宝平, 张庆明, 黄风雷. 爆轰物理学[M]. 北京: 兵器工业出版社, 2006.

The numerical simulation of underwater explosion based on Mie-Grüneisen equation of state

WU Zong-duo1, YAN Jin1, JIANG Jie2, DONG Mei-yu3, HUANG Ji1
(1. Guangdong Ocean University School, College of Ocean Engineering, Zhanjiang 524088, China; 2. China Ship Development and Design Institute, Shanghai 201108, China; 3. Shanghai Zhenghua Heavy Industries Co., Ltd, Shanghai 200125, China)

Due to the complicated form of Mie-Grüneisen equation of state, it is difficult to handle the interface well. We separate several variables from the Euler equation and introduce mass fractions, manage to achieve the numerical simulation of multi-medium compressible flows under Mie-Grüneisen equation of state, while the process of calculation is simplified here. The approach was verified well in our case. Then under the polynomial equation of state, which is a general form of Mie-Grüneisen equation of state, we used this approach to simulate the underwater explosion of a sphere bomb, and obtained the near-field character of gas-water interaction. By introducing the mass fraction, the three states of different region, which include gas, water and gas-water mixture, are revealed clearly.

underwater；numerical calculation；Mie-Grüneisen equation of state；multi-medium

O382.1

A

1672 – 7649(2017)07 – 0029 – 05

10.3404/j.issn.1672 – 7649.2017.07.006

2016 – 07 – 12；

2016 – 11 – 10

广东青年创新人才类资助项目(2014KQNCX086)；渔船渔港设施装备标准化与安全保障能力研究资助项目(GDOU2016050258)

吴宗铎(1984 – )，男，博士，讲师，主要从事激波间断及水下爆炸冲击波数值模拟。