声场中空化气泡的耦合振动及形状不稳定性的研究∗

马艳林书玉徐洁

1)(陕西师范大学,陕西省超声学重点实验室,西安 710062)

2)(宁夏师范学院物理与电子信息工程学院,纳米结构及功能材料工程技术研究中心,固原 756000)

1 引 言

声致发光测量实验表明:气泡在声场中的形状并非全是球形,在非球形气泡由最大半径急剧塌缩到最小半径的崩溃阶段,这种形变扰动会变得十分显著,甚至使得驱动压强没有达到空化阈值时气泡就发生破裂,影响到气泡的稳定性[1−3].因此,泡群中气泡的稳定振动成为空化研究中的重要课题之一.由于大量气泡之间的相互作用复杂,因此过去关于气泡不稳定性的研究主要集中在单个气泡的形状不稳定性[4−9]上.20世纪50年代,Plessset[10]首次推导了单气泡的形状稳定性方程,研究了单个气泡的不稳定性.Brenner等[11]讨论了引起非球形扰动的两种不稳定机制.Bogoyavlenskiy[12]证明了气泡形变扰动的时间导数随着气泡半径的减小显著增长,是气泡反弹阶段不稳定的主要原因.Wang和Chen[13]引入非球对称的驱动声场,成功地解释了单个气泡的稳定非球形脉动.刘海军和安宇[14]考虑了水蒸气在边界上的凝结和蒸发效应,发现考虑水蒸气后对单气泡稳定区域的确定相对有所改进.钱梦騄等[15]把气泡看作以流体为负载的振子,对单个气泡在膨胀、崩溃和回弹过程中的动力学和形状稳定性进行了讨论,得到了声致发光中气泡动力学的简明物理图像.这些结论的得出表明单气泡形状不稳定性的理论研究已经相对比较成熟[16−18].然而,含气泡液体中,空化的发生大多是以多个气泡和气泡云的形式出现,相邻气泡之间存在着声相互作用[19,20],气泡的运动非常复杂,理论处理非常困难,因此目前关于泡群中气泡稳定性的相关研究较少,仍然处于探索阶段[21,22].本文从双气泡模型出发,研究了两个具有非球形扰动气泡之间的相互作用力对气泡形状不稳定性的影响,并提出了具有非球形扰动的两个气泡之间可能存在一种新的耦合模式,一定条件下,这种耦合模式能够对气泡的稳定振动产生一定的影响,希望对多气泡环境中气泡的稳定振动及多气泡空化提供理论基础.

2 含气泡液体系统的能量

如果液体中只含有两个具有非球形扰动的气泡,气泡之间的距离为d.如图1所示,在声场作用下,假设每个气泡在形变扰动中,只出现单一形状模态的形变扰动,不存在不同形状模态的耦合,则气泡1(2)壁面离气泡1(2)中心的距离可表示为

式中,R1(t)和R2(t)是气泡1,2无扰动时的平均半径;an是气泡1壁面的第n阶形状模态振幅,bm是气泡2壁面的第m阶形状模态振幅,其中Ym分别是n,m阶球谐函数,在这里只考虑2的模式.考虑气泡之间的相互作用后,an,bm应该是关于气泡R1,R2和t的函数,如果假设液体为不可压缩的液体,根据Plesset提出的算法[10],可得气泡1,2周围液体的速度势可表示为

气泡系统的动能为

式中积分在整个流体体积进行,ρ为液体的密度,v为液体粒子的振速.液体粒子的速度可以表示为如下形式:

式中v1,v2分别为第1,2个气泡振动引起液体粒子的振动速度,可表示为

假设气泡1,2在整个振动过程中关于轴对称,联立(3),(4),(5),(6),(7)和(8)式,可得液体的动能为

其中,

气泡系统的势能等于周围液体对气泡做的功,即

式中

其中Pa为驱动压强,σ为表面张力系数,γ为绝热指数,P0为液体静压力,Pv为气泡内部的蒸汽压强,为气泡1发生形变后的2个主要曲率半径,其可表示为[10]

将(11)和(12)式代入(10)式中第一部分积分,可得在气泡1的势能,将(11)和(12)式中的下脚标1换成2按照同样的积分方法可得气泡2的势能.气泡1,2所组成的系统的总势能为

图1 声场中两个非球形气泡示意图Fig.1.Two nonspherical bubbles in an acoustic field.

3 含气泡液体系统的Lagrangian方程

气泡系统的无耗散的Lagrange方程为

式中L为Lagrangian函数,定义为系统的动能减去势能;Θ为广义坐标.若分别令两个气泡的径向半径和形状振幅为广义坐标,将(9)和(13)式代入(14)式,可得到两组方程:

方程(15)和(16)正好是考虑了气泡之间声相互作用的Rayleigh-Plesset(R-P)方程,即Doinikov模型方程,如果考虑气泡在振动过程中受到黏滞阻力的影响,可以在方程(15)和(16)等式右边中引入黏滞阻力项作为修正,其中µ为液体的黏滞阻力系数.

对比(19),(20)和(21)式可以看出,两个非球形气泡之间的声相互作用对于每个气泡形状模态振幅的影响体现在An和Bm因子上,表示两个气泡在振动过程中出现形状耦合振动,表示两个气泡无形状耦合振动.而和项则是流体力学观点和能量观点两种方法计算精度不同造成的,当两个气泡之间距离足够大时,气泡之间的相互作用忽略不计,两个气泡的径向振动方程(15)和(16)就回到了单气泡的R-P方程,形状模态振幅方程(19),(20)就回到单气泡形状模态振幅方程.

4 相互作用力对气泡形状不稳定性的影响

对于单个气泡而言,表面张力总是指向气泡内部,其效应是使得气泡能够保持最稳定的球形,另一方面,在声场作用下当气泡收缩到最小体积时,气泡内部的压强也非常大并且指向气泡外部,这两种强的且反方向的力导致了气泡表面的不稳定[9].当两个气泡相互接近时,两个振动的气泡会产生相互作用力,而气泡之间的作用力的方向则取决于两个气泡形状模态阶数、气泡的初始半径、气泡间距和驱动声场,正的相互作用力代表了扩张力,负的相互作用力代表了压缩力,任意时刻,由于气泡之间相互作用力使得具有非球形扰动的气泡表面的受力情况变得复杂,从而使得气泡形状不稳定性也受到一定的影响.

4.1 气泡耦合振动对气泡形状不稳定性的影响

对于具有非球形扰动的气泡而言,其在声场中振动有两种运动趋势:1)其形状模态振幅随时间呈现指数衰减,这种情况下气泡是稳定的;2)其形状模态振幅随时间呈现指数增大,这种情况下气泡越来越偏离球形,此时气泡属于形状不稳定,最终破裂.为了研究气泡之间的耦合振动对气泡形状不稳定性的影响,我们数值模拟了声场中具有相互作用的两个非球形气泡的形状不稳定性相图,模拟条件为:ρ=1000 kg/m3,µ =1×10−3kg/(m·s),σ=0.0725 N/m,γ=1.4,P0=1.05×105Pa.

图2(a)是当驱动频率为20 kHz,液体中只有一个非球形气泡时气泡的不稳定性相图,图2(b)是气泡1和气泡2做径向耦合振动时,气泡1的形状不稳定性相图,图2(c)是气泡1和气泡2做径向耦合振动和形状耦合振动时,气泡1的形状不稳定性相图.图中黑色区域代表气泡形状的稳定区域,白色区域代表气泡形状的不稳定区域.对比可以看出:一定条件下,气泡之间的形状耦合振动和径向耦合振动使得同样驱动条件下气泡的形状稳定区域有了一定程度的增加.换言之,气泡发生形状偏移时也会使得气泡之间产生相互作用,一定条件下,这种作用力能够减弱气泡的塌缩速度的剧烈程度,使得气泡对抗非对称性扰动的能力增加.

为了更直观地研究气泡之间的相互作用力对气泡形状不稳定性的影响,我们数值模拟了不同驱动条件下具有非球形扰动的气泡的2阶形状模态振幅随时间的变化关系.图3(a)—(c)分别是驱动压强振幅为1.265×105Pa、驱动频率为20 kHz时,初始半径为2µm的单个非球形气泡的二阶形状模态振幅随时间的变化关系曲线,初始半径为2µm的气泡1和相距200µm的初始半径为2µm的气泡2做径向耦合振动时,气泡1的2阶形状模态振幅随时间的变化曲线,初始半径为2µm的气泡1与相距200µm的初始半径为2µm的气泡2做形状耦合和径向耦合振动时,气泡1的2阶形状模态振幅随时间的变化曲线.对比可以发现,当气泡开始扩张时,气泡表面就已经激发了2阶形状振动模式,气泡不稳定性出现在气泡的收缩阶段,在上述模拟条件下,单个气泡的非球形扰动随着时间的增加越来越大,在20个驱动周期后,气泡表面越来越不稳定,最终导致气泡的崩溃破裂.当液体中存在两个气泡且两个气泡形状模态不同时,随着时间的增加,气泡1的2阶形状模态振幅越来越大,20个驱动周期后,也会使得气泡崩溃破裂,但相较于单个气泡而言,两个气泡之间径向耦合产生的相互作用力使得每个气泡的形状模态振幅减小.当两个气泡的形状模态相同时,对比可以发现,20个驱动周期后,气泡1的形状模态振幅只有气泡半径的0.02左右,且呈周期性变化,因此,气泡不会因形状扰动而发生破裂,气泡是稳定的,也就是说,一定条件下气泡间的形状耦合振动增加了气泡对抗非球形扰动的能力.

上述结论是建立在两个气泡相距较近的基础上,两个气泡之间的相互作用与两个气泡之间的距离有关,当两个气泡之间的距离超过一定临界距离时,两个气泡之间的相互作用非常微弱,可以近似为两个无耦合的自由气泡,此时便回到单气泡稳定振动问题上.通过数值模拟发现,对于形状耦合而言,小气泡的临界作用距离要大于大气泡的临界作用距离.除此之外,数值模拟结果还表明,当驱动声场的频率继续增加,到达两个气泡中某个气泡的共振频率之上,甚至超过某个气泡的共振频率时,气泡的稳定性区域较低频声场有了较大程度的提高,但气泡之间的形状耦合振动仍然在一定程度上增加了气泡对抗非球形扰动的能力,由于篇幅有限,在下面的讨论中仍然只讨论驱动频率为20 kHz的情况.

图2 气泡的2阶形状不稳定性相图(f=20 kHz)(a)单个气泡;(b)气泡1,,d=200µm,R20=2µm;(c)气泡1,n=m,R20=2µm,d=200µmFig.2.Phase diagrams for shape instability of a bubble with the second shape mode(f=20 kHz)(a)A single bubble;(b)the first bubble,,d=200µm and R20=2µm;(c)the first bubble,n=m,R20=2µm and d=200µm.

图3 气泡的2阶形状振幅随时间的变化(a)单个气泡(R0=2µm);(b)气泡1;(c)气泡1(n=m,R10=2µm)Fig.3.Amplitude of the second shape mode of a bubble with the time:(a)A single bubble(R0=2µm);(b)the first bubble;(c)the first bubble(n=m and R10=2µm).

4.2 气泡形状模态阶数对气泡形状不稳定性的影响

同样驱动条件下,初始形状模态对气泡的形状不稳定性亦会产生一定的影响.我们研究了驱动声场频率为20 kHz时,不同气泡形状模态阶数下气泡的形状不稳定性相图(n=m),图4(a)—(c)分别是n=2,3,4时气泡1做形状耦合和径向耦合振动(n=m,R20=2µm,d=200µm)时的形状不稳定性相图.对比可以发现:在本文的研究条件下,气泡1的2阶形变扰动最不稳定,随着形状模态阶数的增加,气泡1的形状稳定性区域也在增加,气泡1在3,4阶形状模态下的形状稳定性相较同样驱动条件下气泡处于2阶形状模态下的形状稳定性有了明显的提高.

为了更直观地表明上述现象,我们数值模拟了10个驱动周期下气泡1在不同形状模态下的形状振幅(n=m),如图5(a)—(c)所示,气泡1的初始半径为11µm,气泡2的初始半径为2µm,两个气泡之间的距离为200µm,驱动压强振幅为0.85×105Pa.对比可以看出:气泡1的2阶形状模态振幅在10个驱动周期内随着时间的变化逐渐增加,且远远超出了其半径,这表明,在这样的驱动条件下,气泡1的2阶形变扰动是不稳定的,形变扰动的结果是使得气泡1由于形状不稳定性而破裂;而同样的驱动条件下,同样初始半径的气泡1的3,4阶形变扰动是稳定的,在我们研究的10个驱动周期内,其形状模态振幅非常小,其中气泡1的3阶形状模态振幅只有气泡半径的0.04左右,且形变扰动具有周期性特征,气泡1的4阶形状模态振幅的大小为气泡1的半径的0.5左右,这说明,气泡1的3和4阶形状振幅在这样的驱动条件下是非常小的,且基本上稳定不变,不会使得气泡1由于形变而发生破裂,气泡1此时处于形状稳定区域.

图4 不同形变模态下气泡1的形状不稳定性相图(n=m,R20=2µm,d=200µm)(a)n=2;(b)n=3;(c)n=4Fig.4.Phase diagrams for shape instability of the first bubble with different shape mode order(n=m,R20=2µm and d=200µm):(a)n=2;(b)n=3;(c)n=4.

图5 气泡1的形状模态振幅随时间的变化(n=m,R10=11µm,R20=2µm,Pa=0.85×105Pa,d=200µm)(a)n=2;(b)n=3;(c)n=4Fig.5.Amplitude of the second shape mode of the first bubble with the time(n=m,R10=11µm,R20=2µm,Pa=0.85×105Pa and d=200µm):(a)n=2;(b)n=3;(c)n=4.

4.3 相邻气泡初始半径对气泡表面不稳定性的影响

当两个具有非球形扰动的气泡在声场中振动时,同样的声场驱动下,不同尺寸的气泡对之间的相互作用力不同,耦合振动方式也有所不同,因此对每个气泡的形状稳定性的影响也会不同,为了研究不同初始尺寸的相邻气泡对气泡1的形状不稳定性的影响,数值模拟了驱动频率为20 kHz下,两个气泡的形状模态相同(n=m)时,当气泡2的初始半径发生变化时气泡1的形状不稳定性相图.

我们发现,当气泡间距和驱动声场频率均保持不变时,气泡2初始半径的改变对气泡1的形状不稳定性也会产生一定的影响.图6是气泡1与不同初始半径的气泡2进行形状耦合和径向耦合振动(n=m)时的形状不稳定性相图,对比可以看出,气泡2的初始半径的改变会引起气泡1的形状不稳定性的明显改变.如图6(a)—(c)所示,在上述模拟条件下,随着气泡2的初始半径的增加,气泡1的形状不稳定性区域减小,且形状不稳定性区域的减小主要集中在小尺寸气泡范围内;对于大尺寸范围(4µm以上)内的气泡,随着气泡2初始半径的增加,其形状不稳定性区域不发生明显变化.发生这一现象的原因是,当大尺寸的气泡处于形状稳定性区域时,气泡2(初始半径从2µm增加到4µm)也处于形状稳定性区域,因此两个气泡的形状耦合和径向耦合振动所产生的作用力对于4µm以上气泡的形状稳定性影响不大,而对于小范围尺寸的气泡1(2—4µm),当其处于形状稳定性区域时,气泡2的形状不稳定性取决于其初始半径,如果此时气泡2处于不稳定区域,那么气泡2的塌缩破裂会对这个尺寸范围内的气泡1产生巨大的作用力,使得气泡1也迅速破裂,进而使得气泡1的形状不稳定性区域减小.根据上述原因,气泡1的形状稳定性区域会随着气泡2的初始半径的增加而减小,并且减小的区域主要位于与气泡2初始半径相近的区域范围.

图6 气泡1的形状不稳定性相图(n=m,f=20 kHz,d=200µm)(a)R20=2µm;(b)R20=3µm;(c)R20=4µmFig.6.Phase diagram for shape instability of the first bubble with different radius of the second bubble(n=m,f=20 kHz and d=200µm):(a)R20=2µm;(b)R20=3µm;(c)R20=4µm.

5 结 论

声场中的气泡,由于做周期性振动产生声场,因而和相邻气泡之间产生相互作用力,本文计算了两个气泡系统的能量,并基于Lagrange方程得到了存在声相互作用的气泡的动力学方程和形状稳定性方程,数值研究了声场中存在声相互作用的气泡的动力学及相互作用力对非球形气泡形状不稳定性的影响.研究结果表明两个存在非球形扰动的气泡之间存在着两种耦合模式:形状耦合模式和径向耦合模式.当两个气泡的非球形扰动模态相同(n=m)时,气泡之间同时存在着两种耦合模式,气泡之间的相互作用力由这两种耦合模式所产生;当两个气泡的非球形扰动模态不相同时,两个气泡之间只存在径向耦合模式,气泡之间的相互作用力由径向耦合模式产生.

研究发现,同样的驱动条件下,非球形气泡的形状不稳定性区域会随着气泡形状模态阶数的增加而减小.两个具有非球形扰动的气泡之间的相互作用力会对单个气泡的形状不稳定产生一定的影响,声场中气泡与比自身初始半径小或者等于自身初始半径的气泡的径向耦合振动和形状耦合振动能够在不同程度上减少气泡自身的形状不稳定性,反之则能够增加自身的形状不稳定性.

由于实际气泡的非球形扰动模态并非单一形状模态,而是存在不同形状模态的相互耦合,因此实际泡群振动中可能同时具有径向耦合振动和形状耦合振动这两种模式,而形状耦合振动的出现会显著地增加气泡之间的相互作用,改变空化气泡的形状不稳定,从而影响泡群中气泡的空化剧烈程度.

[1]Crum L A 1994 J.Acoust.Soc.Am.95 559

[2]Gaitan D F,Crum L A,Church C C,Roy R A 1992 J.Acoust.Soc.Am.91 3166

[3]Holot R G,Gaitan D F 1996 Phys.Rev.Lett.77 3791

[4]Zhang S G,Duncan J H 1994 Phys.Fluids 6 2352

[5]Leong T,Yasui K,Kato K,Harvie D,Ashokkumar M,Kentish S 2014 Phys.Rev.E 89 043007

[6]Pelekasis N A,Tsamopouslos J A 1990 Phys.Fluids A 2 1328

[7]Feng Z C,Leal L G 1997 Annu.Rev.Fluid Mech.29 201

[8]Reddy A J,Szeri A J 2002 Phys.Fluids 14 2216

[9]Harkin A A,Kaper T J,Nadim A 2013 Phys.Fluids 25 062101

[10]Plesset M S 1954 J.Appl.Phys.25 96

[11]Brenner M P,Lohse D,Dupont T F 1995 Phys.Rev.Lett.75 954

[12]Bogoyavlenskiy V A 2000 Phys.Rev.E 62 2158

[13]Wang W J,Chen W Z 2003 J.Acoust.Soc.Am.114 1898

[14]Liu H J,An Y 2003 Acta Phys.Sin.52 620(in Chinese)[刘海军,安宇 2003物理学报52 620]

[15]Qian M L,Cheng Q,Ge C Y 2002 Acta Acust.27 289(in Chinese)[钱梦騄,程茜,葛曹燕2002声学学报27 289]

[16]Hilgenfeldt S,Lohse D,Brenner M P 1996 Phys.Fluids 8 2808

[17]Godinez F A,Navarrete M 2011 Phys.Rev.E 84 016312

[18]Ueno I,Ando J,Koiwa Y,Saiki T,Kaneko T 2015 Eur.Phys.J.Special Topics 224 415

[19]Lu Y,Katz J,Prosperetti A 2013 Phys.Fluids 25 073301

[20]Ida M,Naoe T,Futakawa M 2007 Phys.Rev.E 76 046309

[21]Zhang W J,An Y 2013 Phys.Rev.E 87 053023

[22]Hens A,Biswas G,De S 2014 Phys.Fluids 26 012105