基于新型轴对称无网格方法的水下爆炸冲击波和气泡运动数值模拟

摘要

采用无网格光滑粒子流体动力学（smooth particle hydrodynamics， SPH）方法，结合粒子体积自适应技术（volume adaptive scheme， VAS）构建高精度、高效率的新型轴对称SPH模型，实现水下爆炸冲击波和气泡运动模拟。在冲击波传播阶段，通过与经验公式及文献理论值对比，验证了该轴对称SPH模型的精度。针对气泡脉动模拟，选取典型的双气泡耦合问题，分析气泡形态和压力载荷的力学特性，讨论了两气泡之间的距离参数对气泡相互作用过程的影响。结果表明，该新型轴对称SPH模型可实现水下爆炸过程压力载荷和气泡形态演化的精确预报。

关键词

水下爆炸; 轴对称SPH模型; 冲击波; 双气泡耦合

随着现代水下攻击型武器的制导精度和攻击速度等的大幅度提升，攻击形式更加多样化，水下爆炸载荷对水面舰船或潜艇等的威胁也不断提高。因此，精准预报水下爆炸过程的载荷变化对舰船等水下兵器的设计优化至关重要^［1］。对此，本文针对水下爆炸过程包括冲击波传播阶段及典型的双气泡相互作用展开研究，分析压力载荷等动力特性，为提升我国相关军事装备的战斗力、生命力提供数据支持。

早期研究主要围绕水下爆炸冲击载荷展开，随着研究的深入，学者们发现气泡脉动是结构损伤的重要因素，为此开展了大量水下爆炸试验^［2- 3］，并通过试验数据拟合得到水下爆炸经验公式，但实尺度的水下爆炸试验将耗费巨大的人力和物力，且可重复性低。目前，学者们采用TNT装药在试验水箱内开展爆炸试验^［3］，但受限于测量仪器，大多只能实现装药较远处的流场压力测量，装药附近和爆炸气泡内的压力变化无法监测；且对于近场冲击波传播阶段，由于持续时间为毫秒级，所以试验数据的获取难度较大。

针对以上试验研究的不足，理论和数值方法逐渐脱颖而出^［4］。张阿漫等提出了气泡统一方程^［5］，解决了不同源、不同环境和不同尺度下的空化气泡、水下爆炸气泡、微纳米气泡等气泡动力学特性的理论预测难题。基于不同数值方法，水下爆炸过程中冲击波传播及气泡脉动特性也得到广泛和深入的研究^［6-8］。为提升国防武器性能，精准有效应对多发、多处导弹同时爆炸，多气泡耦合的研究也日益受到关注^［9］。多气泡耦合涉及气泡膨胀、挤压、融合、撕裂等复杂现象，意味着数值算法的精度和稳定性亟需进一步提升。近年来，因无网格粒子方法具有擅长捕捉多相界面变形、破碎、掺混等优点，在水下爆炸研究领域逐渐发挥突出作用^［10-12］，使冲击载荷预报精度得到提升^［13］。因此，本文采用SPH（smooth particle hydrodynamics）方法，实现水下爆炸冲击波传播、典型的双气泡耦合场景的数值模拟。

目前，工程实际中的水下爆炸问题依赖于三维SPH计算模型进行模拟，计算成本高昂。考虑到如鱼雷、水雷、导弹等水中兵器的装药均为近似圆柱体，水下爆炸时均可以考虑为轴对称问题，实现从三维空间转化到二维平面，以节省计算量。Petschek等^［14］首次提出轴对称SPH模型并运用到多种材料高速冲击问题研究中，但该模型在水动力问题研究上的精度仍没有得到验证。在水动力问题上，Zhang等^［10］采用轴对称SPH模型实现水下爆炸冲击波传播模拟，但随着高压气体的膨胀，气体粒子分布稀疏，仅计算了较短时间，模型稳定性有待进一步提升。Brookshaw^［15］提出的拉格朗日形式的轴对称SPH模型，被学者采纳并应用于物体入水^［16］、溃坝等典型的水动力问题研究中，该模型计算精度得到验证，但其在对称轴处的数值奇点问题没得到有效解决。因此，实现轴对称SPH模型对水下爆炸过程的模拟存在两大难点：一是水下爆炸气泡脉动等强可压缩问题的计算；二是对称轴处存在奇异性，计算精度和稳定性受限。对此，本文将新型轴对称SPH算法拓展应用于气泡脉动和双气泡耦合作用研究，针对传统SPH算法的不足，构建高精度、高效率的新型轴对称SPH模型，实现水下爆炸过程冲击波传播和气泡脉动过程的数值模拟。

1 轴对称无网格计算模型

1.1　控制方程

本文模拟的炸药爆破及高能气泡都是关于中心轴对称的，为节约计算成本，采用新型SPH轴对称模型^［17］实现三维平面水下爆炸问题冲击波传播和气泡运动过程模拟，如图1所示，图中$\mathrm{\Delta }r$为初始粒子间距。控制方程如下：

$$\left\{\begin{array}{l}\frac{\mathrm{D}\rho }{\mathrm{D}t}=-\frac{\rho }{r}\left[\frac{\partial }{\partial r}\left(r{u}_r\right)+\frac{\partial u_{\theta }}{\partial \theta }+r\frac{\partial u_z}{\partial z}\right]\\ \frac{\mathrm{D}{u}_r}{\mathrm{D}t}=-\frac{\mathrm{1}}{\rho }\frac{\partial p}{\partial r}\\ \frac{\mathrm{D}{u}_{\theta }}{\mathrm{D}t}=-\frac{\mathrm{1}}{\rho r}\frac{\partial p}{\partial \theta }\\ \frac{\mathrm{D}{u}_z}{\mathrm{D}t}=-\frac{\mathrm{1}}{\rho }\frac{\partial p}{\partial z}+g\\ \frac{\mathrm{D}e}{\mathrm{D}t}=-\frac{p}{\rho r}\left[\frac{\partial }{\partial r}\left(r{u}_r\right)+\frac{\partial u_{\theta }}{\partial \theta }+r\frac{\partial u_z}{\partial z}\right]\end{array}\right.\underset{\partial f_{\theta }/\partial \theta =\mathrm{0}}{\stackrel{u_{\theta }=\mathrm{0}}{\to }}\left\{\begin{array}{l}\frac{\mathrm{D}\rho }{\mathrm{D}t}=-\rho \frac{u_r}{r}-\rho \left(\frac{\partial u_r}{\partial r}+\frac{\partial u_z}{\partial z}\right)\\ \frac{\mathrm{D}{u}_r}{\mathrm{D}t}=-\frac{\mathrm{1}}{\rho }\frac{\partial p}{\partial r}\\ \frac{\mathrm{D}{u}_z}{\mathrm{D}t}=-\frac{\mathrm{1}}{\rho }\frac{\partial p}{\partial z}+g\\ \frac{\mathrm{D}e}{\mathrm{D}t}=-\frac{p}{\rho }\frac{u_r}{r}-\frac{p}{\rho }\left(\frac{\partial u_r}{\partial r}+\frac{\partial u_z}{\partial z}\right)\end{array}\right.$$

（1）

式中：$\frac{\mathrm{D}f}{\mathrm{D}t}$为物质导数；$\rho$、$p、u$、$e$为密度、压力、速度和单位质量内能；$r\mathrm{、}z、\theta$分别为径向、轴向及周向坐标；$g$为重力加速度，取-9.8 m⋅s^-2。公式（1）左侧给出柱坐标系下无黏和可压缩流体的控制方程，根据轴对称假设，周向速度及沿周向方向变量的偏导数为零，即$u_{\theta }=\mathrm{0}$，$\partial f_{\theta }/\partial \theta =\mathrm{0}$，因此，可得到简化后的轴对称问题的控制方程（公式（1）右侧）。本文轴对称模型原理如图1所示，$z‒r$平面上的每一个拉格朗日粒子都代表$O‒r\theta z$平面的一个圆环，初始质量$m_{\mathrm{i}\mathrm{n}\mathrm{t}}$分布沿粒子径向坐标增大，$m_{\mathrm{i}\mathrm{n}\mathrm{t}}=\mathrm{2}\mathrm{\pi }\rho r_i{V}_{\mathrm{i}\mathrm{n}\mathrm{t}}$（$V_{\mathrm{i}\mathrm{n}\mathrm{t}}$为圆环的截面面积，也是二维SPH模拟的粒子体积），从而得以将柱坐标系$O‒r\theta z$立体空间简化到$z‒r$平面。

图1  三维轴对称问题简化到二维平面的示意图

Fig. 1  Simplification from a 3D cylindrical coordinates system to a 2D axisymmetric plane

将公式（1）右侧方程进行粒子近似，可得本文轴对称SPH模型控制方程^［17］如下：

（2）

式中：$W_{ij}$为核函数；$V$为粒子体积；速度$\mathit{u}=(u_r,u_z)$，$\mathit{r}$为粒子位移；重力加速度$\mathit{g}=\left(g_{\mathrm{r}},g_{\mathrm{z}}\right)=(\mathrm{0},-\mathrm{9.8})$ m·s^-2。为保持计算稳定性，在连续方程中引入密度耗散项$\mathit{D}$，及在动量方程中引入人工黏性项$\mathit{\Pi }$^［17］$\mathrm{。}$同时，为实现计算过程中粒子分布一直均匀，本文采用粒子位移修正技术（公式（2）最后一行公式），具体参数细节可参考文献［18］。速度散度${}_i^L$、压力梯度${}_i^{L_{\mathrm{2}}}$的离散采用重构矩阵^［19］进行修正，以提高计算精度。

（3）

此外，为了对控制方程（公式（2））解耦，引入状态方程。在冲击波传播阶段，本文采用Mie-Gruneisen状态方程^［20］将水模拟为可压缩流体。

当$\mu >\mathrm{0}$时，水处于压缩状态，则有

$$p=\frac{{\rho }_{\mathrm{0}}{c_{\mathrm{0}}}^{\mathrm{2}}\mu \left[\mathrm{1}+\left(\mathrm{1}-\frac{{\gamma }_{\mathrm{0}}}{\mathrm{2}}\right)\mu -\frac{a}{\mathrm{2}}{\mu }^{\mathrm{2}}\right]}{{\left[\mathrm{1}-\left(S_{\mathrm{1}}-\mathrm{1}\right)\mu -S_{\mathrm{2}}\frac{{\mu }^{\mathrm{2}}}{\mu +\mathrm{1}}-S_{\mathrm{3}}\frac{{\mu }^{\mathrm{3}}}{{\left(\mu +\mathrm{1}\right)}^{\mathrm{2}}}\right]}^{\mathrm{2}}}+\left({\gamma }_{\mathrm{0}}+a\mu \right)e$$

（4）

当$\mu <\mathrm{0}$时，水处于膨胀状态，状态方程为

$$p={\rho }_{\mathrm{0}}{c_{\mathrm{0}}}^{\mathrm{2}}\mu +\left({\gamma }_{\mathrm{0}}+a\mu \right)e_{}$$

（5）

式（4）、（5）中：$S_{\mathrm{1}}\mathrm{、}{S}_{\mathrm{2}}\mathrm{、}{S}_{\mathrm{3}}$为材料的$\mathrm{H}\mathrm{u}\mathrm{g}\mathrm{o}\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{t}$常数；$\mu =\rho /{\rho }_{\mathrm{0}}-\mathrm{1}$，${\rho }_{\mathrm{0}}$为流体质点的初始密度；$c_{\mathrm{0}}$为真实声速；${\gamma }_{\mathrm{0}}$为$\mathrm{G}\mathrm{r}\ddot{\mathrm{u}}\mathrm{n}\mathrm{e}\mathrm{i}\mathrm{s}\mathrm{e}\mathrm{n}$系数；$a$为体积修正系数；$e$为单位质量内能；具体参数取值见文献［20］。

TNT爆轰后采用JWL（Jones-Wilkins-Lee）状态方程描述如下：

$$p=A\left(\mathrm{1}-\frac{\omega \eta }{R_{\mathrm{1}}}\right){\mathrm{e}}^{-\frac{R_{\mathrm{1}}}{\eta }}+B\left(\mathrm{1}-\frac{\omega \eta }{R_{\mathrm{2}}}\right){\mathrm{e}}^{-\frac{R_{\mathrm{2}}}{\eta }}+\omega \eta {\rho }_{\mathrm{0}}{e}_{}$$

（6）

式中：$\eta =\rho /{\rho }_{\mathrm{0}}$；$A\mathrm{、}B\mathrm{、}{R}_{\mathrm{1}}\mathrm{、}{R}_{\mathrm{2}}\mathrm{和}\mathrm{\omega }$通过试验拟合得到；具体取值见文献［20］。

在气泡作用阶段，水粒子和气相粒子均采用Tait状态方程^［21］，对于水粒子，状态方程为

$$p=p_{\mathrm{0}}^{\mathrm{\text{'}}}\left[{\left(\frac{\rho }{{\rho }_{\mathrm{0}}}\right)}^{\gamma }-\mathrm{1}\right]+p_{\mathrm{b}},p_{\mathrm{0}}^{\mathrm{\text{'}}}=\frac{c_{\mathrm{0}}^{\mathrm{\text{'}}\mathrm{2}}{\rho }_{\mathrm{0}}}{\gamma }$$

（7）

对于气泡，状态方程写为

$$p=p_{\mathrm{0}}{\left(\frac{\rho }{{\rho }_{\mathrm{0}}}\right)}^{\gamma },p_{\mathrm{0}}=\frac{c_{\mathrm{0}}^{\mathrm{2}}{\rho }_{\mathrm{0}}}{\gamma }$$

（8）

式（7）、（8）中：$p_{\mathrm{0}}$、${\rho }_{\mathrm{0}}$分别为初始压力、密度；$p_{\mathrm{b}}$为背景压力；$c_{\mathrm{0}}^{\text{'}}$为人工声速^［21］，区别于冲击波阶段，该声速取值小于真实声速，在确保计算稳定情况下，可增大计算时间步，以减少计算时间。对于水粒子取$\gamma =\mathrm{7.15}$，对于气相粒子取$\gamma =\mathrm{1.25}$。

1.2　粒子体积自适应技术

在轴对称SPH算法中，粒子质量$m$代表整个流体环的质量，而粒子体积作为二维简化截面的积分微元，为流体环横截面的面积^［17］。假如不采用粒子撕裂/融合技术，在计算过程中每个粒子的质量保持不变（图1），而粒子体积的变化遵循以下规律：

$$V_i(t)=\frac{m_i}{\mathrm{2}\mathrm{\pi }{r}_i(t){\rho }_i(t)}$$

（9）

随着炸药引爆，炸药快速转换为高温高压气体引爆完成，气体/水体发生膨胀/收缩，即气体粒子径向距离会增大/减小，引起粒子体积的大幅变化。在SPH计算框架内，当相互作用的粒子对体积相差较大时，将引起较大的数值误差，因此确保计算过程粒子体积的均匀性是准确模拟气泡形态变化过程的关键。对此，本文引入一种新颖的粒子体积自适应算法^［21］，实现粒子体积的实时监测与调整。该算法的原理如图2所示。具体来说，在每个时间步对粒子体积进行监测，当气体粒子体积过度膨胀后（判断条件为$V\ge \mathrm{1.6}{V}_{\mathrm{0}}$，$V_{\mathrm{0}}$为初始体积），该粒子撕裂为4个小粒子。相反，当其体积过度收缩（$V<\mathrm{2}/\mathrm{3}{V}_{\mathrm{0}}$），该粒子与周围最近的粒子发生融合，形成更大粒子^［21］，根据质量守恒和动量守恒，撕裂/融合生成新粒子的物理量赋值见公式（10）。公式（10）中，下标$\mathrm{m}$代表体积较大的母粒子，$\mathrm{d}$代表体积较小的子粒子， $F_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{e}}^{-\mathrm{1}}$为状态方程的反函数，h为光滑长度。粒子数量随气泡体积而变化，即气体粒子数随着气泡膨胀而撕裂增多，随着气泡收缩而融合变少，从而实现大幅气泡脉动时，平均粒子体积和间距保持近似不变，进而实现了拉格朗日粒子法对气泡脉动的模拟。

图2  粒子体积自适应技术原理示意图

Fig. 2  Sketch of volume adaptive scheme

$$\begin{array}{l}\mathrm{撕裂}\left\{\begin{array}{l}{r}_{\mathrm{d}}=r_{\mathrm{m}}\pm \frac{\sqrt[]{V_{\mathrm{m}}}}{\mathrm{4}},z_{\mathrm{d}}=z_{\mathrm{m}}\pm \frac{\sqrt[]{V_{\mathrm{m}}}}{\mathrm{4}},\\ m_{\mathrm{d}}=\frac{\mathrm{2}\mathrm{\pi }{\rho }_{\mathrm{m}}{V}_{\mathrm{m}}}{\mathrm{4}}{r}_{\mathrm{d}},V_{\mathrm{d}}=\frac{m_{\mathrm{d}}}{\mathrm{2}\mathrm{\pi }{r}_{\mathrm{d}}{\rho }_{\mathrm{d}}},h_{\mathrm{d}}=h_{\mathrm{m}},\\ p_{\mathrm{d}}=p_m+<\nabla p{>}_{\mathrm{m}}^L\left({\mathit{r}}_{\mathrm{d}}-{\mathit{r}}_{\mathrm{m}}\right),{\rho }_{\mathrm{d}}=F_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{e}}^{-\mathrm{1}}\left(p_{\mathrm{d}}\right),{\mathit{u}}_{\mathrm{d}}={\mathit{u}}_{\mathrm{m}},e_{\mathrm{d}}=e_{\mathrm{m}}\end{array}\right.\\ \mathrm{融合}\left\{\begin{array}{l}{\mathit{r}}_{\mathrm{m}}=\frac{m_{\mathrm{d}}^{\mathrm{1}}{\mathit{r}}_{\mathrm{d}}^{\mathrm{1}}+m_{\mathrm{d}}^{\mathrm{2}}{\mathit{r}}_{\mathrm{d}}^{\mathrm{2}}}{m_{\mathrm{d}}^{\mathrm{1}}+m_{\mathrm{d}}^{\mathrm{2}}},{\mathit{u}}_{\mathrm{m}}=\frac{m_{\mathrm{d}}^{\mathrm{1}}{\mathit{u}}_{\mathrm{d}}^{\mathrm{1}}+m_{\mathrm{d}}^{\mathrm{2}}{\mathit{u}}_{\mathrm{d}}^{\mathrm{2}}}{m_{\mathrm{d}}^{\mathrm{1}}+m_{\mathrm{d}}^{\mathrm{2}}},\\ m_{\mathrm{m}}=m_{\mathrm{d}}^{\mathrm{1}}+m_{\mathrm{d}}^{\mathrm{2}},h_{\mathrm{m}}=h_{\mathrm{d}},V_{\mathrm{m}}=\frac{m_{\mathrm{m}}}{\mathrm{2}\mathrm{\pi }{r}_{\mathrm{m}}{\rho }_{\mathrm{m}}},\\ {\rho }_{\mathrm{m}}=\frac{m_{\mathrm{d}}^{\mathrm{1}}{\rho }_{\mathrm{d}}^{\mathrm{1}}+m_{\mathrm{d}}^{\mathrm{2}}{\rho }_{\mathrm{d}}^{\mathrm{2}}}{m_{\mathrm{d}}^{\mathrm{1}}+m_{\mathrm{d}}^{\mathrm{2}}},p_{\mathrm{m}}=F_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{e}}^{-\mathrm{1}}\left({\rho }_{\mathrm{m}}\right),e_{\mathrm{m}}=\frac{m_{\mathrm{d}}^{\mathrm{1}}{e}_{\mathrm{d}}^{\mathrm{1}}+m_{\mathrm{d}}^{\mathrm{2}}{e}_{\mathrm{d}}^{\mathrm{2}}}{m_{\mathrm{d}}^{\mathrm{1}}+m_{\mathrm{d}}^{\mathrm{2}}}\end{array}\right.\end{array}$$

（10）

1.3　对称轴处理

在传统轴对称SPH算法中，对称轴的处理一直是一个棘手的问题，在计算中，靠近对称轴$\left(r➝\mathrm{0}\right)$的粒子在另一侧$(r<\mathrm{0}$）会缺少邻居粒子，导致核函数截断（图3左侧）。当控制方程（公式（2））计算$\mathrm{1}/r_i$项时，对称轴附近$(r\mathrm{➝}\mathrm{0})$的粒子将存在数值奇点，使得计算稳定性和精度极为受限。为解决这个问题，本文新型轴对称SPH算法的另一新型数值技术是对对称轴处的粒子镜像处理。具体来说，在对称轴的另一侧（r$<\mathrm{0}$），在一个核函数半径的距离通过镜像流体粒子建立对称轴，如图3右侧所示。对称轴粒子的物理量通过镜像对应的水粒子得到^［21］，具体见公式（11）。公式（11）中，下标g代表虚粒子，$\mathrm{r}$代表实粒子。

$$\left\{\begin{array}{l}{m}_{\mathrm{g}}=m_{\mathrm{r}},h_{\mathrm{g}}=h_{\mathrm{r}}\\ v_{\mathrm{g}}^{}=-v_{\mathrm{r}}^{}\\ p_{\mathrm{g}}=p_{\mathrm{r}},{\rho }_{\mathrm{g}}={\rho }_{\mathrm{r}}\end{array}\right.$$

（11）

图3  新型轴对称SPH模型对称轴附近镜像虚粒子示意图

Fig. 3  Sketch of mirroring ghost particles near the axis in the novel SPH model

1.4　固壁边界条件施加

本文所有算例的左边界均为对称轴，对称轴的处理方法在上一节已进行论述；上下边界及右边界均为固壁边界。本文固壁边界条件的施加采用固定虚粒子法^［22］，将壁面边界离散为固定虚粒子，并与对称轴处理类似，固定虚粒子的布置宽度不小于一个核函数半径的距离（即2h）。固定虚粒子的位置固定不变，自然满足壁面不可穿透条件，其物理量通过与周围流体粒子插值得到。为避免压力波在壁面处反射，所有壁面边界处均施加海绵层边界^［23］。本文算例的计算域均取到足够大而不显著增加计算量，足以满足气泡长时间的形态发展。

2 水下爆炸冲击波计算和验证

为验证本文轴对称SPH方法预报的水下载荷的精度，本节采用水下爆炸冲击波压力时历曲线与水下爆炸经验公式拟合曲线^［24］进行对比验证。如图4所示，在计算域中心布置8.3 kg的TNT球形药包，半径$r_{\mathrm{0}}=\mathrm{106.4}\mathrm{m}\mathrm{m}$，初始粒子间距$r_{\mathrm{0}}/\mathrm{\Delta }r=\mathrm{0.005}\mathrm{3}$，模拟水域$\mathrm{10}\mathrm{m}\times \mathrm{6}\mathrm{m}$。

图4  冲击波传播过程模拟计算域初始布置示意图

Fig. 4  Schematic of the initial setup of the case of shock wave propagation

2.1　冲击波传播

炸药转化形成高压气泡向各个方向均匀膨胀，辐射产生冲击波。计算结果如图5所示，图中显示冲击波沿炸药爆轰的径向方向传播，压力传播速度极快，同时，炸药附近压力快速降低，并恢复到静水压强（图5d）。观察对称轴（云图中心）处，在不同时刻，流场压力均分布均匀，未出现压力物理性振荡现象。

图5  不同时刻水下爆炸冲击波压力云图

Fig. 5  Snapshots of pressure wave propagation at different times

在计算效率方面，新型的轴对称SPH模型使得三维水下爆炸问题可以在二维平面实现模拟。本算例在台式计算机上使用8核并行Inter（R） Core（TM） i9-11900K CPU完成计算，粒子总数为54万，相较于相同的分辨率下传统的三维SPH计算模型，粒子总数减少将近350倍，计算效率大幅提升，便于实现爆炸装药附近的精细流场模拟。

2.2　水下爆炸载荷验证

距爆心$\mathrm{1}\mathrm{m}$处的压力载荷变化曲线如图6所示。在波前阶段，压力急速增大并在极短时间内跃升到峰值，比较压力载荷峰值，可见轴对称SPH计算结果与文献［5］理论解及经验公式拟合值^［24］吻合良好；在冲击波后阶段，压力按指数规律衰减，变化趋势与经验公式拟合曲线吻合较佳，再次说明本文采用的新型轴对称SPH计算模型具有较高精度，能实现冲击波传播过程中压力载荷的精准预报。

图6  在爆距1 m处冲击波载荷对比

Fig. 6  Pressure pulse caused by shockwave propagation at a standoff distance of 1.0 m

3 双气泡耦合作用

工程中多发/多处炸药同时爆炸、气枪阵列等问题都涉及到多气泡的耦合作用，对此，本文选取典型的双气泡耦合进行探究。双气泡轴对称初始计算域布置见图7。定义量纲一距离参数$\gamma =h_{\mathrm{m}}/r_{\mathrm{b}}$表征双气泡之间的距离。文献［25］表明，在1.10<γ<2.00与 2.00<γ<3.65时气泡行为表现为不同模式，为结合文献试验结果验证本文新型轴对称模型在气泡作用模拟上的精度，并同时探究双气泡耦合作用，本文分别选取上述两个距离参数范围内$\gamma =\mathrm{1.66}$和$\gamma =\mathrm{3.01}$开展数值模拟研究。

图7  双气泡轴对称计算域初始布置示意图

Fig. 7  Snapshots initial setup of two bubbles interacting cases

3.1　双气泡距离参数$\gamma$＝1.66

当选取$\gamma$= 1.66时，计算结果如图8所示。由于初始气泡距离较大，双气泡在膨胀过程中没有发生合并。对比试验结果，SPH模拟的气泡形态变化与试验结果^［25］吻合良好，只是略微有些时间差，该现象在文献［26］中也有所提及。由图8可知，在水下爆炸过程中，气球快速膨胀，体积增大，并在$t=\mathrm{1.44}\mathrm{m}\mathrm{s}$时达到峰值，此时气泡内部压力低于水体压强，气泡开始收缩（$t=\mathrm{2.56}\mathrm{m}\mathrm{s}$）。在不断收缩过程中两气泡之间的水体压力较小，使得气泡产生较大的压力梯度，气泡对逐渐往中间区域运动，最终接触并破裂，可见，SPH计算结果可准确预报双气泡耦合作用下气泡形态与其物理特性变化。此外。在试验中由于气泡收缩导致的射流使得气泡界面不清晰，而数值模拟结果的气泡形态特性更为清晰。

图8  $\gamma$ ＝ 1.66，不同时刻双气泡形态

Fig. 8  Comparison of results for two bubbles interacting at $\gamma$ ＝ 1.66

3.2　双气泡距离参数 $\gamma$＝3.01

由图9可得，在$\gamma$ = 3.01工况下，SPH能够清晰捕捉气泡界面，气泡形态的变化发展过程与试验^［25］吻合，验证了本模型对于气泡脉动的强可压缩问题模拟的精度。压力变化规律与$\gamma$ = 1.66工况算例类似，同样地，气泡在未收缩之前保持规则的球形。由于气泡的间距变大，水体的阻滞作用增强，气泡产生的压力梯度较小，使得双气泡向中间区域移动的速度变慢，从而在收缩阶段，气泡并没发生相互接触，各自收缩形成环形气泡，相对于工况$\gamma$ = 1.66，气泡在发生收缩时两气泡之间的距离增大。

图9  $\gamma$ ＝ 3.01，不同时刻双气泡形态

Fig 9  Comparison of results for two bubbles interacting at $\gamma$ ＝ 3.01

图10进一步讨论两种距离参数条件下气泡总体积$V_{\mathrm{b}\mathrm{u}\mathrm{b}\mathrm{b}\mathrm{l}\mathrm{e}}$变化与气泡内部压强的关系。在计算中，通过累加每一个气体粒子代表的环形体积$V_i$得到气泡总体积为

$$V_{\mathrm{b}\mathrm{u}\mathrm{b}\mathrm{b}\mathrm{l}\mathrm{e}}=\sum _{i\in \mathrm{g}\mathrm{a}\mathrm{s}}\mathrm{2}\mathrm{\pi }{r}_i{V}_i$$

（12）

图10  不同距离参数下气泡最大压力和总体积时历曲线

Fig 10  Time evolution of maximum pressure and total volume of two bubbles at γ ＝ 3.01 and γ ＝ 1.66

对比两种工况下气泡总体积和压力变化曲线进一步说明，气泡膨胀伊始，气泡中心压力急速减小，直到气泡收缩时，气泡内部压力逐渐增大。当气泡收缩到最小体积时，压力跃升到峰值，量级可达到兆帕级，可见气泡脉动的压力加载对于结构物的毁伤不容忽视。随后气泡体积再次增大，气泡开启另一周期压力脉动。两种参数下的压力与体积变化曲线见图10。由图10可知，大距离参数下，气泡先到达最大膨胀半径，并取得较大的压力峰值，说明两气泡中心的距离参数与气泡溃灭时压力加载具有直接联系，在大距离参数下，水体的阻滞作用增强，气泡之间的相互作用减弱。因此，在实际军事应用中，多发制导武器的位置分布与攻击效果密切相关。

4 结论

本文采用光滑粒子流体动力学方法，针对水下爆炸问题中的轴对称特性，引入轴对称SPH方法处理近场水下爆炸，实现在二维算法框架下对三维水下爆炸过程的准确预报，计算效率大幅提高。针对传统轴对称SPH模型的对称轴处理、粒子体积分布不均匀等难点，通过对称轴处粒子镜像处理、粒子体积自适应技术，构建了高稳定性、高效率和高精度的新型轴对称SPH水下爆炸模型。计算结果表明：

（1）本文建立的新型轴对称SPH模型能够实现水下爆炸过程冲击波传播模拟，经对比验证可知，SPH压力载荷预报结果与文献结果吻合良好。

（2）与试验结果对比可知，SPH能够清晰重现双气泡相互作用过程的形态变化，并能实现气泡内部压力监测。

（3）气泡收缩到最小体积时，压力载荷达到最大，且其大小与气泡之间的距离参数$\gamma$密切相关。

（4）本文构建的新型轴对称SPH方法可有效减少计算成本，并实现冲击波传播和气泡脉动过程的全过程预报。

作者贡献声明

黄晓婷：数据分析验证，数值模拟，论文写作和修改。

孙鹏楠：概念提出，模型建立，研究指导，论文修改。

彭玉祥：算法修正，论文修改。

刘念念：数值验证，论文修改。

熊自明：概念提出，研究指导，论文修改。

邱艳宇：概念提出，研究指导，论文修改。