固壁的表面形状影响空泡的生长和破裂^[1]，米建东等^[2-4]开展了针对平坦壁面附近的空泡动力学研究. 在实际工程中，空化多发生在不规则壁面，形成的空泡在内外泡壁压差的作用下发生溃灭的过程被称为空泡溃灭，空泡溃灭产生的冲击致使壁面材料损失的现象被称为空蚀. 在壁面材料受到空蚀破坏形成凹坑后，空泡附近的流场受到新形成的表面几何形状的影响，相较于单一平坦壁面，这种壁面对空泡特性的影响非常强烈.

关于凹坑壁面附近空泡破裂过程的研究很少，空泡动力学理论研究多为对Rayleigh-Plesset方程的修正，空泡溃灭多发生在近壁区，具有非球形、非线性和瞬时性等特点，往往无法求得解析解. Zhang等^[5]基于双渐进分析方法(double progressive analysis method)，引入三维波动方程，将双渐进法扩展为非球形空泡的边界积分法，研究了弱可压缩流场下的空泡运动. Zhang等^[6]建立新的振荡气泡动力学理论，考虑了边界条件、气泡相互作用、流场可压缩性等影响，统一了经典气泡方程. 固壁^[7]和自由表面^[8]是研究较多的两类边界，由于空泡与壁面的相互作用，空泡在近壁面的溃灭过程相比在自由场中的溃灭过程更为复杂. Cui等^[9]使用电火花生成气泡和高速摄像技术(high speed video technology)研究具有圆形开口舱底附近的水下爆炸气泡的振荡特性，发现气泡的演化过程高度依赖于开口距离和开口半径；从180个实验中总结出7种独特的振荡场景，即边缘约束振荡、向内射流、向外射流、凸起和凹陷、准球形振荡、有射流的球形振荡和无射流的球形振荡. Kim等^[10]对有脊状和凹槽图案结构附近的水下空泡溃灭过程开展实验研究，发现脊上的气泡在收缩过程中损失了一些能量，削弱了向表面射流的强度，有效减少了对结构表面的侵蚀. Ma等^[11]以半球为主要研究边界，通过使用400 V的高压电产生电火花空泡研究3种半球形边界与空泡运动特性之间的关系，发现在同一曲率下，空泡与壁面距离越大，壁面受到的冲击载荷越大；在同一壁面距离下，曲率越大，空泡的射流速度越小. Aziz等^[12]使用边界积分法(boundary integral method)对弯曲刚性板间的三维气泡振荡进行数值模拟，发现随着空泡溃灭程度的增加，空泡沿径向拉伸长度逐渐增大并最终分裂为2个空泡，空泡射流速度随偏心距的增大而增大. Manmi等^[13]采用边界积分法研究处于声场中的两弯曲刚性板间的微气泡，发现空泡最大半径及射流速度不受声场影响，空泡破裂产生的高压区与射流方向受声场影响较大.

上述研究多以半球形壁面形式分析空泡在不同曲率壁面下的运动特性，揭示了空泡在不同曲率壁面下射流速度、压力载荷和空泡形态等变化规律. 本研究以单空泡为对象，基于欧拉方法的可压缩两相流求解器CavitatingFoam，给定多种凹坑曲率，探究空泡在不同曲率凹坑壁面附近的空泡溃灭特性和载荷特性.

选择基于欧拉方法的CavitatingFoam求解空泡模型^[14]. 该求解器采用压力和密度的关系封闭方程，通过压力修正密度. 假设计算流体为两相均质可压缩混合物，各向同性的正压关系^[15]为

(1)$ \frac{{{\text{d}}{\rho _m}}}{{{\text{d}}t}} = {\psi _{\text{m}}}\frac{{{\text{d}}p}}{{{\text{d}}t}}. $

式中：ρ_m为混合物密度；p为压力；t为时间；ψ_m为混合物的压缩系数，

(2)$ {\psi _{\text{m}}} = {\rho _{\text{m}}}{\beta _{\text{s}}}. $

β_s为等熵压缩率，

(3)$ {\beta _{\text{s}}} = \frac{1}{{{\rho _{\text{m}}}{c^2}}}. $

c为混合物的声速. 将式(3)代入式(2)得到

(4)$ {\psi _{\text{m}}} = \frac{1}{{{c^2}}}. $

蒸汽和液体的状态方程分别为

(5)$ {\rho _{\text{v}}} = {\psi _{\text{v}}}p, $

(6)$ {\rho _{\text{l}}} = \rho _{\text{l}}^0+{\psi _{\text{l}}}p. $

式中：ρ_v、ρ_l分别为蒸汽密度和液体密度，ψ_v、ψ_l分别为蒸汽和液体的压缩系数，$\rho _{\text{l}}^0$为未压状态下的液体密度. 混合物的质量分数为

(7)$ \gamma_{{\mathrm{m}}}= \frac{{{\rho _{\text{m}}} - {\rho _{{\text{ls}}}}}}{{{\rho _{{\text{vs}}}} - {\rho _{{\text{ls}}}}}}. $

式中：ρ_ls、ρ_vs分别为液体和蒸汽的饱和密度. 液体饱和密度降低到蒸汽饱和密度，出现空化现象. 当$\gamma_{\mathrm{m}} $=0时，为水，当$\gamma_{\mathrm{m}} $=1时，发生空化，即流场局部全部转化为蒸汽. 混合物密度的计算式为

(8)$ {\rho _{\text{m}}} = (1 - \gamma_{\mathrm{m}} )\rho _{\text{l}}^0+\left( {\gamma_{\mathrm{m}} {\psi _{\text{v}}}+(1 - \gamma_{\mathrm{m}} ){\psi _{\text{l}}}} \right){p_{{\text{sa}}}}+{\psi _{\text{m}}}(p - {p_{{\text{sa}}}}). $

式中：p_sa为饱和压力，$\psi _{\mathrm{m}}$为混合物的压缩系数. 采用Linear模型计算混合物的压缩系数^[16]：

(9)$ {\psi _{\text{m}}} = \gamma_{\mathrm{m}} {\psi _{\text{v}}}+(1 - \gamma_{\mathrm{m}} ){\psi _{\text{l}}}. $

定义μ_v、μ_l分别为蒸汽和液体的黏度，则混合物的黏度计算式为

(10)$ {\mu _{\text{m}}} = \gamma_{\mathrm{m}} {\mu _{\text{v}}}+(1 - \gamma_{\mathrm{m}} ){\mu _{\text{l}}}. $

定义U为流体速度矢量，连续方程^[17]为

(11)$ \frac{\partial {\rho }_{\text{m}}}{\partial t}+\nabla ·({\rho }_{\text{m}}{\boldsymbol{U}})=0. $

动量方程为

(12)$ \frac{\partial {\rho }_{\text{m}}{\boldsymbol{U}}}{\partial t}+\nabla · ({\rho }_{\text{m}}{\boldsymbol{UU}})=-\nabla p+\nabla ·\left({\mu }_{\text{eff}}(\nabla {\boldsymbol{U}}+{(\nabla {\boldsymbol{U}})}^{\text{T}})\right). $

根据Mettin等^[18-19]的研究，将作用在空泡上的力表示为

(13)$ {F_{\text{i}}} = {F_{{\text{B}}}}+{F_{{\text{D}}}}. $

式中：F_B、F_D分别为次Bjerknes力和黏性阻力，根据双空泡理论推导凹坑壁面下单空泡溃灭时的次Bjerknes力及阻力，得到

(14)$ {F_{{\text{B}}}} = \frac{1}{3}{\text{π}} R_{\text{i}}^{\text{3}}\rho (2{R_{\text{j}}}{{\boldsymbol{R}}'_{\text{j}}}^2+{{R}}_{\text{j}}^{\text{2}}{{\boldsymbol{R}}''_{\text{j}}}^{\text{2}})\frac{1}{{{\xi ^{\text{2}}}}} {\boldsymbol{e}} , $

(15)$ {F_{{\text{D}}}} = - 12{\text{π}} \mu {R_{\text{i}}}{v_{\text{i}}}. $

式中：R_i为主空泡半径，ρ为流场密度，R_j为镜像空泡半径，${\boldsymbol{R}}_{\mathrm{j}}' $为壁面坍缩速度，${\boldsymbol{R}}_{\mathrm{j}}'' $为壁面坍缩加速度，μ为运动黏度系数，v_i为空泡运动速度. 定义无量纲参数ξ为凹坑相对曲率，γ为空泡近壁系数，即空泡中心到凹坑壁面的无量纲距离：

(16)$ \xi = {{{R_{\text{m}}}}}/{{{R_{\text{s}}}}},\;\gamma = {L}/{{{R_{\text{m}}}}}. $

式中：R_m为空泡初始半径，R_s为凹坑半径，L为空泡中心与凹坑壁面的距离.

由于空泡溃灭时会向周围的流场辐射声波，形成次级声场，当次级声场反作用于空泡时，会产生相互作用力（次Bjerknes力），当脉动中的空泡同向振荡时相互吸引，反向振荡时相互排斥，这种现象被称为次Bjerknes效应. 如图1所示，分析不同凹坑相对曲率下空泡溃灭时受到的次Bjerknes力，当凹坑相对曲率增大时，作用在空泡上的次Bjerknes力逐渐增大. 随着空泡溃灭时间增加，次Bjerknes力先表现为引力并逐渐减小，接着表现为斥力，在空泡溃灭末期次Bjerknes力表现为引力，凹坑相对曲率越大，引力越大. 如图2所示为不同凹坑相对曲率下次Bjerknes力随空泡溃灭过程的变化. 对比图1与图2发现，在空泡溃灭的过程中，次Bjerknes力对空泡溃灭产生的作用力比阻力大一个数量级，说明次Bjerknes效应对空泡溃灭时动力学特性影响较大.

如图3所示为研究空泡在凹坑壁面附近的溃灭过程. 所研究流场的环境压力为1.0×10⁵ Pa，流场介质选用清水，密度为998 kg/m³，空泡初始半径为1 mm，空泡内部初始压力为298 K时的饱和蒸汽压为3 160 Pa，蒸汽饱和密度为0.59 kg/m³，考虑重力和表面张力系数的影响，表面张力系数为0.072 N/m. 设计时，将空泡放置于凹坑中心正上方L处，L_d为凹坑的深度.

如表1所示，设计6种凹坑相对曲率方案探究凹坑相对曲率对单空泡性能特性的影响. 在L_d=0.8 mm的前提下，对每种方案的壁面距离进行无量纲化，得到相应的空泡近壁系数γ. 如图4所示，简化不同凹坑相对曲率的凹坑壁面，凹坑深度固定为0.8 mm，研究空泡溃灭过程.

设定R_m=1 mm，设置二维平面计算域面积为50 mm×50 mm. 将空泡放置于计算域底部凹坑壁面附近，避免计算域其他边界对空泡溃灭过程的影响. 如图5所示，对空泡溃灭区域网格进行加密. 图5（b）中圆形部分体积分数为1，表示为蒸汽，其他部分体积分数为0，表示为水. 使用有限体积法离散控制方程，方程中数值离散格式如下：时间导数项使用欧拉隐式格式，梯度项为高斯积分，插值格式采用中心差分格式，面法向梯度项采用未修正的离散格式，网格体心到网格面心的插值使用线性化插值格式. 压力速度的耦合采用压力-速度耦合算法进行循环求解，为了保证数值稳定性，将最大库朗数q和最大声学库朗数设置为0.5；流场的4个壁面速度边界设置为无滑移壁面，压力边界设置为零梯度边界.

为了检查数值模拟结果对网格的依赖性，设计4套网格划分方案，以自由场空泡溃灭过程作为验证实例. 不同网格方案下空泡半径随时间的变化如图6所示，可以看出，不同网格方案下空泡体积的变化过程差异极小，基本上保持一致，说明4套网格方案均具有足够的精度来计算空泡的溃灭过程. 为了保证足够的分辨率来捕捉空泡溃灭过程中释放的压力波，也为了节约计算资源，采用650×650的网格方案进行后续的数值模拟，最小网格尺寸$ {{x}}_{\text{min}} $=0.019 2 mm，沿空泡半径方向网格数量为52.

Cui等^[20]采用相似理论，在静止水箱中进行电火花实验产生空泡并模拟空泡的溃灭过程. 当凹坑相对曲率ξ=0.5，空泡近壁系数γ=1.5时进行实验，球形空泡由于内部压力远小于外界压力，随着空泡溃灭时间增加，空泡体积不断缩小，但横向收缩速度大于纵向收缩速度，空泡最终表现为纵向尺寸大于横向尺寸的椭球形. 空泡上表面收缩速度大于下表面收缩速度，因此在空泡溃灭末期，空泡上表面曲率大于下表面曲率，导致空泡向内凹陷，最终分裂为2个环形空泡. 如图7所示，由于空泡溃灭时间较短，数值模拟与实验条件完全一致，且实验结果与数值模拟结果基本一致，证明本研究采用的数值模拟方法可靠. 本研究通过数值模拟得到不同凹坑曲率下空泡坍塌到最小半径所需要的时间t^*，进而计算空泡溃灭时间无量纲系数τ：

(17)$ \tau = {t}/{{{t^*}}}. $

式中：t为空泡溃灭过程时间.

推导双空泡理论，得到空泡间次Bjerknes力计算式，进行空泡溃灭时间仿真计算. 如图8所示，随着凹坑相对曲率增大，空泡溃灭时间逐渐增加，在相同凹坑相对曲率下，空泡近壁系数越大，空泡溃灭时间越短，同一空泡近壁系数下，凹坑相对曲率越大，空泡所受到的次Bjerknes效应越大，空泡溃灭时间随之延长，说明凹坑相对曲率越大，对空泡溃灭过程的阻碍作用越大. 空泡近壁系数越大，空泡溃灭时间越小，说明随着近壁系数增大，次Bjerknes效应对空泡的作用逐渐减小，空泡对壁面的作用随之减小. 基于图8建立空泡溃灭时间系数随凹坑相对曲率和空泡近壁系数变化的关联关系：

(18)$ \tau = - 14.58\exp \left({\xi }/{{0.22}}\right)+(14.83 - 0.06\gamma ). $

由式（18）可知，空泡溃灭时间系数与凹坑相对曲率近似呈指数函数关系并单调递增，与空泡近壁系数近似呈一次函数关系.

如图9所示，在同一个空泡近壁系数下，随着凹坑相对曲率逐渐增大，空泡轮毂逐渐向两侧延伸，空泡溃灭形态为环形空泡在二维平面的2个环形轮廓，空泡尾部相互接近且内表面呈一定角度. 如图10所示，进一步研究空泡射流速度与凹坑相对曲率以及空泡近壁系数之间的变化关系，随着凹坑相对曲率增大，空泡溃灭时射流速度逐渐减小，当ξ=0时，空泡溃灭时射流速度最大；随着空泡近壁系数逐渐增大，空泡溃灭时射流速度逐渐增大. 空泡溃灭射流速度变化说明随着凹坑相对曲率和壁面与空泡间的次Bjerknes效应逐渐增大，空泡溃灭对壁面的作用增强，空泡溃灭过程加快. 通过非线性拟合，使用20 ℃清水中声速V*=1 500 m/s将所得射流速度处理为马赫数：

(19)$ Ma = {V}/{{{V^*}}}. $

式中：V为空泡射流速度. 建立空泡射流马赫数随凹坑相对曲率和空泡近壁系数变化的关联关系：

(20)$ Ma = 0.015\exp \left({{ - \xi }}/{{0.36}}\right)+(0.007+0.009\gamma ). $

由式（20）可知，空泡射流马赫数与凹坑相对曲率呈指数关系并单调递减，与空泡近壁系数近似呈一次函数关系.

通过非线性拟合，使用环境压强p=1.0×10⁵ Pa，计算得到压力冲击载荷峰值系数：

(21)$ \eta = {{{p_{\max }}}}/{p}. $

式中：p_max为壁面压力冲击载荷峰值点. 如图11所示，随着空泡溃灭时间的增加，壁面中心的压力冲击载荷均呈先增大到最高点后快速下降直至平缓的趋势. 随着空泡近壁系数的增加，壁面中心的压力冲击载荷逐渐减小，并且在同一空泡近壁系数下，壁面中心的载荷压力峰值最大的工况表现出壁面曲率由0.2向0.6增大的趋势. 随着凹坑相对曲率的增加，壁面中心压力冲击载荷到达峰值的时间逐渐延缓，说明空泡初始位置距离壁面越近，次Bjerknes效应越强烈，对壁面的作用越大；不同的壁面曲率使得空泡受到的作用力不同，改变了空泡溃灭过程中的几何形状，进而影响空泡溃灭对壁面的压力冲击载荷. 在空泡近壁系数由0.5向1.0增加的过程中，空泡与壁面间的距离逐渐增大，而空泡溃灭末期凹坑壁面中心的压力冲击载荷峰值逐渐减小，说明在空泡与壁面距离增大的过程中，空泡溃灭产生的压力波传播时间延长，导致凹坑壁面中心受到的压力减小，空泡对凹坑壁面中心的作用随之减小，使得凹坑壁面中心压力冲击载荷峰值随着空泡近壁系数的增加而逐渐向凹坑相对曲率更大的方向移动.

如图12所示为凹坑壁面中心冲击载荷的时间历程. 由图12得到凹坑壁面附近压力波的传播过程图，如图13所示. 以ξ=0.5、γ=1.2的工况为例，选取a、b、c、d、e、f时刻研究空泡溃灭末期压力波传播特性. 结合对二维平面空泡溃灭过程的研究发现，当空泡由a向b发展时，空泡上表面压力逐渐增大，空泡表面不断向内凹陷，在b时刻射流贯穿空泡，在二维平面中为中空的环形空泡的2个环形轮廓并向斜上方传播2束压力波. 壁面中心所受到的压力急剧上升，在d时出现第一个峰值. 当d向e发展时，在环形空泡尾部出现高压区，壁面中心受到的压力出现短暂下降，说明2个环形空泡间存在相互作用力并对流场产生扰动. 当空泡溃灭发展到f时刻时，空泡完全溃灭，同时产生2束环形压力波向周围传播，压力峰值的同向叠加在底部靠近壁面的位置产生高压区，壁面中心受到的压力载荷出现远大于其他时刻的第二个压力峰值，说明2束压力波相互叠加，对凹坑壁面中心载荷压力的影响剧烈，对凹坑壁面的作用力增强.

在环境压力为1个大气压，20 ℃ 清水的不同凹坑相对曲率和空泡近壁系数下，用CavitatingFoam对初始直径为1 mm的单空泡进行数值模拟，得到以下结论. 1）在同一空泡近壁系数下，凹坑相对曲率与空泡溃灭时间呈指数关系，在同一凹坑相对曲率下，空泡近壁系数与空泡溃灭时间呈一次函数关系. 2）在凹坑相对曲率逐渐增大的过程中，空泡溃灭产生的射流速度逐渐增大，空泡射流马赫数与凹坑相对曲率呈指数关系，与空泡近壁系数呈一次函数关系. 在相同的空泡近壁系数下，空泡壁面压力中心载荷峰值呈现凹坑相对曲率由0.2向0.6变化的趋势. 3）对空泡溃灭时向周围流场辐射的两束压力波进行监测，发现随着空泡溃灭时间的增加，压力波首先在空泡溃灭的上表面产生高压区，接着在空泡溃灭的下表面产生高压区，当空泡完全溃灭后，2束压力波峰值同向叠加，出现压力峰值，此时高压区的压力载荷最大. 本研究揭示了空泡溃灭时压力波的传播特性，得到单空泡在凹坑附近溃灭时的动力学特性，总结出的普适性规律能够为抑制空化破坏及利用空化的工程应用提供理论依据和技术参考. 由于多空泡和不规则壁面的动力学特性研究还存在不足，未来计划探索多空泡与不同形状壁面的动力学特性.