在船舶与海洋工程领域，爆炸^[1-6]问题具有极其重要的研究意义。随着水下高技术武器的发展，舰船和潜艇等海洋装备的安全受到越来越严峻的挑战，水下爆炸^[7-11]对海洋装备造成的结构毁伤十分严重。因此，研究舰艇在水下爆炸作用下的毁伤特性十分必要，掌握不同边界条件下的舰艇毁伤机理对提高现代海军的作战能力具有重要意义。

水下爆炸载荷作用较为复杂，主要分为瞬态的冲击波载荷和脉动的气泡载荷。对于近场水下爆炸，高强度的瞬态冲击波会造成结构的塑性破坏，气泡的脉动将会进一步造成舱室涌流以及结构的全局破坏。从水下爆炸载荷的作用机理出发，冲击波过后会形成气泡和不完整边界的耦合作用，科学家们在该问题上进行了很多研究。刘云龙等^[12]对背空船底板在气泡载荷下的破坏机理进行了研究，刘念念等^[13]对有限水域破口问题进行了研究，刘润泉等^[14]进行了舰船单元结构模型水下接触爆炸破口的实验研究，李金河等^[15]通过实验分析了近场冲击波的传播规律，杨棣等^[16]对接触爆炸载荷作用下船体板架破口大小进行了预测。对于复杂边界条件下的水下爆炸问题来说，理论解的获得是极其困难的，而实验研究^{[11, 14-15]}又会花费大量的人力物力，所以介绍实尺寸实验的文献特别少。随着计算机技术的不断发展，数值模拟方法^[17-18]逐渐成为研究水下爆炸破坏机理及舰艇结构抗冲击性能的有效工具，其中边界积分方法^[19]、有限体积方法^[20]和耦合欧拉-拉格朗日方法^[21]等广泛应用于水下爆炸问题的模拟。

本文中基于欧拉有限元方法对近场水下爆炸气泡与双层破口结构的相互作用机理进行研究。通过实验验证欧拉有限元理论数值模型的准确性，通过气泡形态、流场压力、舱室进水量等参数对不同破口尺寸、不同起爆位置和不同壳间水位条件下的水下爆炸过程进行系统分析。

1.   问题描述

对于装备双层外壳的舰艇结构，当冲击波载荷将双层外壳破坏后，其内部结构会受到因气泡和复杂不完整边界的耦合作用产生的水射流和涌流的冲击，而由于破口结构的卸载作用，涌流效应更加明显，高速涌流会对内层壳结构造成毁伤，其作用过程本质上是近场水下爆炸气泡与双层破口结构的相互作用问题，如图1所示。

图  1  近场水下爆炸气泡与双层破口结构相互作用示意图

Figure  1.  Schematic of interaction between bubble and double breaken structure in near field underwater explosion

下载: 全尺寸图片 幻灯片

2.   理论数值模型

2.1   可压缩欧拉有限元

为了在保证求解精度的前提下提高数值模拟的求解效率，将复杂的工程问题简化为轴对称的数值模型。可压缩欧拉有限元方法^[22-24]在处理流体瞬态大变形问题中具有明显的优势，如图2所示。忽略黏性项，在求解过程中每个单元的控制方程如下：

图  2  欧拉有限元方法示意图

Figure  2.  Schematic of Eulerian finite element method

下载: 全尺寸图片 幻灯片

式中：${{u}}=\{ {u_r},{u_z}\}$为速度矢量，$\varphi$为待定变量，$S$为源项。式(1)～(2)是在柱坐标系下建立的，$r$、z分别表示径向和轴向。式(1)～(2)通过算子分离的方式分成拉格朗日步和欧拉步，这样的处理方式结合了拉格朗日描述和欧拉描述的优势。对于质量、动量和能量方程来说，$S$和$\varphi$表现为不同的形式：

式中：$\rho$为流体的密度，g为重力加速度,${e_{\rm{in}}}$为质量内能，${p'_i}$表示在$i$方向上的偏导数，速度散度$\nabla \cdot {{u}}$=$\displaystyle\sum\limits_{i=1}^2 {\dfrac{{\partial {{{u}}_i}}}{{\partial {x_i}}}} + \frac{1}{r}{u_r}$，下标$i=1,2$分别表示径向r和轴向z。

拉格朗日计算步是对方程(1)的求解，本文中采用显式有限元方法处理。在柱坐标下，对动量方程进行积分及分步变换，可得：

式中：$\varOmega$表示离散的二维正交计算域，$\varGamma$表示计算域的边界，${ n}$表示指向计算域外侧的法向量，$\phi$表示形函数。式（4）可用于离散求解。文中采用显式积分格式更新流体速度，并且网格随流体材料的改变而改变，即:

式中：上标表示时间步数，$\Delta t$表示时间步长。系统的材料密度和能量方程通过下式求解：

此时每个计算时间步的拉格朗日过程处理完毕，欧拉的计算步骤开始执行。本文中欧拉有限元主要通过体积分数法和单调迎风输运格式确定单元之间的输运量，进而完成动量等变量的更新。下一个计算步的初始压力等变量则通过状态方程进行更新。拉格朗日计算步和欧拉计算的反复迭代，实现了复杂问题的求解。另外，本文中选择Tammann方程^[25]作为水和空气的状态方程，即：

式中：$\chi$表示绝热率，${p_{\rm{w}}}$表示参考压力。对于水，${p_{\rm{w}}}$=330.9 MPa，$\chi$=7.15，$\rho$=1 000 kg/m³；对于空气，${p_{\rm{w}}}$=0 MPa、$\chi$=1.25、$\rho$=1.25 kg/m³。

另外，爆炸产物的状态方程为Jones-Wilkens-Lee (JWL)方程^[26]，即：

式中：A=371.2 GPa, B=3.2 GPa, ${R_1}$=4.15, ${R_2}$=0.95，w=0.3，是TNT炸药爆轰产物参数；${\rho _0}$=1 630 kg/m³，为初始炸药密度。

2.2   相关参数与无量纲处理

为了对物理问题进行深层次的理解和认识，并拓展得到的结论和规律，对相关参数进行无量纲化：无量纲时间${t^{\rm{*}}}{\rm{=}}t/\left({R_{\rm{m}}}\sqrt {{{{\rho _{\rm{l}}}}}/{{{p_\infty }}}} \right)$，其中${\rho _{\rm{l}}}$为水的密度，${p_\infty }$为无穷远处的压力，${R_{\rm{m}}}$为气泡最大半径（通过${R_{\rm{m}}}=3.38{\left(\dfrac{W}{{10 + d}}\right)^{{1}/{3}}}$求解^[23]）；无量纲破口尺寸${R_{\rm{h}}^{\rm{*}}}{\rm{=}}{R_{\rm{hi}}}/{R_{\rm{b}}}$，其中${R_{\rm{hi}}}$为破口尺寸（${R_{{\rm{h}}1}}$为内板破口、${R_{{\rm{h}}{\rm{2}}}}$为外板破口），${R_{\rm{b}}}$为两层板之间的距离（文中为1.0 m）；无量纲起爆距离${R_{\rm{d}}^{\rm{*}}}{\rm{=}}{R_{\rm{d}}}/{R_{\rm{m}}}$，其中${R_{\rm{d}}}$为起爆距离；无量纲涌流量${{V_{\rm{OL}}^{\rm{*}}}}{\rm{=}}{V_{\rm{OL}}}/{R_{\rm{m}}^3}$，其中${V_{\rm{OL}}}$为涌流量（${V_{{\rm{OL}}{\text{总}}}}$为总涌流量、${V_{{\rm{OL}}{\rm{1}}}}$为内舱室涌流量）；无量纲压力${p^{\rm{*}}}{\rm{=}}p/{p_\infty }$；无量纲压载水系数${R_{\rm{s}}^{\rm{*}}}{\rm{=}}{R_{\rm{s}}}/{R_{\rm{b}}}$，其中${R_{\rm{s}}}$为壳间的水位。

3.   实验验证

利用气泡放电实验对所建立的水下爆炸数值模型进行验证。针对所研究的问题，本文选择气泡与自由液面及破口平板相互作用的实验工况进行对比验证。实验在50 cm×50 cm×50 cm的水箱内完成，其中气泡由220 V高压放电产生，高速摄像机用于拍摄整个实验过程，具体的实验操作可参考文献^[11]，实验装置如图3所示。另外，带有30.0 mm破口的2.0 mm厚平板放置于静止水面，气泡初始位置距离水面12.0 mm，气泡最大半径12.6 mm。

图  3  实验装置图

Figure  3.  Picture of experimental device

下载: 全尺寸图片 幻灯片

图4为不同时刻的数值结果与实验结果对比图。

图  4  数值结果和实验结果的对比图

Figure  4.  The comparison between numerical and experimental results

下载: 全尺寸图片 幻灯片

在数值模拟中，计算模型与实验模型中破口尺寸、板厚等基本参数一致。另外，计算域尺寸为0.15 m×0.05 m，最小网格尺寸为2.5×10⁻⁴ m，计算域边界为数值的无反射边界^[27]。由图4可知，数值模拟中的气泡膨胀、收缩和射流等气泡形态和水冢现象与实验现象基本一致。其中，由气泡能量传递所形成的自由面处水冢并诱导结构产生的涌流现象被有效地模拟，其多峰效应也被很好地展示。另外，如图4(d)所示，本文的数值模拟方法能较好地模拟气泡运动后期的撕裂破碎和破口附近复杂的自由液面演化等过程。由图4可知，量化的涌流高度及气泡形态典型参数基本一致。综上所述，本文的数值结果与实验结果吻合较好，整个气泡运动过程都被很好地模拟出来，从而验证了算法的有效性和准确性。下文中将讨论近场水下爆炸气泡与双层破口结构的相互作用机理。

4.   结果讨论

4.1   基本物理现象

近场水下爆炸气泡与双层破口结构的相互作用是一个复杂的瞬态非线性问题，影响该物理问题的参数较多。本文中以${R_{\rm{m}}}=5.0\;{\rm{ m}}$、${R_{{\rm{h}}{\rm{1}}}^*}{\rm{=3}}{\rm{.0}}$、${R_{{\rm{h}}{\rm{2}}}^*}=1.5$和${R_{\rm{d}}^*}=0.4$为参考工况分析该物理过程的相关现象，为下文分析多种因素影响提供对比。主要研究该物理过程中流场的动态变化，所以假设破口边界为刚性边界。图5为不同时刻的气泡的形态和流场压力变化图。

图  5  不同时刻气泡的形态和流场压力变化图

Figure  5.  The bubble shape and flow field pressure change at different moments

下载: 全尺寸图片 幻灯片

由于气泡初始内压远大于周围流场压力，因此气泡会迅速膨胀。当$t=0.{\rm{027 }}\; {\rm{s}}$ 时，随着自由液面进一步运动，涌流会超过双层内壳继续向内舱室涌入，破口附近空气流动较剧烈，流速可以达到百米每秒量级。当$t=0.{\rm{037 }}\; {\rm{s}}$ 时，可以明显地看到破口对气泡形态的影响,气泡形态开始产生畸变。另外，由于破口的形成，气泡的运动会带动自由液面向上运动，进而形成舱室气泡涌流现象。在舱室进水的过程中，空气压缩导致了舱室内空气的回流。在内部空气、流体惯性以及破口诱导的联合作用下，当$t=0.{\rm{062 }}\; {\rm{s}}$ 时，气泡演化过程中出现了气泡分割现象，初始的爆炸气泡分成了几个小气泡和一个较大的主气泡。在自由液面和破口结构的诱导作用下，当$t=0.{\rm{090 }}\; {\rm{s}}$ 时，在气泡收缩过程中产生了一个远离自由液面的细长射流，射流穿透气泡形成环状气泡，并进入回弹阶段。当$t=0.{\rm{309 }}\; {\rm{s}}$ 时，由于流体质点的羽化运动，结构内产生液面飞溅现象。考虑舰艇的生命力，大量的水会随着气泡的脉动涌入到破口中，导致结构加速沉没。

4.2   内外板破口尺寸的影响

水下爆炸瞬态冲击载荷和结构材料属性不同，结构破坏后形成的破口尺寸也不同，破口尺寸不同会改变物理模型的耦合边界，对舱室涌流量和气泡运动状态产生较大的影响。另外，根据冲击波的作用机理，内层板和外层板的破口尺寸也会有所不同。在保持${R_{{\rm{h}}{\rm{1}}}^*}{\rm{=}}3.0$不变的情况下，选择${R_{{\rm{h}}{\rm{2}}}^*}{\rm{=}}0.5$，${\rm{1}}{\rm{.0}}$，${\rm{1}}{\rm{.5}}$ 时分别进行数值模拟。图6为${R_{{\rm{h}}{\rm{2}}}^*}={\rm{0}}{\rm{.5}}$时不同时刻气泡的形态和流场压力变化图。

图  6  不同时刻气泡的形态和流场压力变化图

Figure  6.  The bubble shape and flow field pressure change at different moments

下载: 全尺寸图片 幻灯片

由图6可知，当保持外板破口尺寸不变、只改变内板破口尺寸时，初始阶段气泡及流场运动特性区别不大。而当涌流与内板接触形成新的边界条件后，气泡和流场的运动特性发生明显改变。相比4.1节中的工况，${R_{{\rm{h}}{\rm{2}}}^*}={\rm{0}}{\rm{.5}}$ 时涌流的形态较稳定。当$t=0.{\rm{027 }}\; {\rm{s}}$ 时，涌流开始穿过内板破口进入舱室内，涌流尺寸和破口尺寸近似相等，此时自由液面与破口结构之间会形成小的液膜，这是由于空气压缩导致的。由流场速度变化图可以看到，此时内舱室中出现一个较大的高流速区，这是内舱室破口抑制空气的扩散导致的。当$t=0.{\rm{070 }}\; {\rm{s}}$ 时，内舱室破口同时抑制了水的局部流动，气泡产生了宽射流现象。当$t=0.{\rm{592 }}\; {\rm{s}}$ 时，涌流与内板接触形成新的边界条件，诱导分离后的气泡产生进一步破碎，同时内舱室内形成了二次涌流，涌流的形状呈现“纺锤”状。总之，破口尺寸对整个流场演化影响很大。

图7中对比了三种不同破口尺寸下，总涌流量和内舱室涌流量随时间的变化情况。由图7可知，舱室涌流大致可以分成以下主要阶段：初始瞬间急速涌流阶段、涌流稳定阶段和滞后回落阶段。其中，初始急速涌流阶段是气泡膨胀带动流体的瞬态运动，稳定阶段是舱室内的液体在多种因素诱导下达到一个相对动态的平衡状态，水流回落阶段是流体受重力作用后的向下运动。当保证结构外底板尺寸不变时，不同工况下舰艇的总涌流量以及水进入内舱室所花费的时间基本一致。另外，在初始内板破口参数${R_{{\rm{h2}}}^*}{\rm{ {\text{＜}} 1}}{\rm{.5}}$的范围内，随着破口尺寸的增大，舱室总涌流量减小，而最终的内舱室涌流量却变大，这是诱导舰船和潜艇发生沉没现象的重要原因。由此可知，当破口产生后，破口尺寸越大其受涌流进水的影响越明显，提高内层壳的结构强度能充分地提高舰艇抗冲击的能力。

图  7  不同破口尺寸下舱室涌流量随时间的变化曲线

Figure  7.  Changes of cabin inrush flow with time under different breach size

下载: 全尺寸图片 幻灯片

4.3   炸药起爆初始位置的影响

水中武器对舰艇进行攻击时，炸药的起爆位置并不相同。起爆位置距离结构越近，气泡和复杂边界的耦合作用越明显，气泡后期运动及涌流特性越复杂。本文中选择${R_{\rm{d}}^*}=0.{\rm{1}}$，${\rm{0}}{\rm{.2}}$，${\rm{0}}{\rm{.4}}$ 分别进行数值模拟，得出炸药起爆初始位置对该物理过程的影响。图8为${R_{\rm{d}}^*}=0.{\rm{1}}$ 时不同时刻气泡的形态和流场压力变化图。

图  8  不同时刻气泡的形态和流场压力变化图

Figure  8.  The bubble shape and flow field pressure change at different moments

下载: 全尺寸图片 幻灯片

由图8可知，炸药起爆位置距离结构较近时，自由液面处出现了破碎和重闭合现象。针对${R_{\rm{d}}^*}=0.1$ 的工况，在初始气泡膨胀阶段，当$t=0.{\rm{009 }}\; {\rm{s}}$ 时，气泡迅速涌入到结构内侧，气泡内气体运动速度可以达到500 m/s。破口的存在对气泡的演化具有明显的诱导作用，在双层结构的诱导下，当$t=0.{\rm{046 }}\; {\rm{s}}$ 时，自由液面出现分流趋势，上表面在内板的作用下重新闭合。当$t=0.{\rm{075 }}\; {\rm{s}}$ 时，由于流体运动较为剧烈，气泡的羽化破碎和液滴飞溅现象会出现在双层结构壳内部，这种较为精细的数值模拟用传统的数值模型是较难实现的。同时，自由表面的相对运动使气泡产生了向下的“皇冠”型射流。当$t=0.{\rm{227 }}\; {\rm{s}}$ 时，在气泡回弹阶段，气泡进一步破碎形成一些无规则的小气泡，这也是由于局部流体速度不均匀造成的。相比于4.1中的工况，涌流形态整体上呈现“飞溅”型，这使得随着时间推移只有少量流体会涌入到内舱室。

图9为不同起爆位置下，舱室涌流量随时间的变化曲线。由图9可知，在初始位置${R_{\rm{d}}^*} {\text{＜}} {\rm{0}}{\rm{.4}}$ 的范围内，随着起爆位置与自由液面位置距离的减小，初始瞬间涌流的速度变大（即曲线的斜率较大），水进入舱室所花费的时间越短。而爆炸气泡脉动诱导的进水量最大值和最终舱室的涌流量却越小，这是由于初始时刻在气泡上方的水较少造成的。由此可知，当破口产生后，爆炸产物距离结构越近，其受涌流进水的影响反而越小。由于双层结构之间间隔的存在，初始阶段内舱室的进水量所占的比例较小，水不能完全进入内舱室，进而存在一定的反应时间，这也有利于舰船和潜艇的紧急应急，避免结构沉没。

图  9  不同起爆位置下舱室涌流量随时间的变化曲线

Figure  9.  Changes of cabin inrush flow with time under different detonation position

下载: 全尺寸图片 幻灯片

4.4   双层壳间初始水位的影响

舰艇双层壳之间的水密舱可以用来调节结构自身重力，从而改变其在海里的平衡浮力，由于自由液面的存在，水密舱内水位对本文所涉及的多相流动问题存在较大影响。因此应分析水密舱内水位对近场水下爆炸气泡与双层破口结构相互作用的影响。文中选择${R_{\rm{s}}^*}=0$，$0.5$，$1.0$分别进行数值模拟，得出双层壳间初始水位对该物理问题的影响。图10为${R_{\rm{s}}^*}={\rm{1}}{\rm{.0}}$ 时不同时刻气泡的形态和流场压力变化图。

图  10  不同时刻气泡的形态和流场压力变化图

Figure  10.  The bubble shape and flow field pressure change at different moments

下载: 全尺寸图片 幻灯片

由图10可知，当舱室壳内水位较高时，在静水压作用下，整个流场流体特性相对较为稳定。在气泡膨胀阶段，舱室壳内的水就开始涌入内舱室，形成舱室涌流。由于自由液面处水冢速度较大（流速可以达到50 m/s），所以水冢高度不断升高并带动周围空气产生无规则运动。伴随着水冢的演化，破口附近液面逐渐收缩，这是由于局部压力不平衡和涌入双层壳结构的水与进入舱室内的水的体积差造成的。当$t=0.267{\rm{}}\; {\rm{s}}$ 时，气泡受浮力作用较明显，从而产生向下和向上的对射流，其中向上的射流是由于局部较大的静水液压导致的。随着气泡的演化，当$t=0.{\rm{572 }}\; {\rm{s}}$ 时，气泡穿透进入环状阶段，并进一步破碎形成多个气泡环。整体上来说，此时破口对气泡运动影响不明显，而舱室涌流特性变化较剧烈，接下来将对其进行定量分析。

图11为不同水位位置下，舱室涌流量随时间的变化曲线。由图11可知，当${R_{\rm{s}}^*}={\rm{1}}{\rm{.0}}$时（即液压舱满舱），涌流量较大，涌流量迅速增加直至平衡状态，也就是说，此时舰艇结构极易沉没。而当${R_{\rm{s}}^*} {\text{＜}} 0.5$时，随着${R_{\rm{s}}^*}$的增大，涌流速度逐渐减小，内舱室进水时间逐渐增加，导致内舱室的水量增大。当${R_{\rm{s}}^*}$较小时，水域流场静水液压较小，其在气泡脉动后期较易受到能量传递作用而产生向上的运动，导致后期涌流增加。总之，壳内水位对舱室涌流量的影响作用较复杂，具有较强的非线性，适当对其进行调整可以有效减轻水下爆炸涌流对结构的破坏。

图  11  不同水位下舱室涌流量随时间的变化曲线

Figure  11.  Changes of cabin inrush flow with time at different water levels

下载: 全尺寸图片 幻灯片

5.   结　论

基于欧拉有限元方法对近场水下爆炸气泡与双层破口结构的相互作用机理进行了研究，对不同破口尺寸，不同起爆位置条件下的不同工况进行了较详细的对比和分析。具体结论如下：

（1）当${R_{\rm{d}}^*}=0.4$、${R_{{\rm{h}}{\rm{1}}}^*}{\rm{=3}}{\rm{.0}}$、${R_{{\rm{h}}{\rm{2}}}^*}=1.5$时，在内部空气、流体惯性以及破口诱导的联合作用下，气泡演化过程中出现了气泡分割现象，初始的爆炸气泡分成了几个小气泡和一个较大的主气泡，多个羽化气泡继续演化。

（2）内层破口尺寸越小，内板结构对涌流的抑制作用越大。当${R_{{\rm{h}}{\rm{2}}}^*}=0.5$时，内舱室内会出现二次涌流现象，且涌流形态较细长。在初始内板破口参数${R_{{\rm{h2}}}^*}{\rm{ {\text{＜}}1}}{\rm{.5}}$的范围内，随着破口尺寸的增大，内舱室的涌流量增大。因此加强内层板的强度可以有效提高舰艇生命力。

（3）炸药起爆位置距离结构较近时，自由液面处会出现破碎和重闭合现象。在初始位置${R_{\rm{d}}^*}{\text{＜}} {\rm{0}}{\rm{.4}}$的范围内，随着起爆位置与自由液面位置距离的减小，初始瞬间涌流的速度变大（即曲线的斜率较大），水进入舱室所花费的时间越短，爆炸气泡脉动诱导的进水量最大值和最终舱室的涌流量却越小。

（4）当舱室壳内水位较高时，在静水压作用下，整个流场流体特性较稳定。当${R_{\rm{s}}^*}{\text{＜}} 0.5$时，随着${R_{\rm{s}}^*}$增大，涌流速度逐渐减小，内舱室进水时间逐渐延长，导致舱室的水量增大。当水位满舱时，急速涌流会缩短船艇的应急时间。所以适当调节舱内水位也可以提高舰艇生命力。