摘要:基于摄动理论和广义伯努利方程, 推导出单气泡在超声场中径向振动方程、平移方程和气泡形变方程. 数值计算这3个方程, 可以得到气泡半径、气泡中心的位移和气泡形变随时间的演化图. 计算结果表明: 当气泡初始半径和驱动声压不变时, 气泡中心初始平移速度增大, 气泡径向振动几乎不变, 但气泡中心位移和形变量增大, 气泡非球形振动愈加明显. 当初始平移速度比较小时, 气泡的${R_{0}}\text {-}p_{\rm a}$相图中, 不稳定区域仅集中在高驱动声压区域. 随着气泡中心初始平移速度不断增大, 半径和驱动声压均较小的区域开始呈现不稳定性, 且整体不稳定空间范围逐渐增大. 另外, 气泡在声驻波场中不同位置呈现出不同的振动特征. 离波腹点越近的气泡, 其径向振动幅度越大, 但气泡的平移和形变量变化很小, ${R_{0}}\text {-} p_{\rm a}$相图中不稳定性区域平面分数之间的误差小于4%.
关键词: 空化泡/
平移/
非球形振动/
非稳定性

English Abstract

全文HTML

--> --> -->

当强超声通过液体时, 液体中微小气泡(即空化核)会成长为肉眼可见的、微米量级的气泡, 这种现象称为声空化^[1]. 形成的气泡叫做空化泡. 空化泡在声场的作用下, 其径向运动表现出先缓慢膨胀, 然后急剧塌缩, 最后反弹的周期性脉动现象. 空化泡塌缩至最小半径前后, 其内部会产生极端的高温高压, 气泡内部物质可能发生电离, 导致光辐射, 这种现象称为声致发光^[2]. 空化泡内部的高温高压已被应用到超声清洗^[3]、纳米材料制备^[4]和催化化学反应^[5]等领域.
声致发光是空化泡急剧运动导致的结果, 空化泡动力学一直以来是声空化和声致发光研究的课题之一. 早在1917年, Rayleigh^[6]就建立了著名的球形空化泡动力学模型. 1954年, Plesset^[7]从流体力学基本方程组出发, 研究了单个球形空化泡的形状不稳定性. 1980年, Keller和Miksis^[8]建立了压缩液体中空化泡的动力学模型. 1996年, Hilgenfeldt等^[9]从扩散平衡不稳定性方面系统地研究了单个空化泡的非稳定性. 2005年, An等^[10]考虑了非球形声压对单泡声致发光气泡不稳定性的影响, 计算了气泡的不稳定性相图. 2013年, Zhang和An^[11]研究了空化泡链中气泡的不稳定性, 具体分析了气泡的形状不稳定性、扩散平衡不稳定性以及位移不稳定性. 2012年, Liang等^[12]研究了两个气泡在球形声场非球形振动; 2017年, 他们研究了两个气泡在非球形声场中的动力学行为^[13]. 在这些研究中, 气泡均被假定在一个固定的位置.
但实验观测结果表明, 气泡在超声驱动下的运动是非常复杂的. 1998年, Madrazo等^[14]观测到单气泡在超声场中的非球形脉动. 1999年, Barbat等^[15]观察到两个气泡在声场中周期性平移运动. 2007年, Flannigan和Suslick^[16]观测到硫酸溶液中单个快速平移气泡的声致发光光谱. 2012年, Cui等^[17]观察到含不同酒精量水溶液中单个空化泡平移、抖动行为. 2019年, Wu等^[18]研究了强超声场中多气泡环境下单个气泡的振动和平移. 2021年, Wu等^[19]用高速相机观测到空化泡在油-水界面处的动力学行为. 这些工作说明气泡在超声场中的动力学行为包含了气泡的径向脉动、气泡的平移和气泡表面的形变. 为了研究气泡脉动、平移和形变的动力学机制, 一些****做了很多有价值的工作. 2001年, Doinikov^[20]运用拉格朗日理论体系研究了两个空化泡在超声场中径向脉动和平移运动; 2002年, Doinikov^[21]研究了平面驻波声场中单个空化泡平移运动和径向振动. 2015年, 沈壮志^[22]研究了声驻波场中空化泡在不同位置处的动力学特性. 2018年, 马艳等^[23]基于拉格朗日方程得到两个空化气泡的动力学方程, 分析了两个气泡的耦合振动及其形状的不稳定性. 2020年, Zhang等^[24]详细研究了两个气泡在声场中的平移动力学机制. 这些工作只考虑了气泡的径向脉动和平移或气泡的径向脉动和形变两种情况, 没有把气泡的径向脉动、平移和形变进行耦合. 1995年, Feng和Leal^[25]首次建立了单个气泡在超声场的径向脉动、气泡平移以及气泡形变的动力学耦合模型. 2002年, Reddy和Szeri^[26]研究了单个气泡的径向脉动、气泡平移随时间的演化关系. 2004年, Doinikov等^[27]研究了单个气泡在声驻波场中径向脉动、气泡平移和形变三者的非线性耦合运动. 2009年, Mettin和Doinikov^[28]研究了单个气泡在高频声场中不同位置的稳定性特征. 但气泡径向脉动、平移和形变之间的相互作用以及声场分布对气泡稳定性的影响还不十分清楚. 本文建立了单个气泡径向脉动、平移和形变相互耦合的动力学模型, 数值计算这3个方程，得到了气泡初始平移速度和驻波声场分布分别对气泡的脉动、平移、形变、非球形振动和稳定性的影响.

考虑不可压缩液体中的1个气泡, 且液体流动是无旋的. 以气泡中心为球坐标原点, 如图1所示. 气泡在超声波驱动和液体黏滞力作用下, 可能会沿x方向平动. 根据势流理论, 液体中气泡附近的速度势(?)满足拉普拉斯方程$ \nabla^{2}\phi $ = 0, 则液体中单个空化泡附近的速度势可以假定为^[12,20]
图 1 含有气泡脉动、平移和形变的几何图形
Figure1. Geometry for single bubble with pulsation, translation and shape perturbation.

$ \phi = \sum\limits_{n = 0}^{\infty}\textstyle\dfrac{A_{n}}{r^{n+1}}{\rm P}_{n}({\rm {cos}}\theta), $

气泡表面的方程可以表述为

$ F(r,\theta,\phi,t) = r-S(\theta,\phi, t) = 0, $

如果气泡的形状在Z轴方向是旋转对称的, 则气泡表面函数可以近似表示为

$ S(\theta,t)\approx R+\epsilon a_{2}{\rm{P}}_{2}(\mu), $

其中, R表示气泡振动的球形半径; $ a_{2} $表示气泡非球形振动部分的振幅; $\epsilon $是摄动参数, 其值设定小于1, 保证$ a_{2} $小于R. 为了简化, 这里只取第2阶勒让德函数($ {\rm{P}}_{2} $(μ))描述气泡的形变, μ = cosθ.
根据方程(1)—(3), 气泡附近总的速度势近似为

$ \phi\approx\frac{A_{0}}{r}{\rm{P}}_{0}(\mu)+\frac{A_{1}}{r^{2}}{\rm{P}}_{1}(\mu)+\epsilon \frac{A_{2}}{r^{3}}{\rm{P}}_{2}(\mu), $

考虑气泡的径向脉动、平移和形变, 气泡表面的法向速度近似满足方程

$ \frac{\partial F}{\partial t}-\dot{x}{\rm{cos}}\theta+(\nabla\phi)\cdot{{\boldsymbol{e}}_{\boldsymbol{r}}} \frac{\partial F}{\partial r} = 0, $

其中, $ F\equiv F(r, \theta, t) $, $ {{\boldsymbol{e}}_{\boldsymbol{r}}} $是速度法向方向单位矢量,$\nabla$表示沿 r方向的梯度, $ \dot{x} $表示x对时间的一阶导数. 把方程(2)—(4)代入方程(5), 然后对方程(5)进行?泰勒级数展开, 可以得到

$ \epsilon^{0}:\; \left( {\frac{{{A_0}}}{{{R^2}}} + \dot R} \right){{\rm{P}}_0} + \left( {\frac{{2{A_1}}}{{{R^3}}} + \dot x} \right){{\rm{P}}_1} = 0, $

$ \epsilon^{1}:\;\left( { - \dot{a_{2}} - \frac{{3{A_2} + 6{A_1}{a_2}{{\rm{P}}_1}}}{{{R^4}}} + \frac{{2{A_0}{a_2}}}{{{R^3}}}} \right){{\rm{P}}_2} = 0, $

其中$ {\rm{P}}_{0}\equiv {\rm{P}}_{0}(\mu) $, $ {\rm{P}}_{1}\equiv {\rm{P}}_{1}(\mu) $, $ {\rm{P}}_{2}\equiv {\rm{P}}_{2}(\mu) $, $ \dot{R} $表示R对时间的一阶导数. 根据勒让德函数的正交性^[29], 由方程(6)和方程(7)解得:

$ A_{0} = -R^{2}\dot{R}, $

$ A_{1} = -\frac{1}{2}R^{3}\dot{x}, $

$ A_{2} = -\frac{1}{3}R^{3}(R\dot{a_{2}}+2a_{2}\dot{R}+3a_{2}\dot{x}P_{1}). $

在不可压缩的液体中, 如果流体是无旋的, 则液体中1个平动气泡周围的速度势满足广义伯努利方程^[30]:

$ p - {p_\infty } = - \rho \left[ {\frac{{\partial \phi }}{{\partial t}} - \dot x \cdot \nabla \phi + \frac{1}{2}{{(|\nabla \phi |)}^2}} \right], $

其中, ρ是液体密度, p和$ p_{\infty} $分别表示液体中r处和无穷远处的压强. 把方程(4)和方程(8)—(10)代入方程(11), 然后对方程(11)进行?泰勒级数展开, 利用勒让德多项式的性质^[29], 可以得到

$ \epsilon^{0}:\; R\ddot{R}+\frac{3}{2}\dot{R}^{2}-\frac{\dot{x}}{4} = \frac{p-p_{\infty}}{\rho}, $

$ \frac{3}{2}\dot{R}\dot{x}+\frac{1}{2}R\ddot{x} = 0, $

$ \epsilon^{1}:\; \dot{a}_{2}\dot{R}-\frac{a_{2}\dot{x}^{2}}{14R}+\frac{1}{3}R\ddot{a}-\frac{1}{3}a_{2}\ddot{R} = 0, $

其中, $ \ddot{R} $, $ \ddot{x} $分别表示R, x对时间的二阶导数. $ p\equiv $$ p $(r)的表达式为^[12]

$ p(r) = p_{\rm g}-p_{\rm d}-\frac{2\sigma}{R}-\frac{4\eta\dot{R}}{R}+\frac{R}{c}\frac{\rm d}{{\rm d}t}(p_{\rm g}-p_{\rm d}), $

$ {p_{\rm g}} = \left( {{p_0} + \frac{{2\sigma }}{{{R_0}}}} \right){\left( {\frac{{R_0^3 - {h^3}}}{{{R^3} - {h^3}}}} \right)^\gamma }, $

$ p_{\rm d} = -p_{\rm a}{\rm {sin}}(wt), ~~w = 2\pi f, $

其中, σ是气泡的表面张力, $ p_{0} $是气泡内初始压强, η表示液体的黏度系数, $ R_{0} $表示气泡初始半径, h表示范德瓦耳斯排斥半径, h = 8.5/$ R_{0} $, γ表示气体的绝热指数, $p_{\rm a}$和f分别表示驱动声压幅度和超声频率. 从方程(12)—(14)能明显看出, 若不考虑气泡的平移, 方程(12)就退化为著名的单气泡动力学Rayleigh-Plesset (RP)方程^[6,31], 方程(14)也退化为文献[32]中的单气泡形变方程.
值得注意的是, 在推导方程(12)—(14)的过程中, 因在假定速度势时, $ \epsilon^{0} $项只考虑了气泡的径向振动和平移, 所以只选取了$ \epsilon^{0} $中的$ {\rm{P}}_{0} $和$ {\rm{P}}_{1} $项, 分别用来描述气泡的径向振动和平移, 省略$ {\rm{P}}_{2} $项. 同理, 对于含有$ \epsilon^{1} $的式子, 类比速度势中的微扰项, 只取$ \epsilon^{1} $中的$ {\rm{P}}_{2} $项来描述气泡的形变, 省略$ {\rm{P}}_{0} $, $ {\rm{P}}_{1} $, $ {\rm{P}}_{3} $和$ {\rm{P}}_{4} $项.

为了分析单个气泡在液体超声场中的脉动、平移和形变等动力学特性, 采用龙格-库塔法数值计算方程(12)—(14). 初始条件为$ t = 0 $时, $ R = R_{0} $,$ {\rm d}R/{\rm d}t = 0 $, $ x = 0 $, ${\rm d}x/{\rm d}t = v_{x0}$, $ a_{2} = 1\times10^{-10} $ m, $ {\rm d}a_{2}/{\rm d}t = 0 $. 本文使用的物理量参数见表1.

物理量	单位	数值
液体密度ρ	kg/m3	1000
液体中声速c	m/s	1481
液体表面张力系数σ	N/m	0.072
液体粘度系数η	Pa·s	0.001
超声频率f	Hz	2.5 × 104
气泡内初始压强$ p_{0}$	Pa	1.013 × 105
环境压强$ p_{\infty}$	Pa	1.013 × 105
绝热指数γ		1.4

表1数值计算中使用的物理量参数
Table1.Physical parameters used in the numerical calculation.

2

-->

3.1.初始平移速度对气泡的脉动、平移位移和形变的影响

为了研究初始平移速度对气泡的脉动、平移和形变的影响, 联合方程(12)—(17), 分别计算模拟气泡在不同初始平移速度($v_{x0}$ = 0, 0.1, 1, 3, 5 m/s)下, 其半径、中心平移位移和形变随时间的演化图像, 结果如图2所示. 可以看出, 随着初始平移速度的增加, 即从0增加到5 m/s时, 气泡半径的大小及其周期性演化几乎没有变化, 但位移和形变的大小不断增大, 形变的增大最终导致气泡的破裂.
图 2 不同初始平移速度$v_{x0}$下, (a)气泡径向半径、(b)气泡中心平移以及(c)气泡表面形变随时间的演化. $R_{0}$ = 4.5 μm, $p_{\rm a}$ = 1.15 × 10⁵ Pa
Figure2. Evolutions of (a) radius, (b) translation and (c) deformation with time for a cavitation bubble with different initial translational velocity ($v_{x0}$), respectively. $R_{0}$ = 4.5 μm, $p_{\rm a}$ = 1.15 × 10⁵ Pa

2 -->

3.2.初始平移速度对气泡非球形振动的影响

为了模拟气泡的非球形振动, 把气泡的形状近似看成1个轴对称的椭球. 椭球两个轴的伸长范围由$ L = {a_{2}}/{R} $决定, 第3个轴看成是恒定对称的. 当$ L = 0 $时, 气泡的形状是完全球形的. 当$ 0 < L < 1 $时, 气泡变成1个椭球. 当$ L\geqslant1 $时, 气泡被认为破裂. 考虑椭球的对称性(即气泡在每个时刻都被看成是对称的椭球), 应用文献[13, 33, 34]中方法, 把气泡投影到XY平面, 可以得到气泡振动时的形状. 图3(a)—(d)分别表示气泡中心初始平移速度$v_{x0}$分别等于0, 1, 3和5 m/s时, 气泡在不同时刻的振动形状. 图中t = 0表示气泡初始时刻的形状, t = 13.42 μs表示气泡膨胀到最大半径时刻的形状, t = 18.61 μs表示气泡坍塌到最小半径时刻的形状. 当$v_{x0}$分别等于0, 1, 3和5 m/s时, 在t等于58.04, 56.11, 31.76和24.76 μs (1个声周期40 μs)时, 气泡形状严重变形, 容易破灭. 这说明气泡初始平移速度越大, 气泡存活的时间越短.
图 3 不同初始平移速度条件下, 不同时刻气泡振动的形状　(a) $v_{x0}$ = 0; (b) $v_{x0}$ = 1 m/s; (c) $v_{x0}$ = 3 m/s; (d) $v_{x0}$ = 5 m/s. 其中$R_{0}$ = 4.5 μm, $p_{\rm a}$ = 1.15 × 10⁵ Pa
Figure3. Simulations of shapes of a gas bubble’s oscillation at different times under different initial translational velocity: (a) $v_{x0}$ = 0; (b) $v_{x0}$ = 1 m/s; (c) $v_{x0}$ = 3 m/s; (d) $v_{x0}$ = 5 m/s. $R_{0}$ = 4.5 μm, $p_{\rm a}$ = 1.15 × 10⁵ Pa

2 -->

3.3.初始平移速度对气泡非稳定性的影响

气泡在超声场中的非稳定性空间范围是表征气泡动力学的重要指标之一. 本文采用Rayleigh-Taylor (RT)判据$ |a_{2}/R|\geqslant1 $来确定气泡的不稳定性空间^[9]. RT判据是指随着气泡的周期性振动, 其形变量$ a_{2} $可能会不断增大. 当$ a_{2} $大于气泡球形半径R时, 气泡就会破灭. 根据该判据, 计算得到4种不同初始平移速度的$R_{0}\text {-} p_{\rm a}$相图, 如图4所示. 图中蓝色部分表示不稳定性区域, 灰色部分表示固定不变区域. $S_{\rm a}$表示不稳定区域的面积百分数, 其值等于蓝色区域面积与灰色区域面积(恒定)的百分比.
图 4 不同初始平移速度下, 气泡的${R_{0}}\text {-}$$p_{\rm a}$相图　(a) $v_{x0}$ = 0; (b) $v_{x0}$ = 1 m/s; (c) $v_{x0}$ = 3 m/s; (d) $v_{x0}$ = 5 m/s. $p_{\rm a}= $$ 1.15$× 10⁵ Pa, $R_{0}$ = 4.5 μm
Figure4. ${R_{0}}\text {-}$$p_{\rm a}$ phase diagram of a gas bubble under different initial translational velocity: (a) $v_{x0}$ = 0; (b) $v_{x0}$ = 1 m/s; (c) $v_{x0}$ = 3 m/s; (d) $v_{x0}$ = 5 m/s. $p_{\rm a}=1.15$ × 10⁵ Pa, $R_{0}$ = 4.5 μm

从图4(a)—(d)可以看出, 当$ v_{x0} $较小时(图4(a)和图4(b)), 不稳定区域仅集中在高驱动声压区域 ($p_{\rm a}\geqslant$ 0.9 × 10⁵ Pa); 随着$ v_{x0} $的不断增大(图4(c)和图4(d)), 初始半径和驱动声压均在较小的区域开始呈现不稳定性, 整体的不稳定性空间范围逐渐增大.
2 -->

3.4.声驻波场中不同位置处气泡的非球形振动和非稳定性

为研究声场分布对单气泡的径向振动、平移、形变、非球形振动和稳定性的影响, 计算模拟气泡在声驻波场中不同位置处的动力学行为. 声驻波场可以表示为$ p_{\rm d} = -p_{\rm a} \sin (wt) \cos (kd) $. 这里的k表示声波的波数, $ k = 1/\lambda = f/c $, d表示气泡到波腹点的距离, λ表示声波波长.
图5给出波腹点不同位置处, 气泡的半径、平移位移和形变随时间的演化图像. 从图5可以看出, 当d从0变化到λ/2时, 气泡振动过程中径向部分的最大半径压缩比$ R/R_{0} $从2.5减小到1.5, 如图5(a)所示, 这说明气泡在波腹处的振动比较剧烈. 气泡中心的位移和气泡形变变化较小 (如图5(b)和图5(c)).
图 5 不同位置处, 气泡半径、平移位移和形变随时间的演化图像($p_{\rm a}=1.15$ × 10⁵ Pa, $R_{0}$ = 4.5 μm, $v_{x0}$ = 5 m/s)
Figure5. Evolution of bubble’s radius, translation and deformation at different distance d from the antinodal point of acoustical wave. $p_{\rm a}=1.15$ × 10⁵ Pa, $R_{0}$ = 4.5 μm, $v_{x0}$ = 5 m/s.

图6是对应d在0, λ/5, λ/4和λ/2处, 气泡的非球形振动图像. t = 0是初始时刻的图像. t = 10.63 μs是气泡膨胀到最大半径时的图像. t = 23.60 μs是气泡坍缩到最小半径时的图像. t = 24.02, 24.60, 34.60 μs是气泡反弹阶段的图像. 从图6可以看出, 气泡在d = 0 (波腹处), 气泡的形变较小. 随着d 值增大到λ/2处, 气泡的最大半径明显变小. 当t = 35.60 μs时, 气泡扭曲最大, 说明离波腹点的距离越大, 气泡越容易破裂.
图 6 不同位置处, 不同时刻气泡振动的形状($p_{\rm a}=1.15$ × 10⁵ Pa, $R_{0}$ = 4.5 μm, $v_{x0}$ = 5 m/s)　(a) d = 0; (b) d = λ/5; (c) d = λ/4; (d) d = λ/2
Figure6. Simulations of shapes of a gas bubble’s oscillation at different times at different distances from the antinodal of acoustical wave: (a) d = 0; (b) d = λ/5; (c) d = λ/4; (d) d = λ/2. $p_{\rm a}=1.15$ × 10⁵ Pa, $R_{0}$ = 4.5 μm, $v_{x0}$ = 5 m/s.

图7是根据RT判据$ |a_{2}/R|\geqslant1 $, 数值计算得到的气泡在声驻波场中不同位置处的${R_{0}}{\text {-}}p_{\rm a}$相图. 从图7(a)—(d)可以看出, 4个不同位置处的不稳定性空间区域的面积分数非常接近, 相互之间的差值小于4%. 这说明声驻波场分布对气泡的整体不稳定性空间范围影响较小.
图 7 不同位置处, 气泡的$R_{0}$-$p_{\rm a}$相图($p_{\rm a}=1.15$ × 10⁵ Pa, $R_{0}$ = 4.5 μm,$v_{x0}$ = 5 m/s)　(a) d = 0; (b) d = λ/5; (c) d = λ/4; (d) d = λ/2
Figure7. $R_{0}$-$p_{\rm a}$ phase diagrams of a gas bubble at different distances from the antinodal of acoustical wave: (a) d = 0; (b) d = λ/5; (c) d = λ/4; (d) d = λ/2. $p_{\rm a}=1.15$ × 10⁵ Pa, $R_{0}$ = 4.5 μm, $v_{x0}$ = 5 m/s.

基于摄动理论和广义伯努利方程, 建立了描述单气泡在超声场中的径向振动、平移和形变的数学模型. 数值模拟结果表明, 当气泡的初始半径和驱动声压幅度不变时, 随着气泡中心初始速度增大, 气泡的非球形振动更加明显, R₀-p_a相图中非稳定性空间范围增大. 分析了离波腹点不同位置处气泡的径向振动、平移、形变、非球形振动和稳定性. 在1个声波长内, 声场位置对气泡的平移和不稳定性空间范围影响较小, 但对气泡的径向振动影响较大. 气泡在波腹点的振动幅度最大, 距离波腹点越远, 振动幅度越小. 本文的研究结果对研究空化泡动力学有一定的指导意义.