单泡超声空化仿真模型的建立及其动力学过程模拟*

徐 珂 许 龙

(中国计量大学理学院 杭州 310018)

0 引言

超声空化通常是指液体中存在的微气泡在超声波作用下振荡、生长、收缩和崩溃等一系列非线性动力学过程[1]。1895年，Thorneycroft 等[2]在观察到经过短时间航行的潜艇螺旋桨推进器的空蚀现象之后，第一次提出了空化的概念。1917年，Rayleigh[3]给出了Rayleigh 气泡动力学模型，Plesset[4]、Noltingk 等[5]、Neppiras 等[6]和Poritsky[7]考虑了液体的表面张力、黏滞性和驱动声场等因素，给出Rayleigh-Plesset 方程。后面关于泡及泡群的理论研究[8-14]主要基于球形空化泡在考虑更复杂情况下的Rayleigh-Plesset方程的修正。

以上理论研究都是在基于空化泡球形变化基础上做的分析，实际情况下的空化泡运动过程并不会一直保持球形。为了观察实际情况下的动力学过程，许多学者进行了大量实验研究[15-22]。但实验研究受实验设备的限制很难准确观察流体中各位置相关物理参数。随着计算流体动力学(Computational fluid dynamics，CFD)技术的发展，观察空化泡运动过程中空化泡形态以及流体各位置基本物理量成为可能。目前基于CFD 的空化泡的仿真模拟研究，多数工作[23-29]都是对近壁面以及静态流场中半径比较大的空化泡的溃灭形变研究，很少有人研究声场作用下微小空化泡的形态变化。

本文基于有限元仿真方法，构建了流体环境中单泡的有限元仿真模型，采用流体动力学控制方程和流体体积分数(Volume of fluid，VOF)模型，模拟了超声驱动下水中单泡的演化过程并与实验拟合结果进行对比，进一步探讨了单泡形态变化过程中泡内及泡外单泡附近压强和气体密度分布。用有限元软件对单泡的动力学行为进行仿真，仿真分析结果与实验结果符合良好。本文建立的单泡超声空化有限元仿真模型为后续双泡以及泡群超声空化动力学过程模拟提供有益参考。

1 仿真模型

1.1 控制方程

在计算空化泡演变过程中引入基本假设条件：(1)液体为不可压缩牛顿流体，且为湍流流动；(2)不考虑重力作用影响；(3)忽略体积力的作用。

质量连续方程为[30]

式(1)中，ρ为流体密度，μ为流体速度。

Navier-Stokes(N-S)方程为[30]

式(2)中，ε为流体运动黏度，P表示流体压力。

空化气泡状态变化整个过程涉及到气-液两相流。为追踪水中气泡的边界变化，引入体积分率函数αq，采用VOF[31]进行气-液交界面的模拟。VOF模型中流体没有相互穿插，体积分数的连续方程(3)及约束条件(4)如下：

其中，q为流体的相，体积分率函数αq表示流体在网格中所占空间的比例。

1.2 有限元模型

为了仿真超声驱动下水中单气泡的演化过程，建立了如图1所示的二维轴对称模型，图中流体区域为长0.6 mm、宽0.3 mm的矩形区域，其中半圆形区域为初始半径为R的空泡。

图1 计算模型Fig.1 Model for the calculation

计算区域中的5 个边界条件分别为：(1)左边axis为轴对称边界条件；(2)右边界和下边界wall为刚性壁面；(3)空泡界面interior为内部边界；(4)上边pressure inlet为声压入口边界。为了模拟压力入口超声变化，使用UDF方法编写输入驱动声压函数P=Pasin(2πft)，其中Pa为驱动声压幅值，f为超声频率。

采用有限元软件求解器中基于压力的求解器进行计算，过程为瞬态，相关参数设置如表1所示。

表1 有限元法相关参数设置Table 1 Parameter setting of finite element method

1.3 仿真模型验证

为了验证本文建立的有限元仿真模型的可靠性，依据文献[19]中超声空化的相关实验数据，对初始半径R0= 6.18 μm 的超声空化泡的动力学行为进行了仿真模拟。根据文献[19]，计算过程中驱动声压幅值Pa= 1.29 atm (1 atm = 101325 Pa)，频率f= 25 kHz，根据超声频率可知驱动信号的周期为40 μs，设置时间步长为0.02 μs，步数为2000 步。仿真计算的空化泡的形态变化过程如图2所示。

图2(a)～(j)为气泡膨胀阶段，约占驱动信号的二分之一周期，其中图2(j)为单泡最大半径时图像；图2(j)～(n)为气泡塌缩阶段，约占驱动信号的十分之一周期，21 μs以后为单泡膨胀压缩做往复运动阶段。由此可以看出气泡形态随时间的演化规律是先缓慢膨胀，到达最大值后，迅速塌缩后消失。由仿真计算的空化泡的形态变化图可知，球形气泡在超声波驱动下其空化泡的变化形状并非理想的球形，而是沿声压激励方向向两边分裂。这也与相关实验观察[20-21]所得结论一致。

图2 不同时刻气泡的形态变化情况Fig.2 The shape change of the bubble at different time

基于上述超声驱动下气泡演变过程的有限元仿真模型，通过计算获得不同时刻空化泡的形态图像后，用图像处理算法，通过读取不同时刻的气泡图像，将气泡区域的像素进行计数并转换成气泡面积，根据气泡面积计算出不同时刻的气泡半径，从而获得如图3蓝色点画线所示的气泡归一化半径随时间的演化曲线。为了验证本文仿真模型的可靠性，将本文仿真计算结果与文献[19]中用实验测试数据拟合的R-P方程的计算结果(图3中的红色点画线)进行对比。

图3 有限元仿真与实验拟合的R-P 方程计算曲线Fig.3 The calculation curve of R-P equation fitted with the experiment and Finite element simulation

由图3可知，仿真计算所得空化泡的膨胀与溃灭曲线与经实验测试数据拟合的R-P 方程计算曲线基本一致。不同之处在于：(1)仿真结果中空化泡无法观测到反弹阶段，这是因为单泡初始半径过小，在微泡缩小时无法适应网格大小；(2)实验拟合的R-P 方程计算曲线略高于仿真曲线，其原因可能是仿真计算过程中考虑了湍流，而用R-P 方程计算时未考虑湍流影响，从而带来二者的差别。

1.4 监测点设置

单泡空化过程中泡内和泡外附近流场中各物理量会发生剧烈变化。在仿真计算过程中，为了监测超声激发单泡空化过程中泡内及泡外流体域内压强和密度随时间的变化，在泡内和泡外设置了9 个监测点，如图4所示。点p0位于单泡中心，点p1、p2、p3、p4在声压激励方向，分别距离单泡中心10 μm、20 μm、30 μm 和40 μm，点p5、p6、p7、p8在声压垂直方向，分别距离单泡中心10 μm、20 μm、30 μm和40 μm。

图4 监测点位置示意图Fig.4 Schematic diagram of monitoring point location

2 结果与讨论

2.1 超声作用下单泡泡内及附近压强变化

对单泡泡内和单泡附近流域压强变化进行了研究。仿真过程中单泡初始半径R0= 6.18 μm，驱动声压幅值Pa= 1.29 atm、频率f= 25 kHz。图5为本文图4中9个监测点处压强随时间变化曲线。

图5 监测点压强变化Fig.5 Pressure change of monitoring point

如图5中p0曲线所示，单泡泡内压强先缓慢减小后迅速增大，在t= 20.60 μs 达到最大值87 atm后迅速减小并做往复变化。结合图2和图3发现，泡内压强变化与单泡体积变化成反比。

当t小于16.30 μs 时，单泡处于膨胀阶段。如图5所示，在t= 6.40 μs 之前，p1与p5曲线、p2与p6曲线、p3与p7曲线和p4与p8曲线分别重合且逐渐递减，表明在单泡膨胀阶段，泡外压强在声波激励方向和声波垂直方向上与泡中心相同距离处压强相同，两个方向的泡外压强随着与泡中心距离的增加逐渐递减。当t= 6.40 μs，单泡半径达到R= 10 μm，此时单泡的表面正好膨胀到p1和p5监测点的位置，随着时间的推移单泡进一步膨胀，p1、p5进入泡内，p1和p5曲线加速向p0曲线靠拢，并在7.2 μs时与p0曲线重合，从此刻开始p1、p5和p0点的压强相同。在t= 10.6 μs 时，单泡半径R=20 μm，此时p2和p6监测点也正好位于单泡表面处，随着时间的增加单泡进一步膨胀，p2和p6监测点也进入泡内，p2、p6曲线也开始加速向p0曲线靠拢，并在12 μs 时与p0曲线重合，从此刻开始p2、p6与p1、p5和p0点的压强都相同。由此可见，在单泡膨胀过程中，当监测点刚进入泡内时，相关检测点的压强与泡中心的压强有一定差别，但随着时间的推移，相关监测点处压强会加速与泡中心点压强趋于一致。

当16.30 μs＜t ＜20.60 μs 时，泡外压强大于泡内压强，单泡处于压缩阶段，随着泡的压缩，泡内压强开始增加。在单泡压缩阶段，泡外声压在声波激励方向和垂直方向压强变化出现明显差异，从而导致在超声空化过程中气泡呈非规则球形变化。

当20.60 μs＜t ＜40 μs 时，如图5所示，泡外压强随距离单泡中心点的距离增大而减小。同时，在声压垂直方向的压强值要高于声压激励方向。对比单泡形态变化，就是在单泡压缩阶段以后，单泡不再是规则的球形，并且沿声压激励方向向两边分裂[20]，导致声压垂直方向的压强值要高于声压激励方向。

2.2 超声作用下单泡泡内及附近气体密度变化

在与上述相同边界条件下对单泡泡内及泡外附近气体密度变化进行了研究。图6为图4中9 个监测位置处气体密度随时间变化曲线。

如图6中p0曲线所示，单泡泡内气体密度先缓慢减小后迅速增大，在t= 20.60 μs 达到最大值103 kg/m3，然后迅速减小并做往复变化。结合图2和图3发现，泡内气体密度变化与单泡体积变化成反比。由图6可知，单泡在缓慢膨胀阶段，泡内密度各处相同且逐渐减小，泡外各处密度随距泡心距离的增大而减小，且泡外密度在声压激励方向和声压垂直方向始终相同。单泡在压缩阶段，单泡附近气体密度增大，并且声压激励方向气体密度变化相对垂直方向较大，在声压垂直方向气体密度要大于声压激励方向气体密度，对比单泡形变，说明单泡压缩过程中发生了沿声压激励方向的形变，单泡在压缩阶段沿声压激励方向向两边分裂，与文献[20]实验观测结果吻合。

图6 监测点的气体密度变化曲线Fig.6 Variation curve of gas density at monitoring points

3 结论

建立了超声空化泡的有限元仿真模型，利用有限元软件模拟了超声驱动下水中单泡的空化动力学过程。基于仿真模拟获得的不同时刻空化泡的演化图像，采用图像处理算法，获得了气泡半径随时间的变化曲线。在此基础上，分析了单泡形态变化过程中泡内及泡外附近压强和气体密度分布。结论如下：

(1)仿真计算所得气泡半径随时间的演化规律是先缓慢膨胀到最大后迅速塌缩，与实验拟合的R-P 方程计算所得气泡半径随时间的变化趋势一致，从而验证了有限元仿真模型的可靠性。

(2)单泡形变过程中，泡内压强与气体密度变化与单泡体积变化成反比。

(3)单泡膨胀阶段，单泡泡内压强与气体密度大于泡外附近压强与气体密度，并且压强与气体密度以单泡中心为球心向外递减。

(4)单泡压缩阶段，单泡附近的压强与气体密度不再是球形均匀变化，而是在声压垂直方向要大于声压激励方向，从而导致在超声空化过程中气泡呈非规则球形变化。

(5)单泡压缩阶段后，在声压垂直方向压强与气体密度要大于声压激励方向气体密度，对比单泡形变，说明单泡压缩阶段过后，沿声压激励方向向两边分裂。

本文研究结果将为模拟更加复杂的超声空化泡及泡群动力学过程提供参考借鉴，相关研究工作将在后续进行。