PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of procaThe Royal Society PublishingProceedings AAboutBrowse by SubjectAlertsFree Trial
 
Proc Math Phys Eng Sci. 2016 April; 472(2188): 20160031.
PMCID: PMC4892283

Lamb-type waves generated by a cylindrical bubble oscillating between two planar elastic walls

Abstract

The volume oscillation of a cylindrical bubble in a microfluidic channel with planar elastic walls is studied. Analytical solutions are found for the bulk scattered wave propagating in the fluid gap and the surface waves of Lamb-type propagating at the fluid–solid interfaces. This type of surface wave has not yet been described theoretically. A dispersion equation for the Lamb-type waves is derived, which allows one to evaluate the wave speed for different values of the channel height h. It is shown that for ht, where λt is the wavelength of the transverse wave in the walls, the speed of the Lamb-type waves decreases with decreasing h, while for h on the order of or greater than λt, their speed tends to the Scholte wave speed. The solutions for the wave fields in the elastic walls and in the fluid are derived using the Hankel transforms. Numerical simulations are carried out to study the effect of the surface waves on the dynamics of a bubble confined between two elastic walls. It is shown that its resonance frequency can be up to 50% higher than the resonance frequency of a similar bubble confined between two rigid walls.

Keywords: bubble, ultrasound, surface wave, solid–fluid interface, microfluidics

1. Introduction

This study is inspired by the use of bubbles in acoustically activated microfluidic systems [1]. An example of such a system is described by Rabaud et al. [2]: micrometre-size cylindrical bubbles are positioned in a fluid channel, 25 μm in depth, confined by elastic walls which are set into oscillation by a vibrating glass rod moulded in the upper wall. Under such conditions, the bubble dynamics exhibit a number of interesting effects, such as the generation of surface waves propagating along the solid–fluid interfaces and the self-organization of bubbles in crystal-like structures where equilibrium distances between bubbles are much smaller than distances that are predicted by the theory of secondary Bjerknes forces [3]. Another phenomenon of interest is acoustic microstreaming generated by bubbles in microfluidic devices [4]. This phenomenon is used for micromixing of fluids [5] and sorting of microparticles and cells [6,7]. A variety of experimental data on the dynamics of bubbles in narrow planar channels can be found in a series of recent papers [813], where both single bubbles and pairs of interacting bubbles are investigated.

Available theoretical models on the dynamics of cylindrical bubbles in planar channels consider only bulk scattered waves generated by the bubble oscillation in the fluid gap [1416]. These models are appropriate if the walls of the channel are rigid. There are no models that allow for surface waves, which are generated at the solid–fluid interfaces if the channel walls are elastic, whereas experimental observations show that surface waves can play a decisive role in the motion of the fluid layer [2].

In acoustics, several types of surface waves are encountered [17]. The surface waves relevant to our study are Rayleigh, Scholte and Lamb waves, so it is pertinent to remind their definitions. Rayleigh waves are surface waves that travel at an interface between an elastic solid and a vacuum. The existence of these waves was theoretically predicted by Lord Rayleigh in 1885 [18]. The speed of Rayleigh waves is slightly less than that of shear waves propagating in the body of an elastic solid. Scholte waves are surface waves that propagate at an interface between an elastic solid and a fluid. They are named after Scholte, who discovered them in 1947 [19]. Scholte waves are slower than Rayleigh waves. Lamb waves propagate in an elastic plate placed in a vacuum or a fluid [20]. The mathematical analysis of these waves was first developed by Lamb in 1917 [21]. Lamb waves are divided into symmetric and antisymmetric types depending on whether the motion of the elastic medium is symmetric or antisymmetric about the middle plane of the elastic plate.

Although the type of surface waves considered in our study is intimately related to the types described above, it is different from them. In a certain sense, the case under study is opposite to the case considered by Lamb because we deal with surface waves propagating in a fluid channel embedded in an elastic solid. To the best of our knowledge, this case has not been considered as yet. Our study provides a dispersion equation for this type of surface waves, which allows one to evaluate the speed of the surface waves for different values of the channel height. The ultimate purpose of our study is to develop a theoretical model for the volume oscillation of a cylindrical bubble placed between two planar elastic walls, allowing for both bulk and surface waves generated in this system. It should be emphasized that we deal with a special case of bubble confinement. Indeed, the channel height is so small with respect to the bubble size that the bubble always remains in contact with the walls, which leads to a shape close to a cylinder, as depicted in figure 1. Therefore, there is no possibility for the bubble to take the spherical shape or be at a distance from the walls. This fact is a result of the design of microfluidic devices of our interest [2], which makes our case different from those implied by Ilinskii et al. [15] and Leighton et al. [2224].

Figure 1.
A cylindrical bubble in a fluid channel between two elastic plane walls. (Online version in colour.)

2. Problem formulation and solution

Let us consider a cylindrical gas bubble located in a fluid channel of height h between two planar elastic walls (figure 1). The bubble undergoes volume oscillation in response to an imposed acoustic pressure field. The bubble oscillation is considered to be linear and dependent on time as exp(−iωt). For simplicity, the time factor will be omitted from equations. Quantities related to the upper and lower walls will be denoted by subscripts 1 and 2, respectively. The action of the bubble on the fluid–solid interfaces is modelled as a normal harmonic point load. Because the problem under consideration is axially symmetric, cylindrical coordinates are used, whose origin is located in the middle of the bubble and the z-axis is perpendicular to the surface of the walls, as shown in figure 1. It should be emphasized that the wave field in the fluid channel consists of two components. One component is the bulk scattered wave caused by the bubble radial oscillation. The second component is caused by the surface waves that propagate in the elastic walls and which, in the process of this propagation, induce perturbations in the fluid. This component will be called surface waves in the fluid channel.

The calculation consists of the following stages. The equations for the surface waves in the elastic walls and in the fluid channel are derived in §2a and §2b, respectively. The derivation is based on Hankel transforms. The boundary conditions on the solid–fluid interfaces are obtained in §2c. The dispersion equation for the surface wave speed is considered in §2d. The resulting expressions for the surface waves in the space domain are calculated by the inverse Hankel transforms in §2e. The equations for the bulk scattered wave in the fluid channel are given in §2f. The boundary conditions on the side bubble surface are obtained in §2g. The equation of bubble oscillation is derived in §2h.

(a) Waves in the elastic walls

As shown in figure 1, the walls occupy the space with |z|>h/2. The displacement vector of the jth wall (j=1,2) can be written as

uj=φj+×ψj,
2.1

with [var phi]j and ψj being called the scalar and vector potentials, respectively. Equation (2.1) is used as the most general mathematical representation for an arbitrary vector field (e.g. [17] or [25]). In the view of axial symmetry, the dependence on the angular coordinate θ is absent and hence [var phi]j and ψj can be written as [var phi]j=[var phi]j(r,z) and ψj=ψj(r,z)eθ, respectively, where eθ is the unit azimuth vector. The function [var phi]j and ψj obey the following equations [25]:

1rr(rφjr)+2φjz2+kl2φj=0
2.2

and

1rr(rψjr)+2ψjz2+(kt21r2)ψj=0,
2.3

where kl=ω/cl and kt=ω/ct are the wavenumbers of the longitudinal and transverse waves in the solid, respectively. The speeds of these waves, cl and ct, are given by [17]

cl=λ+2μρsandct=μρs,
2.4

where λ and μ are the Lamé coefficients of the walls (μ is also called shear modulus) and ρs is the wall density. The components of the displacement vector uj are calculated by [25]

ujr=φjrψjzandujz=φjz+1r(rψj)r.
2.5

Solutions to equations (2.2) and (2.3) can be found by applying Hankel transforms with respect to the radial variable r. The Hankel transform of order n of a function f(r) is defined as [26]

fHn(k)=0f(r)Jn(kr)rdr,
2.6

where k is the parameter of the Hankel transform and Jn is the Bessel function of the first kind of order n. The inverse transform is given by

f(r)=0fHn(k)Jn(kr)kdk.
2.7

Applying the Hankel transforms of zero and first orders to equations (2.2) and (2.3), respectively, one obtains

d2φjH0dz2ql2φjH0=0
2.8

and

d2ψjH1dz2qt2ψjH1=0,
2.9

where

ql2=k2kl2andqt2=k2kt2.
2.10

Solutions of equations (2.8) and (2.9) which remain finite for large values of |z| are

φ1H0=Φ1(k)eql(zh/2),
2.11

ψ1H1=Ψ1(k)eqt(zh/2),
2.12

φ2H0=Φ2(k)eql(z+h/2)
2.13

andψ2H1=Ψ2(k)eqt(z+h/2),
2.14

where Φj(k) and Ψj(k) are functions to be determined from the boundary conditions on the solid–fluid interfaces at z = ±h/2.

Application of the Hankel transform to equations (2.5) yields

ujrH1=kφjH0dψjH1dzandujzH0=dφjH0dz+kψjH1,
2.15

Substituting equations (2.11)–(2.14) into equations (2.15), one finds

u1rH1=kΦ1(k)eql(zh/2)+qtΨ1(k)eqt(zh/2),
2.16

u1zH0=qlΦ1(k)eql(zh/2)+kΨ1(k)eqt(zh/2),
2.17

u2rH1=kΦ2(k)eql(z+h/2)qtΨ2(k)eqt(z+h/2)
2.18

andu2zH0=qlΦ2(k)eql(z+h/2)+kΨ2(k)eqt(z+h/2).
2.19

(b) Surface waves in the fluid channel

The fluid area corresponds to −h/2≤zh/2. The displacement vector in the fluid can be written as

uf=φf,
2.20

where the potential [var phi]f obeys the equation [27]

1rr(rφfr)+2φfz2+kf2φf=0,
2.21

in which kf = ω/cf and cf is the speed of sound in the fluid. The fluid velocity corresponding to uf is vf=−iωuf. The components of uf are calculated by

ufr=φfrandufz=φfz.
2.22

Applying the Hankel transform of zero order to equation (2.21), one obtains

d2φfH0dz2qf2φfH0=0,
2.23

where

qf2=k2kf2.
2.24

The general solution of equation (2.23) is given by

φfH0=Φf1(k)eqf(zh/2)+Φf2(k)eqf(z+h/2),
2.25

where Φf1(k) and Φf2(k) are functions to be determined from the boundary conditions at z = ±h/2.

Application of the Hankel transform to equations (2.22) yields

ufrH1=kφfH0andufzH0=dφfH0dz.
2.26

Substituting equation (2.25) into equations (2.26), one finds

ufrH1=k[Φf1(k)eqf(zh/2)+Φf2(k)eqf(z+h/2)]
2.27

and

ufzH0=qf[Φf1(k)eqf(zh/2)Φf2(k)eqf(z+h/2)].
2.28

(c) Boundary conditions on the solid–fluid interfaces

The boundary conditions on the surfaces of the elastic walls are written as follows [25,28]:

u1z=ufzatz=h2,u2z=ufzatz=h2,
2.29

T1rz=0atz=h2,T2rz=0atz=h2
2.30

andT1zz=pfQδ(r)2πratz=h2,T2zz=pfQδ(r)2πratz=h2,
2.31

where Tjrz and Tjzz are the shear and normal stresses in the walls, pf is the fluid pressure corresponding to uf, Q is the magnitude of the point load produced by the bubble at the points with r=0 and zh/2 and δ(r) is the Dirac delta function. Equations (2.29) follow from the continuity of the normal displacement at the solid–fluid interfaces, equations (2.30) mean that the shear stress vanishes at the interfaces and equations (2.31) describe the normal stress balance at the interfaces.

The expressions for Tjrz, Tjzz and pf are given by [25,27]

Tjrz=μ(ujrz+ujzr),
2.32

Tjzz=(λ+2μ)ujzz+λr(rujr)r
2.33

andpf=ρfω2φf,
2.34

where ρf is the fluid density. Applying the Hankel transforms to equations (2.29)–(2.34), one obtains

u1zH0=ufzH0atz=h2,u2zH0=ufzH0atz=h2,
2.35

T1rzH1=0atz=h2,T2rzH1=0atz=h2,
2.36

T1zzH0=pfH0Q2πatz=h2,T2zzH0=pfH0Q2πatz=h2,
2.37

TjrzH1=μ[2kdφjH0dz+(2k2kt2)ψjH1],
2.38

TjzzH0=μ[(2k2kt2)φjH0+2kdψjH1dz]
2.39

andpfH0=ρfω2φfH0.
2.40

Substitution of equations (2.11)–(2.14), (2.17), (2.19), (2.25), (2.28) and (2.38)–(2.40) into equations (2.35)–(2.37) yields the following equations for the functions Φj(k), Ψj(k), Φf1(k) and Φf2(k):

kΨ1qlΦ1=qf(Φf1Φf2eqfh),
2.41

kΨ2+qlΦ2=qf(Φf1eqfhΦf2),
2.42

2kqlΦ1(2k2kt2)Ψ1=0,
2.43

2kqlΦ2+(2k2kt2)Ψ2=0,
2.44

(2k2kt2)Φ12kqtΨ1=ρfω2μ(Φf1+Φf2eqfh)Q2πμ
2.45

and(2k2kt2)Φ2+2kqtΨ2=ρfω2μ(Φf1eqfh+Φf2)Q2πμ.
2.46

The waves produced by the bubble in the elastic walls are symmetric, because the walls are assumed to be identical and the action of the bubble on them is symmetric about the point z=0. Therefore, the following conditions should be met [17]:

u1rH1(z=h2)=u2rH1(z=h2)andu1zH0(z=h2)=u2zH0(z=h2).
2.47

From equations (2.16)–(2.19), it follows that conditions (2.47) are met if

Φ2(k)=Φ1(k)andΨ2(k)=Ψ1(k).
2.48

Substitution of equations (2.48) into equations (2.41)–(2.46) shows that

Φf2(k)=Φf1(k)=Φf(k),
2.49

and that only three of these equations remain independent. As a result, we obtain the following set:

qlΦ1kΨ1+qf(1eqfh)Φf=0,
2.50

2kqlΦ1(2k2kt2)Ψ1=0
2.51

and(2k2kt2)Φ12kqtΨ1+ρfω2μ(1+eqfh)Φf=Q2πμ.
2.52

The solutions of this set are

Φ1(k)=2k2kt2D(k)Q2πμ,
2.53

Ψ1(k)=2kqlD(k)Q2πμ
2.54

andΦf(k)=kt2qlqf(1eqfh)D(k)Q2πμ,
2.55

where

D(k)=(2k2kt2)24k2qlqt+ρfkt4qlρsqf1+eqfh1eqfh.
2.56

(d) Dispersion equation for surface waves

The theoretical analysis of Rayleigh waves produced by a point load, based on the use of Hankel integral transforms, is given in a book by Achenbach [25]. An analogous consideration of Scholte waves is presented by Zhu et al. [28]. These works show that at large distances from the loading point, where general solutions similar to equations (2.53) and (2.54) turn into far-field expressions giving the Rayleigh and Scholte waves, the dispersion equations for these waves are identical to the denominators of the general solutions similar to equations (2.53) and (2.54). This follows from the fact that a vanishing denominator gives rise to a dominant term in the spectrum of emitted waves. Analogously, in our case the dispersion equation for the surface waves of our interest is given by

D(k)=(2k2kt2)24k2qlqt+ρfkt4qlρsqf1+eqfh1eqfh=0.
2.57

The roots of this equation are the wavenumbers of Lamb-type symmetric waves that propagate in the microfluidic system under study. If the fluid density tends to zero, ρf→0, equation (2.57) reduces to the equation for Rayleigh waves [17,25], while for infinite channel height, h, it reduces to the equation for Scholte waves [17,28]. Equation (2.57) allows one to evaluate the surface wave speed as a function of h. Examples of such evaluations are presented in §3a. It is noteworthy that equation (2.57) is derived for the first time so it can be considered as one of the main results of this study.

(e) Expressions for u1, uf and pf in the space domain

Applying the inverse Hankel transform to equations (2.16), (2.17), (2.27), (2.28) and (2.40) and using equations (2.53)–(2.55), one obtains

u1r=Qkt2πμIsr(r,z),
2.58

u1z=Qkt2πμIsz(r,z),
2.59

ufr=QktπμIfr(r,z),
2.60

ufz=QktπμIfz(r,z)
2.61

andpf=Qρfω2πμIp(r,z),
2.62

where

Isr(r,z)=01D~(ξ)[(2ξ21)eq~lkt(zh/2)2q~lq~teq~tkt(zh/2)]J1(ξktr)ξ2dξ,
2.63

Isz(r,z)=0q~lD~(ξ)[(2ξ21)eq~lkt(zh/2)2ξ2eq~tkt(zh/2)]J0(ξktr)ξdξ,
2.64

Ifr(r,z)=0q~leq~fkth/2cosh(q~fktz)q~fD~(ξ)(1eq~fkth)J1(ξktr)ξ2dξ
2.65

Ifz(r,z)=0q~leq~fkth/2sinh(q~fktz)D~(ξ)(1eq~fkth)J0(ξktr)ξdξ
2.66

andIp(r,z)=0q~leq~fkth/2cosh(q~fktz)q~fD~(ξ)(1eq~fkth)J0(ξktr)ξdξ.
2.67

Equations (2.63)–(2.67) are expressed in terms of normalized quantities that are defined as follows:

ξ=kkt=ctc,ξf=kfkt=ctcf,ξl=klkt=ctcl,
2.68

q~f=ξ2ξf2,q~l=ξ2ξl2,q~t=ξ21
2.69

andD~(ξ)=(2ξ21)24ξ2q~lq~t+ρfq~lρsq~f1+eq~fkth1eq~fkth,
2.70

where c stands for the speed that corresponds to the wavenumber k. The normalization by kt provides the simplest form of the equations. The method of numerical calculation of integrals (2.63)–(2.67) is described in appendix A.

(f) Bulk wave in the fluid channel

According to Ilinskii et al. [15], the bulk scattered wave generated by the bubble oscillation in the fluid channel can be represented by

vb=φb,
2.71

where vb is the fluid velocity and the velocity potential [var phi]b is given by the equation

φb=AH0(1)(kfr),
2.72

in which H0(1) is the Hankel function of the first kind of order zero [29] and A is a constant to be determined from the boundary conditions on the bubble surface. Equations (2.71) and (2.72) satisfy the wave equation for the fluid motion in the case of cylindrical symmetry. The fluid pressure corresponding to vb is calculated by

pb=iωρfφb.
2.73

It should be mentioned that the wavelength of the bulk wave, λf=2π/kf, is large compared with the wavelength of the surface waves as well as to the dimensions of the microfluidic system under study. This means that the spatial scale of the variations of pb is much greater than that of the perturbations produced by the surface waves. For this reason, pb does not contribute to the above perturbations and hence does not appear in equation (2.31).

(g) Boundary conditions on the bubble surface

Let us assume that the time-varying bubble radius can be represented as

R(t)=R0(1+aeiωt),
2.74

where R0 is the equilibrium bubble radius. Considering that the net fluid velocity is the sum of the velocities of the bulk and surface waves, the boundary condition for the normal velocity at the side bubble surface can be written as

vbr+vfr=iωR0aatr=R0,
2.75

where vbr is the radial component of vb,vfr is the radial component of the surface wave velocity, given by vfr=−iωufr, and the right-hand term is the amplitude of the velocity of the bubble oscillation. Assuming that h is small compared with the transverse wavelength λt=2π/kt, the value of ufr can be taken at z = 0; see equations (2.60) and (2.65). Substitution of equations (2.60), (2.71) and (2.72) into equation (2.75) yields

A=iωkfH1(1)(kfR0)[R0aktQπμIfr(R0,0)].
2.76

The pressure balance at the bubble surface is defined as

pg=P0+pac+pst+(pb+pf)eiωtatr=R0,
2.77

where pg is the gas pressure within the bubble, P0 is the hydrostatic pressure in the fluid, pac is the driving acoustic pressure and pst is the surface tension pressure. Calculating the bubble volume as

V=πR2hπR02h(1+2aeiωt)=V0(1+2aeiωt),
2.78

one obtains

pg=Pg0(V0V)γPg0(12γaeiωt),
2.79

where Pg0 is the equilibrium gas pressure and γ is the ratio of specific heats of the gas. The acoustic pressure pac is specified by

pac=Paeiωt,
2.80

where Pa denotes the acoustic pressure amplitude. For pst, one has

pst=σfRσfR0(1aeiωt),
2.81

where σf is the surface tension coefficient for the fluid–gas interface. Substituting equations (2.62), (2.72), (2.73) and (2.79)–(2.81) into equation (2.77) and separating constant and time-varying terms, one obtains

Pg0=P0+σfR0
2.82

and

a(2γPg0σfR0)=ρfω2QπμIp(R0,0)iωρfAH0(1)(kfR0)Pa.
2.83

Substituting equations (2.76) and (2.82) into equation (2.83) and using the second of equations (2.4) to eliminate μ, one finds the following expression for a:

a=(Pa+C)B,
2.84

where

B=2γP0+(2γ1)σfR0ρfω2R02H0(1)(αf)αfH1(1)(αf)
2.85

and

C=ρfkt2Qπρs[Ip(R0,0)αtH0(1)(αf)αfH1(1)(αf)Ifr(R0,0)],
2.86

αf=kfR0 and αt=ktR0.

(h) Equation for bubble oscillation

To derive an equation for the amplitude of the bubble oscillation a, we need an expression for the load magnitude Q in terms of a. Clearly, Q should be proportional to the variation of the gas pressure in the bubble, given by equation (2.79). It is admissible to assume that Q is also proportional to the area of the lateral section of the bubble, πR02. Therefore, Q can be represented as

Q=2πγPg0R02εa,
2.87

where ε is a proportionality coefficient that allows for the fact that the effective area that describes the action of the bubble on the wall can be different from πR02. In addition, ε can be a complex quantity because a phase shift can exist between the variation of the gas pressure and the reaction of the walls to it. This occurs because the walls may behave as a viscoelastic material with complex elastic and shear moduli, which leads to a phase lag between stress and strain [30].

Substitution of equation (2.87) into equations (2.84) and (2.86) yields

a=PaB+S,
2.88

where

S=2γPg0ρfαt2ερs[Ip(R0,0)αtH0(1)(αf)αfH1(1)(αf)Ifr(R0,0)].
2.89

Equation (2.88) allows one to calculate the response of the bubble as a function of ε for different values of the material and acoustic parameters. It should be mentioned that for ε=0, equation (2.88) reduces to the result obtained by Ilinskii et al. [15] for a cylindrical bubble between two rigid walls.

3. Numerical simulations and experimental verification

(a) Speed of surface waves

The speed of the surface waves propagating at the solid–fluid interfaces is given by solutions of equation (2.57). Numerical evaluations were performed for the following values of the physical parameters: Young’s modulus E=1.6 MPa, Poisson’s ratio σ=0.499,ρs=970 kg m−3,cf=1481 m s−1 and ρf=998 kg m−3. These values are typical of the elastic walls of a microfluidic channel made of polydimethylsiloxane (PDMS) and filled with water. The transverse wave speed for these parameters is ct=23.456 m s−3. This quantity sets the upper limit for the speed of surface waves [17]. The evaluations show that equation (2.57) has one real root which gives a frequency-dependent speed cs<ct. Figure 2 exemplifies the frequency dependence of cs at h=25 μm. The dependence of cs on h is demonstrated by figure 3. Note that the values of cs are normalized by ct, and those of h, by λt=ct/f, where f is the driving frequency. With this normalization, the plot shown in figure 3 remains the same for all frequencies. Figure 3 reveals that cs increases with increasing h and tends to the speed of the Scholte wave, cSch, when h tends to λt. However, at small h, the presence of the second boundary makes cs much smaller than cSch.

Figure 2.
Surface wave speed as a function of frequency at h=25 μm.
Figure 3.
Surface wave speed as a function of h.

(b) Solutions at small and large distances from the bubble

The theory of surface waves predicts that the behaviour of integrals (2.63)–(2.67) should coincide with the behaviour of the Hankel functions of zero and first orders at large distances from the loading point [17]. The calculation of the integrals by the method described in appendix A conforms to this expectation. As an example, figure 4 compares the behaviour of the integral Isz, which describes the vertical displacement of wall 1, and the Hankel function H0(1)(ksr), where ks=ω/cs is the wavenumber of the surface wave, given by equation (2.57). The calculations were made at f=30 kHz and h=25 μm, the material parameters being the same as in the preceding subsection. Figure 4a shows Re{Isz(r,h/2)} (solid) and Re{iH0(1)(ksr)} (dashed) as functions of radial distance r which is normalized by λs=2π/ks. The phases of these functions are compared in figure 4b. Note that the vertical lines in this figure do nothing but show that the sign of the phase is changed by jump. The comparison reveals that the functions give the same result for distances exceeding 0.4λs. However, for smaller distances, there is a considerable difference in their behaviour. Therefore, the approximation by the Hankel function is not suitable for use in equations (2.88) and (2.89) which describe the bubble oscillation.

Figure 4.
(a,b) Comparison of the integral Isz(r,h/2) (solid line) and the Hankel function H0(1)(ksr) (dashed line).

(c) Resonance curves

Comparison with experimental measurements (presented in the next subsection) shows that agreement with the theoretical model is achieved if the parameter ε is considered as a complex quantity whose phase is proportional to the bubble radius R0. Therefore, examples of resonance curves were calculated assuming that

ε=ε1exp(iε2ksR0),
3.1

where ε1 and ε2 serve as variation parameters.

Resonance curves shown in figure 5 were calculated by equation (2.88) for h=25 μm and f=30 kHz, varying R0 in the range from 16 to 83 μm. The material parameters for the solid and the fluid were the same as in the preceding subsections. The parameters related to the gas behaviour and the acoustic excitation were as follows: σf=0.072 N m−1, γ=1.4, P0=101.3 kPa and Pa=10 kPa. Figure 5a shows resonance curves for increasing values of ε1, assuming ε2=0. It is seen that at ε1=0 (no surface waves), the resonance peak is at R0f=1.1 m s−1, which corresponds to the result obtained by Ilinskii et al. [15]. When ε1 is increased, the resonance peak shifts to larger values of R0 and increases in magnitude. At ε1=1, the resonance peak is found to be at R0f=1.5 m s−1 and its magnitude is maximal. Note that for ε1=1, the effective surface area by which the bubble acts on the wall is equal to the area of the lateral section of the bubble. Further increase of ε1 decreases the magnitude of the resonance peak, moving its position to larger values of R0f, which, however, do not exceed 1.7 m s−1.

Figure 5.
Examples of resonance curves given by the theoretical model for different values of the parameters ε1 and ε2. (a) ε1=0,0.5,0.8,1,1.5,1.75,2; ε2=0. (b) ε1=1; ε2=0,−0.1,−0.2,−0.4,−0.8. ...

Figure 5b shows the evolution of resonance curves with increasing |ε2| at ε1=1. Increasing |ε2| can be related to an increase of damping in the wall material, which leads to flattening of resonance curves and decreasing resonance bubble radius. This behaviour is observed if ε2 is negative.

(d) Experimental verification

The experimental data shown by circles in figures 68 were adopted from work by Mekki-Berrada et al. [31], where they were acquired by a microfluidic set-up described by Rabaud et al. [2]. This set-up has been already mentioned in the Introduction. The experimental data show the normalized amplitude of bubble oscillation for cylindrical bubbles of different radii, placed in a fluid channel, 25 μm in depth, confined by elastic walls made of PDMS. The bubble oscillation is induced by a vibrating glass rod moulded in the upper channel wall about 150 μm above the solid–fluid interface. The experimental data were obtained at three values of the driving frequency: 30, 40 and 50 kHz. For each frequency, the measurements were made for a certain range of bubble radii, which was dictated by experimental conditions.

Figure 7.
(a,b) Comparison of experimental (circles) and theoretical (solid line) results for bubble oscillation at 40 kHz.

Figure 6.
(a,b) Comparison of experimental (circles) and theoretical (solid line) results for bubble oscillation at 30 kHz.
Figure 8.
(a,b) Comparison of experimental (circles) and theoretical (solid line) results for bubble oscillation at 50 kHz.

The solid lines in figures figures668 show the fitting of the experimental points by equation (2.88) with the fitting parameter ε given by equation (3.1). The material parameters were taken as pointed out in the preceding subsections. The theoretical curves were calculated at the following values of ε1 and ε2: for 30 kHz, ε1=1.76,ε2=−0.451; for 40 kHz, ε1=2.15,ε2=−0.465; for 50 kHz, ε1=1.95,ε2=−0.756. The parts (a) of the figures show fitting in the experimental measurement range and the parts (b) show the view of the full theoretical curve. The normalization of the data is performed by a maximal value in a depicted range. For this reason, the normalization is different in figure 8a,b. The experimental data indicate that there is strong damping in the PDMS walls since the resonance curves are very flat. As one can see, agreement between the theoretical and experimental data is achieved if the values of the fitting parameters are changed with frequency. This result suggests that both effective surface area of the bubble–wall contact and the damping characteristics of the PDMS walls depend on frequency by a more complicated way than that specified by equation (3.1). This problem calls for further consideration.

4. Conclusion

A theoretical model has been developed that describes the volume oscillation of an ultrasound-activated cylindrical bubble situated in a fluid channel between two planar elastic walls. The model includes both the bulk scattered wave, which propagates in the fluid gap, and the surface waves of Lamb-type, which propagate at the solid–fluid interfaces. The force exerted by the bubble oscillation on the channel walls, which is the source of the Lamb-type waves, was treated as a normal harmonic point load. A dispersion equation for the above surface waves was derived, which allows one to examine the dependence of the speed of these waves on the channel height h. Such an equation was calculated for the first time, which makes it one of the main results of the present study. It shows that for ht, where λt is the wavelength of the transverse wave in the elastic walls, the speed of the Lamb-type waves decreases with decreasing h, while for h→λt, their speed tends to the Scholte wave speed.

The solutions for the wave fields produced by the Lamb-type waves in the elastic walls and in the fluid channel were obtained using the Hankel transforms and expressed in terms of improper integrals. A method was proposed for numerically evaluating the integral solutions. Numerical simulations were carried out to examine the effect of the surface waves on the bubble dynamics. It was shown that the resonance frequency of a cylindrical bubble confined between two identical elastic walls could be up to 50% higher than the resonance frequency of a similar bubble confined between two rigid walls.

Experimental verification of the theoretical model has been carried out using measured values of the bubble oscillation amplitude, which were obtained for bubbles of different radii at three values of the driving frequency. The fitting of experimental data was performed by using two parameters. The first parameter describes the magnitude of the load exerted on the channel walls by the variation of the bubble gas pressure. The second parameter describes the phase shift between the variation of the bubble gas pressure and the reaction of the walls. It was found that agreement between the theoretical and experimental data is achieved on condition that the fitting parameters are considered frequency-dependent quantities.

Appendix A. Calculation of the integral expressions

The denominator functions of the improper integrals (2.63)–(2.67) have zeros, which makes direct numerical integration impossible. This difficulty can be overcome by using contour integration and the residue theorem as proposed by Vrettos [32]. The calculation method will be demonstrated here with integrals (2.65) and (2.67) as these integrals are used in the equation of bubble oscillation.

Expressing the Bessel function J0 as a sum of the Hankel functions of the first and second kind

2J0(ξktr)=H0(1)(ξktr)+H0(2)(ξktr),
A 1

substituting (A1) into equation (2.67) and replacing the real variable ξ with the complex variable s=ξ+, one obtains

2Ip(r,z)=I1+I2,
A 2

where

Ij=0F(s,z)H0(j)(sktr)sds
A 3

and

F(s,z)=q~leq~fkth/2cosh(q~fktz)q~fD~(s)(1eq~fkth).
A 4

Note that the variable ξ is also replaced with s in the radical functions q~f,q~l and q~t. As a result, these functions become complex valued and branch points occur at s=ξf, ξl,1, where ξf<ξl<1 because it is assumed that cf>cl>ct; see equations (2.68). The pole of the integrand F(s,z) is given by the root of the equation D~(ξ)=0. This root, which will be denoted by ξs, corresponds to a Scholte-type wave and ξs>1 as the speed of the Scholte wave is smaller than ct; see the first of equations (2.68). To make the integrand single valued and take into account the presence of the pole, the integration contours Γ1 and Γ2 are taken as shown in figure 9.

Figure 9.

An external file that holds a picture, illustration, etc.
Object name is rspa20160031-g9.jpg

Integration contours in the complex wavenumber plane.

Integrating I1 along the contour Γ1 and using the fact that H0(1)(sktr) vanishes along the infinite arc, one obtains

I1=I11+I12+I13+I14+2πiResp(ξs),
A 5

where

I11=0F+(iτ,z)H0(1)(iτktr)τdτ,
A 6

I12=0ξf[F+(ξ,z)F(ξ,z)]H0(1)(ξktr)ξdξ,
A 7

I13=ξfξl[F+(ξ,z)F(ξ,z)]H0(1)(ξktr)ξdξ
A 8

andI14=ξl1[F+(ξ,z)F(ξ,z)]H0(1)(ξktr)ξdξ,
A 9

the superscripts + and − denote that the function F(s,z) is calculated for the positive and negative values of the imaginary part of the radicals q~f,q~l and q~t, respectively, and Resp(ξs) is the residue at the pole ξs which is calculated by [33]

Resp(ξs)=[q~leq~fkth/2cosh(q~fktz)ξH0(1)(ξktr)[q~fD~(ξ)(1eq~fkth)]/ξ]ξ=ξs.
A 10

Integrating I2 along the contour Γ2 and using the fact that H0(2)(sktr) vanishes along the infinite arc, one obtains

I2=0F(iτ,z)H0(2)(iτktr)τdτ.
A 11

Equations (A 7)–(A 9) can be simplified using the fact that F(ξ,z)=[F+(ξ,z)]*, where the asterisk denotes the complex conjugate. Further, the sum of the integrals I11 and I2 can be transformed using the following identities:

F(iτ,z)=[F+(iτ,z)],
A 12

H0(1)(ix)=2πiK0(x)
A 13

andH0(2)(ix)=2πiK0(x),
A 14

where Kn(x) is the modified Bessel function of the second kind of order n [29]. As a result of these operations, the final expression for the integral Ip(r,z) is found to be

Ip(r,z)=2π0Im{F+(iτ,z)}K0(τktr)τdτi0ξfIm{F+(ξ,z)}H0(1)(ξktr)ξdξiξfξlIm{F+(ξ,z)}H0(1)(ξktr)ξdξiξl1Im{F+(ξ,z)}H0(1)(ξktr)ξdξ+πiResp(ξs),
A 15

where Im{x} means the imaginary part of x. These integrals are nonsingular and convergent and hence they can be calculated numerically without any difficulty.

Application of the same method to integral (2.65) yields

Ifr(r,z)=2π0Im{F+(iτ,z)}K1(τktr)τ2dτi0ξfIm{F+(ξ,z)}H1(1)(ξktr)ξ2dξiξfξlIm{F+(ξ,z)}H1(1)(ξktr)ξ2dξiξl1Im{F+(ξ,z)}H1(1)(ξktr)ξ2dξ+πiResfr(ξs),
A 16

where the residue at the pole ξs is calculated by

Resfr(ξs)=[q~leq~fkth/2cosh(q~fktz)ξ2H1(1)(ξktr)[q~fD~(ξ)(1eq~fkth)]/ξ]ξ=ξs.
A 17

Authors' contributions

A.D. carried out the development of the theory and drafted the manuscript; F.M.B. carried out the experimental measurements and helped draft the manuscript; P.T. participated in the experimental measurements and the interpretation of the results and helped draft the manuscript; P.M. participated in the development of the theory and the experimental measurements, coordinated the study and helped draft the manuscript. All authors gave final approval for publication.

Competing interests

We have no competing interests.

Funding

This research has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007–2013)/ERC grant agreement no. 614655 ‘Bubbleboost’.

References

1. Friend J, Yeo L 2011. Microscale acoustofluidics: microfluidics driven via acoustics and ultrasonics. Rev. Mod. Phys. 83, 647–704. (doi:10.1103/RevModPhys.83.647)
2. Rabaud D, Thibault P, Mathieu M, Marmottant P 2011. Acoustically bound microfluidic bubble crystals. Phys. Rev. Lett. 106, 134501 (doi:10.1103/PhysRevLett.106.134501) [PubMed]
3. Doinikov AA. 2005. Bjerknes forces and translational bubble dynamics. In Bubble and particle dynamics in acoustic fields: modern trends and applications (ed. Doinikov AA, editor. ), pp. 95–143. Kerala, India: Research Signpost.
4. Marmottant P, Raven JP, Gardeniers H, Bomer JG, Hilgenfeldt S 2006. Microfluidics with ultrasound-driven bubbles. J. Fluid Mech. 568, 109–118. (doi:10.1017/S0022112006002746)
5. Ahmed D, Mao X, Shi J, Juluri BK, Huang TJ 2009. A millisecond micromixer via single-bubble-based acoustic streaming. Lab Chip 9, 2738–2741. (doi:10.1039/b903687c) [PubMed]
6. Wang C, Jalikop SV, Hilgenfeldt S 2011. Size-sensitive sorting of microparticles through control of flow geometry. Appl. Phys. Lett. 99, 034101 (doi:10.1063/1.3610940)
7. Franke T, Braunmüller S, Schmid L, Wixforth A, Weitz D 2010. Surface acoustic wave actuated cell sorting (SAWACS). Lab Chip 10, 789–794. (doi:10.1039/b915522h) [PubMed]
8. Zwaan E, Le Gac S, Tsuji K, Ohl C-D 2007. Controlled cavitation in microfluidic systems. Phys. Rev. Lett. 98, 254501 (doi:10.1103/PhysRevLett.98.254501) [PubMed]
9. Quinto-Su PA, Lim KY, Ohl C-D 2009. Cavitation bubble dynamics in microfluidic gaps of variable height. Phys. Rev. E 80, 047301 (doi:10.1103/PhysRevE.80.047301) [PubMed]
10. Quinto-Su PA, Ohl C-D 2009. Interaction between two laser-induced cavitation bubbles in a quasi-two-dimensional geometry. J. Fluid Mech. 633, 425–435. (doi:10.1017/S0022112009008064)
11. Sankin GN, Yuan F, Zhong P 2010. Pulsating tandem microbubble for localized and directional single-cell membrane poration. Phys. Rev. Lett. 105, 078101 (doi:10.1103/PhysRevLett.105.078101) [PMC free article] [PubMed]
12. Gonzalez-Avila SR, Klaseboer E, Khoo BC, Ohl C-D 2011. Cavitation bubble dynamics in a liquid gap of variable height. J. Fluid Mech. 682, 241–260. (doi:10.1017/jfm.2011.212)
13. Yuan F, Sankin G, Zhong P 2011. Dynamics of tandem bubble interaction in a microfluidic channel. J. Acoust. Soc. Am. 130, 3339–3346. (doi:10.1121/1.3626134) [PubMed]
14. Prosperetti A. 2004. Bubbles. Phys. Fluids 16, 1852–1865. (doi:10.1063/1.1695308)
15. Ilinskii YA, Zabolotskaya EA, Hay TA, Hamilton MF 2012. Models of cylindrical bubble pulsation. J. Acoust. Soc. Am. 132, 1346–1357. (doi:10.1121/1.4730888) [PubMed]
16. Wang C, Rallabandi B, Hilgenfeldt S 2013. Frequency dependence and frequency control of microbubble streaming flows. Phys. Fluids 25, 022002 (doi:10.1063/1.4790803)
17. Viktorov IA. 1967. Rayleigh and Lamb waves: physical theory and applications. New York, NY: Plenum Press.
18. Rayleigh L. 1885. On waves propagated along the plane surface of an elastic solid. Proc. Lond. Math. Soc. s1–17, 4–11. (doi:10.1112/plms/s1-17.1.4)
19. Scholte JG. 1947. The range and existence of Rayleigh and Stoneley waves. Geophys. J. Int. 5, 120–126. (doi:10.1111/j.1365-246X.1947.tb00347.x)
20. Ewing WM, Jardetzky WS, Press F 1957. Elastic waves in layered media. New York, NY: McGraw-Hill.
21. Lamb H. 1917. On waves in an elastic plate. Proc. R. Soc. Lond. A 93, 114–128. (doi:10.1098/rspa.1917.0008)
22. Leighton TG, White PR, Schneider MF 1998. The detection and dimension of bubble entrainment and comminution. J. Acoust. Soc. Am. 103, 1825–1835. (doi:10.1121/1.421374)
23. Leighton TG, White PR, Morfey CL, Clarke JWL, Heald GJ, Dumbrell HA, Holland KR 2002. The effect of reverberation on the damping of bubbles. J. Acoust. Soc. Am. 112, 1366–1376. (doi:10.1121/1.1501895) [PubMed]
24. Leighton TG. 2011. The inertial terms in equations of motion for bubbles in tubular vessels or between plates. J. Acoust. Soc. Am. 130, 3333–3338. (doi:10.1121/1.3638132) [PubMed]
25. Achenbach JD. 1973. Wave propagation in elastic solids. Amsterdam, The Netherlands: North-Holland Publishing Company.
26. Piessens R. 2000. The Hankel transform. In The transforms and applications handbook (ed. Poularikas AD, editor. ), pp. 9.1–9.30. Boca Raton, FL: CRC Press LLC.
27. Landau LD, Lifshitz EM 1987. Fluid mechanics. Oxford, UK: Pergamon Press.
28. Zhu J, Popovics JS, Schubert F 2004. Leaky Rayleigh and Scholte waves at the fluid-solid interface subjected to transient point loading. J. Acoust. Soc. Am. 116, 2101–2110. (doi:10.1121/1.1791718)
29. Abramowitz M, Stegun IN 1965. Handbook of mathematical functions. New York, NY: Dover Publications.
30. Meyers M, Chawla K 2009. Mechanical behavior of materials. Cambridge, UK: Cambridge University Press.
31. Mekki-Berrada F, Thibault P, Marmottant P 2016. Acoustic pulsation of a microbubble confined between elastic walls. Phys. Fluids 28, 032004 (doi:10.1063/1.4942917)
32. Vrettos C. 2008. Green’s functions for vertical point load on an elastic half-space with depth-degrading stiffness. Eng. Anal. Bound. Elem. 32, 1037–1045. (doi:10.1016/j.enganabound.2007.10.017)
33. Churchill RV. 1973. Operational mathematics. New York, NY: McGraw-Hill.

Articles from Proceedings. Mathematical, Physical, and Engineering Sciences are provided here courtesy of The Royal Society