0 引言

空化是水动力学中重要的现象，通常发生在高速航行的回转体、水翼、以及高速运转的螺旋桨叶片的表面低压区。空化相变生成大量气泡，构成空泡群。空泡群的非稳态演化特别是溃灭现象，不仅严重影响航行体水动力性能，产生噪声与振动，甚至能够剥蚀结构表面导致破坏。水中压力波与气泡相互作用是空泡群溃灭的基本形式，而压力波作用下气泡的运动响应研究是解决空泡群溃灭问题的核心基础。

针对压力波与气泡相互作用，Ding和Gracewski^[1]通过数值模拟研究了在较弱压力波作用下气泡的响应过程，由于表面张力相对较大，直径100 mm的气泡可以保持球状进行体积振荡。Abe等^[2]采用显微镜扩大观测也发现直径小于50 mm的气泡，在压力波作用下更容易保持球状振荡。但是，当气泡尺寸相对较大时，在压力波作用下，气泡会发生非球形演化，其中，生成射流是非球形演化的最显著特征。Dear等^[3-4]通过实验观测的手段，分析了压力波加载二维柱状圆形气泡动态响应，发现射流的方向与压力波的传播方向相同。Bourne和Field^[5-6]分析了压力波加载下不同形状界面的演化特征，指出射流的方向与界面的曲率有关。随后, 国内外学者利用实验观测和数值模拟手段分析了压力波与气泡耦合问题。在数值模拟研究中，为了提高压力波和气/液两相界面的捕捉，建立了许多创新的数值计算方法，如Shyue的半保守模型^[7-8]、高分辨率的波前追踪方法^[9]、基于自适应特征匹配^[10]、高精度算法^[11-12]、自由-拉格朗日方法^[13]等，分析了压力波在界面的波系变化、界面涡量分布和射流生成等关键过程。但是，射流生成的内在机理尚不十分明确。

从现象上看，射流的生成与RM不稳定现象类似，当压力波加载不同密度的物质界面时，由于界面处压力梯度与密度梯度不共线，导致斜压机制，进而界面获得加速度，界面上的扰动会逐渐演化发展^[14-16]。RM不稳定性是一类十分复杂的多尺度/强非线性问题，是高能量密度物理研究的主要内容之一，广泛应用于惯性约束聚变、尖端武器和天体物理等领域。目前在RM不稳定性的实验研究中，主要利用冲击波加载不同密度气体的界面。本研究中的液/气界面作为不同密度流体介质的一个特例，界面的演化现象与气/气界面的RM不稳定性相似，主要体现在，当压力波加载界面时，两者都涉及各种流体动力学现象的强烈耦合问题，包括压力波在界面处的反射、透射、绕射等波系变化，涡量的产生与分布。但是由于界面物理特性方面的差异，两者也有不同。因此，水中压力波作用气泡界面的演化在学术研究中有着重要的研究价值与应用背景。

本研究以水中压力波加载不同倾斜角度的二维椭圆形气泡为主要对象，通过数值模拟方法，分析压力波作用下气泡界面的非球形演化过程，获得了气泡倾斜角度对射流生成与发展的影响规律，探索了水中压力波作用下射流生成的内在力学机理。在数值模拟中，基于OpenFOAM开源程序的二次开发，利用可压缩的VOF方法和LES方法分别捕捉气液运动界面和流场细节信息。

1 数值模拟方法 1.1 控制方程

对于水中压力波作用下气泡的动态响应问题的研究，本文基于OpenFOAM开源软件，采用可压缩的VOF (Volume of Fluid)方法捕捉气液运动界面的演化，质量方程有如下形式：

$ {\frac{{\partial {\alpha _l}}}{{\partial t}} + \nabla \cdot ({\alpha _l}\mathit{\boldsymbol{u}}) = - \frac{{{\alpha _l}}}{{{\rho _l}}}\frac{{D{\rho _l}}}{{Dt}}} $

(1)

$ {\frac{{\partial {\alpha _g}}}{{\partial t}} + \nabla \cdot ({\alpha _g}\mathit{\boldsymbol{u}}) = - \frac{{{\alpha _g}}}{{{\rho _g}}}\frac{{D{\rho _g}}}{{Dt}}} $

(2)

其中α_l和α_g分别为水和空气的体积分数。

将式(1)和(2)相加可得：

$ \nabla \cdot \mathit{\boldsymbol{u}} = - \frac{{{\alpha _g}}}{{{\rho _g}}}\frac{{D{\rho _g}}}{{Dt}} - \frac{{{\alpha _l}}}{{{\rho _l}}}\frac{{D{\rho _l}}}{{Dt}} $

(3)

由于式(3)中的对流项中含有非零项▽· u，故无法直接进行求解。所以将方程(1)展开后代入式(3)可得：

$ \frac{{\partial {\alpha _l}}}{{\partial t}} + \nabla \cdot ({\alpha _l}\mathit{\boldsymbol{u}}) = {\alpha _g}{\alpha _l}\left( {\frac{1}{{{\rho _g}}}\frac{{D{\rho _g}}}{{Dt}} - \frac{1}{{{\rho _l}}}\frac{{D{\rho _l}}}{{Dt}}} \right) + {\alpha _l}\nabla \cdot \mathit{\boldsymbol{u}} $

(4)

为了使气液界面更加尖锐，利用了Weller提出的界面压缩算法，引入人工对流项抵消数值耗散带来的误差。因此可压缩VOF方法的质量方程如下：

$ \begin{array}{*{20}{l}} {\frac{{\partial {\alpha _l}}}{{\partial t}} + \nabla \cdot ({\alpha _l}\mathit{\boldsymbol{u}}) + \nabla \cdot ({\alpha _l}(1 - {\alpha _l}){\mathit{\boldsymbol{u}}_c})}\\ { = {\alpha _g}{\alpha _l}\left( {\frac{1}{{{\rho _g}}}\frac{{D{\rho _g}}}{{Dt}} - \frac{1}{{{\rho _l}}}\frac{{D{\rho _l}}}{{Dt}}} \right) + {\alpha _l}\nabla \cdot \mathit{\boldsymbol{u}}} \end{array} $

(5)

其中u_c为界面压缩速度，u_c=f(▽α/|▽α|)。

另外，可压缩VOF方法的动量方程形式和混合物方程可写为：

$ \frac{{\partial \rho \mathit{\boldsymbol{u}}}}{{\partial t}} + \nabla \cdot (\rho \mathit{\boldsymbol{uu}}) - \nabla \cdot \tau = - \nabla P + {F_\sigma } $

(6)

$ \begin{array}{*{20}{l}} {\frac{{\partial \rho h}}{{\partial t}} + \nabla \cdot (\rho h\mathit{\boldsymbol{u}}) + \nabla \cdot (\rho K\mathit{\boldsymbol{u}}) - \frac{{\partial p}}{{\partial t}}}\\ { = \nabla \cdot q - \nabla \cdot (\tau \cdot \mathit{\boldsymbol{u}})} \end{array} $

(7)

其中，u混合相速度，ρ混合相密度，τ黏性力，F_σ根据CSF模型给出，即F_σ=σκ▽α_l，σ表面张力系数，κ界面曲率，h焓，K动能，q热通量。

1.2 湍流模型

为了能够获得两相界面流场的细节，本文采用LES (Large Eddy Simulation)方法。基本原理是，先在选定区域内采用滤波操作，将涡分为大尺度涡和小尺度涡，对于大涡直接通过求解N-S方程得到，对于小尺度的涡通过引入亚格子模型，体现其对于大尺度涡的影响。在大涡模拟中，经过滤波操作以后的物理量以横杠标记：

$ \bar \varphi (x, t) = \int G (|x - {x^\prime }|)\varphi ({x^\prime }, t){\rm{d}}{V^\prime } $

(8)

这里G(|x-x′|)代表滤波函数，本文中采用盒式滤波：

$ \begin{array}{l} G|x - {x^\prime }| = \\ \left\{ {\begin{array}{*{20}{c}} {\frac{1}{{\Delta {x_1}\Delta {x_2}\Delta {x_3}}}}&{, |{x^\prime }_i - {x_i}| \le \frac{{\Delta {x_i}}}{2}, \;\:i = 1, 2, 3}\\ 0&{, |{x^\prime }_i - {x_i}| > \frac{{\Delta {x_i}}}{2}, \;\:i = 1, 2, 3} \end{array}} \right. \end{array} $

(9)

其中，Δx₁、Δx₂和Δx₃代表网格各个方向的尺寸。对于可压缩流动，为了避免滤波后产生的非线性应力项，使方程封闭，本文还需采用Favre滤波，操作后的物理量用波浪线表示。滤波后控制方程中产生的亚格子应力张量τ_SGS通过Boussinesq假设封闭：

$ {\tau _{{\rm{SGS}}}} = - 2{\mu _{{\rm{SGS}}}}\left( {{{\tilde S}_{ij}} - \frac{{{\delta _{ij}}}}{3}{{\tilde S}_{kk}}} \right) + \frac{2}{3}\bar \rho {k_{{\rm{SGS}}}}{\delta _{ij}} $

(10)

亚格子黏性μ_SGS和亚格子动能k_SGS采用一方程模型计算：

$ \begin{array}{l} \frac{{\partial \bar \rho {k_{{\rm{SGS}}}}}}{{\partial t}} + \nabla \cdot (\bar \rho {{\tilde u}_j}{k_{{\rm{SGS}}}}) = \\ \nabla \cdot \left[ {(\frac{{{\mu _{{\rm{SGS}}}}}}{{P{r_t}}} + \mu )\nabla {k_{{\rm{SGS}}}}} \right] + {P_{{\rm{SGS}}}} - {\varepsilon _{{\rm{SGS}}}} \end{array} $

(11)

其中ε_SGS为亚格子黏性耗散率。

1.3 计算域设置

图 1为本文计算域的示意图，建立了120 mm×180 mm×2 mm的长方体计算域，用来呈现二维几何形状。在长方体中心处设置了一个圆柱形气泡。气泡右侧设定了矩形高压区，用于模拟水中平面压力波，计算域的上下两个面设置为无滑移壁面。另外四个面定义为出口，水可以自由流入流出，通过采用消波边界条件来避免压力波的反射。坐标系的原点建立在长方体的中心，也就气泡中心，以长度方向为x轴，宽度方向为y轴，高度方向为z轴。计算开始后，压力波沿y轴方向传播。以压力波阵面恰好接触到气泡的瞬间为0时刻。

1.4 计算方法验证

计算域采用结构化网格划分，为了排除网格数密度的影响，进行了网格无关性验证。图 2展示了利用三种密度网格计算气泡溃灭过程的对比。高、中等及低密度网格数目分别为533万、420万和310万。从结果看，对于气泡界面的演化，不同密度网格的计算结果基本一致。下文数值模拟采用的高密度网格，网格数为533万。

为了验证本文计算方法的正确性与可行性，将数值模拟的结果与文献中的实验结果^[17]进行了对比。参考Bourne和Field研究，气泡初始直径是12 mm，在他们的研究中，平面压力波由于炸药冲击青铜板生成的，是一个矩形波，压力峰值是260 MPa。在本文数值模拟中，利用矩形的高压区来模拟水中平面压力波，再根据张阿漫等研究^[18]，经过大量数值计算后选取与实验结果最接近的矩形高压区的宽度，2 mm。图 3展示了气泡界面演化的实验观测与数值模拟的对比。每个图片的时间间隔是10 μs。压力波的传播方向是从下向上，加载气泡后，其下表面逐渐变平，然后向气泡中心凸起，形成高速射流，最终可以打穿气泡壁。总体来说，气泡的形态变化基本与实验观测相吻合，从而证明了数值模拟方法的可靠性。

2 结果与讨论

本文主要分析水中压力波加载不同倾斜角度的二维椭圆形气泡的界面演化规律。椭圆气泡的倾角定义为气泡长轴和压力波传播方向的夹角。在计算中，椭球形气泡的长轴半径设置为4 mm，短轴半径为2 mm，模拟水中平面压力波的矩形高压区的宽度为2 mm，初始温度设置为320 K。

首先研究了椭圆形气泡倾斜角度为0°的情况。图 4展示了压力波加载后气泡界面形态变化。水中压力波传播方向从右向左。与圆形气泡的情况(图 3所示)类似，在压力波传播后，与其先接触的气泡界面开始变形，从凹面变成平面，再逐渐向气泡中心运动，形成射流(主射流)，与压力波传播方向相同。随着气泡界面整体运动，在t=50 μs时，在气泡长轴的另一端也形成一个射流(次射流)，与射流1方向相反。在t=55 μs时，主射流击穿气泡壁，造成气泡分裂。分裂后的两部分在主射流和次射流的共同作用下，界面不断收缩，形成类似旋涡的结构(t=70 μs)。从图中可以看出，在水中压力波加载椭圆形气泡后，在气泡长轴两端各生成一个射流，与压力波传播方向相同的射流占主导作用，因此，射流角度为0°。其中射流角度为主射流运动方向与水中压力波传播方向的夹角。

图 5展示了倾斜角度为35°情况下气泡动态响应。在压力波传播后20 μs时，发现气泡界面明显变化的地方是椭圆形的长轴的一端，并不是与压力最先接触的地方。因此，发现气泡界面的运动除了和压力波的压强有关，还依赖于界面的曲率。在曲率大的地方，界面运动速度也快。由此可得，射流是在气泡长轴一端生成的，但是并不是指向气泡中心的。在压力波作用下，射流运动的方向与压力波传播方向的夹角为31.4°，因此，该工况下射流角度为31.4°。同时，在气泡长轴的另一端又生成了一个射流(次射流)。当两个射流相聚并撞击后，气泡分裂成两部分，界面继续变形，之后的运动变得没有规律。

图 6展示了倾斜角度为90°时气泡界面响应情况。在水中压力波传播过后，随着界面的运动，先接触压力波的气泡界面不再保持平滑，形成了一些小的扰动，而是在气泡长轴两端气泡界面出现了明显的向内收缩(t=30 μs)，形成两个射流。在该工况下，主射流和次射流的速度相同，对气泡界面变形的作用相同。他们的速度方向与压力波传播方向的夹角均为53.9°，因此，该工况下射流角度为53.9°。同时，运动前期生成的小扰动逐渐增大，形成了有两个波峰的扰动。但是，与射流相比，小扰动的幅值和速度都是非常小的。随着两个射流继续向气泡内运动，在t=55 μs时，上下的射流汇聚并碰撞，造成气泡分裂，两个波峰变成了一个。随后，主体气泡界面继续运动，形成了一个与压力波传播方向相反的水平射流。该射流最终穿透气泡壁。

为了分析压力波作用下气泡射流形成的内在机制，提取了不同时刻下的界面的涡量场信息，如图 7所示。图中的黑色实线表示气泡的轮廓，这里定义液相体积分数0.9的等值面为气泡边界。本文主要针对是二维几何，所以与x和y轴相比，z轴方向距离很小，所以气泡变形和涡量主要存在x-y平面上。图 7展示的是x-y平面上的涡量场。在t=25 μs时刻，气泡上部的射流是在一对涡量的作用下生成的，右边是逆时针方向的旋涡，左边是顺时针的。同时，界面处生成的小扰动也是在涡对作用下生成的。在t=55 μs时刻之后，两个射流相聚并穿透气泡壁后，因此，作用在两个射流上的两对涡量变成了一对(如t=55 μs)，上面旋涡的方向为逆时针，下面的为瞬时针。在这样的涡对作用下，气泡生成了与压力波传播方向相反的水平射流。随着射流在泡内运动，涡量强度逐渐减弱。

为了明确诱导射流和界面扰动的涡量生成机制，提取了气泡界面处的斜压场(如图 8所示)。图中的黑色实线表示气泡的轮廓，这里定义液相体积分数0.9的等值面为气泡边界。从图中可以看出，斜压分布与涡量分布相似，在压力波传播后，气泡界面处密度梯度指向形心，压强梯度与界面曲率垂直，所以两者在界面处不共线，进而诱导了斜压机制，产生旋涡。从而也揭示了气泡界面曲率大的地方，界面运动速度快。因为曲率越大的地方，压强梯度与密度梯度夹角越大，进而产生涡量的强度越大。

为了进一步揭示斜压机制在涡量中的主导作用，分析了涡量方程中其他分项的作用。涡量方程式如下所示：

$ \frac{{{\rm{D}}\omega }}{{{\rm{D}}t}} = \frac{{\nabla \rho \times \nabla P}}{{{\rho ^2}}} + (\omega \cdot \nabla )u - \omega (\nabla \cdot \mathit{\boldsymbol{u}}) + \nu {\kern 1pt} {\nabla ^2}\omega $

(12)

其中, 式中ω、u、ν、▽ρ和▽p分别表示涡量矢量、速度矢量、黏度系数、密度梯度和压力梯度。式子右边第一项就是斜压项，(▽ρ×▽p)/ρ²，表示由于密度梯度和压力梯度不共线时诱导的涡度。第二项涡拉伸项，(ω·▽)u，发生在速度场不均匀的情况，特别是在三维情况下发生在两相界面。第三项涡膨胀项，ω(▽·u)，主要是由于高压冲击强度而产生的压缩性。第四项黏性项，ν▽²ω，表示由于流动中的黏性效应而产生的涡度。为了定量对比上述作用项对涡量的影响，将式(13)在x、y和z方向上展开，可以得到：

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {\frac{{\partial {\omega _x}}}{{\partial t}}}\\ {\frac{{\partial {\omega _y}}}{{\partial t}}}\\ {\frac{{\partial {\omega _z}}}{{\partial t}}} \end{array}} \right] = \frac{1}{{{\rho ^2}}}\left[ {\begin{array}{*{20}{l}} {\frac{{\partial \rho }}{{\partial y}}\frac{{\partial P}}{{\partial z}} - \frac{{\partial \rho }}{{\partial z}}\frac{{\partial P}}{{\partial y}}}\\ {\frac{{\partial \rho }}{{\partial z}}\frac{{\partial P}}{{\partial x}} - \frac{{\partial \rho }}{{\partial x}}\frac{{\partial P}}{{\partial z}}}\\ {\frac{{\partial \rho }}{{\partial x}}\frac{{\partial P}}{{\partial y}} - \frac{{\partial \rho }}{{\partial y}}\frac{{\partial P}}{{\partial x}}} \end{array}} \right] + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {\frac{{\partial {\omega _x}}}{{\partial x}} + \frac{{\partial {\omega _y}}}{{\partial y}} + \frac{{\partial {\omega _z}}}{{\partial z}}} \right)\left[ {\begin{array}{*{20}{l}} u\\ v\\ w \end{array}} \right] - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\begin{array}{*{20}{l}} {{\omega _x}}\\ {{\omega _y}}\\ {{\omega _z}} \end{array}} \right]\left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}} + \frac{{\partial w}}{{\partial z}}} \right) + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \nu \left[ {\begin{array}{*{20}{c}} {\frac{\partial }{{\partial x}}\left( {\frac{{\partial {\omega _x}}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {\frac{{\partial {\omega _x}}}{{\partial y}}} \right) + \frac{\partial }{{\partial z}}\left( {\frac{{\partial {\omega _x}}}{{\partial z}}} \right)}\\ {\frac{\partial }{{\partial x}}\left( {\frac{{\partial {\omega _y}}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {\frac{{\partial {\omega _y}}}{{\partial y}}} \right) + \frac{\partial }{{\partial z}}\left( {\frac{{\partial {\omega _y}}}{{\partial z}}} \right)}\\ {\frac{\partial }{{\partial x}}\left( {\frac{{\partial {\omega _z}}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {\frac{{\partial {\omega _z}}}{{\partial y}}} \right) + \frac{\partial }{{\partial z}}\left( {\frac{{\partial {\omega _z}}}{{\partial z}}} \right)} \end{array}} \right] \end{array} $

(13)

在这里，主要讨论x-y平面上的气泡变形和涡量。

图 9展示了在气泡响应过程中涡量方程式中各项云图的对比。时间分别为t=20 μs、40 μs、50 μs和65 μs。在图中，涡量方程式中各项颜色条的范围是不同的，前三项，斜压项、涡拉伸项和涡膨胀项的范围是-5×10⁷~5×10⁷，而黏性项是-5×10⁶~5× 10⁶。由此可以看出，在气泡响应过程中，黏性对涡量演化的影响最小。从图中前三项的云图中，可以看出斜压项的作用最大，第二项涡拉伸项和第三项涡膨胀项的差不多，但是根据涡量方程式，第二项和第三项方向相反。由此可见，斜压项是占主导作用的，尤其是在气泡响应后期阶段(t=65 μs)。因此，可以得出气泡射流的生成是由于RM不稳定性中斜压机制导致的，压力波传播后，界面处的压力场梯度与密度梯度不共线引起的。

图 10展示了倾斜角度等于145°情况下，气泡界面演化。相比倾斜角度35°工况(图 5所示)结果，两者的气泡界面演化在垂直面上对称。在压力波扫过气泡后，在气泡长轴上端点位置开始出现明显向内收缩，生成主射流，其运动方向与压力波传播方向夹角为31.4°，因此，该工况下射流角度为31.4°。

最后，分析了不同气泡倾斜角度对射流角度的影响，如图所示11。横轴是气泡倾斜角度，纵轴是射流角度。气泡倾斜角度为长轴与压力波传播方向之间的夹角，射流角度为主射流的运动方向与压力波传播方向之间的夹角。当气泡长轴与压力波传播方向平行的情况，射流运动方向与其相同。当气泡长轴垂直于压力波传播方向的情况，射流的角度等于53.9°。气泡倾斜角度在0°~90°的范围，射流角度随倾斜角度增大而增加。在90°~180°之间范围，射流角度随倾斜角度增大而减小。而且，从图中还发现，射流角度在这两个范围内是对称的。

3 结论

本文基于OpenFOAM开源程序的二次开发，利用可压缩的VOF方法和LES方法分别捕捉气/液运动界面和流场细节信息。以水中平面压力波加载不同倾斜角度的二维椭圆形性气泡为主要研究对象，分析了气泡界面演化规律，揭示了射流生成的内在机理。获得的主要结论如下：

1) 射流生成位置与气泡倾斜角度没有关系，都是从气泡长轴两个端点上生成的，然后向气泡内运动。这可能与界面曲率有关，长轴端点的界面曲率大，因而界面运动速度快；

2) 射流角度与气泡倾斜角度密切相关。当气泡倾斜角度等于0°的情况，主射流的运动方向与压力波传播方向相同，当气倾斜角度等于90°的情况，射流的运动方向与压力波传播方向夹角为53.9°，当倾斜角度在0°~90°之间，射流角度随着气泡倾斜角度增大而增大。而且还发现，射流角度在倾斜角度90°~180°的工况与0~90°的工况是对称的；

3) 通过定量对比涡量方程中斜压项与其他项作用，发现斜压是诱导涡量产生的主导机制。当压力波扫过气泡后，气泡界面处的压力梯度和密度梯度不共线，诱导了斜压机制，进而生成射流和界面的扰动。另一方面也解释了界面曲率大的地方运动速度快，曲率越大，界面处的压强梯度和密度梯度方向差别越大，斜压项作用越大，从而诱导的涡量强度越大。