2.1 General formulation

In this section, we derive the mathematical formulation for the dynamics of two bubbles suspended in an unbounded fluid. We establish an axisymmetric coordinate system with the symmetric axis $Oz$ connecting the centres of the two bubbles, as shown in figure 1. In this coordinate system, the bubbles translate along $Oz$ . The distance of separation between the two centres is denoted by $L$ . The radial motions and deformations of the bubbles are described by the local polar coordinates ( $r_{1},\unicode[STIX]{x1D703}_{1}$ ) or ( $r_{2},\unicode[STIX]{x1D703}_{2}$ ). In the following equations, the subscript $i$ denotes the one bubble, and $j=3-i$ denotes the other one. Only the equations for bubble $i$ are written, and the equations for bubble $j$ are equivalent, with $i$ and $j$ exchanged.

Figure 1.

The coordinate system of two bubbles.

The velocity $\boldsymbol{u}$ and pressure $p$ of the external field are governed by the incompressible Navier–Stokes equations:

where $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D707}_{l}$ are the liquid density and viscosity, respectively. $\unicode[STIX]{x1D735}_{i}$ denotes the gradient, of which the origin is attached to the centroid of the bubble $i$ . $v_{i}$ is the translational velocity of bubble $i$ in the direction along $Oz$ . We neglect the flow field within the bubbles due to the negligible density and viscosity of gas. The internal pressure $p_{gi}$ is assumed to be spatially uniform and varies in time according to the polytropic process:

where $V_{bi}$ is the volume of bubble $i$ , and the subscript 0 indicates the initial value. $\unicode[STIX]{x1D6E4}$ is the polytropic index.

In the presence of the encapsulating membrane, the bubble interface $S_{i}$ is subjected to a no-velocity-jump condition, yielding continuities of the velocities in both the normal and tangential directions:

where the subscripts $m$ and $l$ represent the material points on the membrane and in the liquid, respectively. The continuity of velocity in the normal direction yields the kinematic boundary condition:

The dynamic boundary condition provides the balance of the pressure difference across the interface, the viscous stress of the fluid, and the membrane stress, i.e.

where $\unicode[STIX]{x1D644}_{i}$ and $\boldsymbol{n}_{i}$ are the unit tensor and unit normal vector with respect to the local polar coordinate, $\unicode[STIX]{x1D640}_{i}=(\unicode[STIX]{x1D735}_{i}\boldsymbol{u}+\unicode[STIX]{x1D735}_{i}\boldsymbol{u}^{\text{T}})/2$ the strain rate tensor, and $\boldsymbol{F}_{i}$ the membrane stress. The surface tension with coefficient $\unicode[STIX]{x1D6FE}$ is included here so that our model can be used for gas bubbles when membrane stress is set to be zero. For encapsulated bubbles, the surface tension is set to be zero. At infinity, we assume that

We decompose velocity and pressure into the potential (with subscript $p$ ) and the viscous correction (with subscript $v$ ) parts:

We will discuss the solution procedures for these two parts in detail in the later sections.

2.2 Geometry of the bubble

We assume that the dominant motion is in the radial direction with perturbations due to translational motion and deformations. The perturbation to a spherical bubble can be expanded in terms of the spherical harmonics. Since we consider the axisymmetric model for simplicity in this paper, the spherical harmonics reduces to the Legendre polynomials, as the azimuthal mode is constant. For a gas bubble, the normal stress jump across the interface due to the surface tension is a function of the curvature, computed using the normal component of the interfacial displacement. By contrast, for an encapsulated bubble, the in-plane stress occurs as a resistance to the tangential displacement of the membrane. Therefore, in the present study, the tangential displacement must be computed, and is written using the associated Legendre polynomial $P_{k}^{1}(\cos \unicode[STIX]{x1D703}_{i})(=\text{d}P_{k}(\cos \unicode[STIX]{x1D703}_{i})/\text{d}\unicode[STIX]{x1D703}_{i})$ . The position of the material point on the membrane is therefrom expressed in radial and zenith directions:

The series of Legendre polynomials $P_{k}(\cos \unicode[STIX]{x1D703}_{i})$ represents the perturbation in the radial direction with amplitudes of shape modes $a_{ki}$ . Notice $a_{1i}$ is the translational distance with the bubble shape unchanged and is related to the translational velocity $v_{i}$ , so that the summation in (2.9) is from $k=2$ . The series of associated Legendre polynomials $P_{k}^{1}(\cos \unicode[STIX]{x1D703}_{i})$ is the perturbation in the tangential direction due to deformation. For a special case, in which the volume and the shape are unchanged, as illustrated in figure 2, $b_{1i}P_{1}^{1}(\cos \unicode[STIX]{x1D703}_{i})$ denotes the extension or compression along the surface, inducing strain and thus stress in the membrane.

Figure 2.

The movements of material points on the membrane at $b_{1i}P_{1}^{1}(\cos \unicode[STIX]{x1D703}_{i})$ .

We assume an infinitesimal element $AB$ in the reference configuration is unstressed. As sketched in figure 3, $AB$ transforms to $A^{\prime }B^{\prime }$ after some time. According to the above definition of surface position, the coordinates of the material points $A$ , $B$ , $A^{\prime }$ and $B^{\prime }$ are

respectively. The arclength between $A$ and $B$ is

and the arclength between $A^{\prime }$ and $B^{\prime }$ is

Figure 3.

Sketch of movements of material points.

The two principal strains $\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D703}i}$ and $\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D719}i}$ are defined by the extension ratio in the arclength direction and the change ratio of the distance from the axis of symmetry, respectively:

On the other hand, the shape of the bubble is independent of the mode $b_{ki}P_{k}^{1}(\cos \unicode[STIX]{x1D703}_{i})$ . Therefore, we can define the surface function as

Other geometric quantities are subsequently obtained by the surface function (Liu et al. Reference Liu, Sugiyama, Takagi and Matsumoto2012), including the unit normal vector $\boldsymbol{n}_{i}$ pointing into the fluid:

the tangential vector $\boldsymbol{t}_{i}$ :

and the two principal curvatures in the $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D719}$ directions:

2.3 Potential part

In this section, we will derive the solution to the potential part of the velocity $\boldsymbol{u}_{pi}$ and pressure $p_{pi}$ . The velocity potential is expressed as the superposition of two Legendre polynomials (Takahira et al. Reference Takahira, Akamatsu and Fujikawa1991):

where $q_{ki}$ and $q_{kj}$ are the coefficients to be determined, in which the first subscript $k$ denotes the order of mode. Using a relation of the Legendre polynomials (Hobson Reference Hobson1931):

(2.21) becomes:

in which the second term in the right-hand side represents the influences of the nearby bubble. Its effects diminish as the mode number increases (as $l$ increases). For simplicity, we consider only the effects of radial and translational motions from the nearby bubble by selecting $l=0$ and 1, and rewrite (2.23) as

It is noted that among the coefficients related to the nearby bubble $j$ , only $q_{0j}$ and $q_{1j}$ are included. The velocity of the potential part is obtained subsequently by

The pressure of the potential part is obtained by integrating the Navier–Stokes equation:

Substituting the relationship between velocity and its potential (2.25) and the surface function (2.16) into the kinematic boundary condition (2.5), and making use of the orthogonal relations of the Legendre polynomials (appendix A), we obtain the coefficients for each mode:

where $v_{ri}=v_{i}-q_{0j}/L^{2}$ . Substituting these coefficients into the velocity potential (2.24), the velocity (2.25) and pressure (2.26), we obtain the formula for the potential part of the flow field.

2.4 Viscous correction

The viscous correction in (2.8) is solved by following Prosperetti (Reference Prosperetti1977). We assume the boundary layers of the two bubbles do not overlap, which will be justified in § 2.6, so that the viscous correction can be obtained separately for each bubble. Under this assumption, the separation distance of the two bubbles should be much larger than the thicknesses of the boundary layers. We estimate the thickness of the boundary layer by (Hilgenfeldt, Lohse & Brenner Reference Hilgenfeldt, Lohse and Brenner1996):

where $\unicode[STIX]{x1D714}_{d}$ is the driving frequency.

Here we briefly review the theory proposed by Prosperetti (Reference Prosperetti1977). The viscous correction is obtained by decomposing vorticity into poloidal and toroidal fields. Only the toroidal field $T_{ki}(r_{i},t)$ is relevant since we consider an axisymmetric model, which yields:

The velocity and pressure are written in terms of the toroidal field as:

where

The partial differential equation (2.31) will be solved numerically based on the boundary conditions given below. The velocity $\boldsymbol{u}_{vi}$ and pressure $p_{vi}$ are subsequently obtained using (2.32) and (2.33).

2.5 Stress balance on the interfaces

We write the dynamic boundary condition (2.6) in the tangential and normal directions, respectively:

The fluid pressure $p$ is composed of the potential part (2.26) and the viscous correction part (2.33). The pressure inside the bubble $p_{gi}$ obeys the polytropic law (2.3). The surface tension for gas bubbles is derived by taking the divergence of the unit normal vector (2.17):

We model the membrane as a viscoelastic material, employing the neo-Hookean law (Mooney Reference Mooney1940) for the in-plane stress, a linear law for the membrane viscosity, and the Love law (Love Reference Love1888) for the bending moments. The membrane stresses $F_{ti}$ and $F_{ni}$ are given by

where $G_{si}$ is the surface elasticity modulus, $G_{bi}$ the bending modulus, $\unicode[STIX]{x1D707}_{si}$ the membrane viscosity and $\unicode[STIX]{x1D710}_{i}$ the Poisson ratio. The detailed derivation of the membrane stresses is given in appendix B or Liu et al. (Reference Liu, Sugiyama, Takagi and Matsumoto2012). The difference from the single bubble case is the stresses from the mode of $P_{1}(\cos \unicode[STIX]{x1D703}_{i})$ and $P_{1}^{1}(\cos \unicode[STIX]{x1D703}_{i})$ due to the tangential movements of the material points on the undeformed interfaces.

The rate of strain can be expressed in the component form:

where

Thus far, all the terms in (2.35) and (2.36) have been derived.

2.6 The final derived equations

Expanding the stress boundary condition in the tangential direction (2.35) with respect to the associated Legendre polynomials, and applying the orthogonality of the Legendre polynomials, we obtain the following equations:

Specifically, when $k=1$ , we have

Equation (2.45) or (2.46) will be used as one of the boundary conditions of the toroidal equation (2.31). Another boundary condition for (2.31) is the vanishing vorticity beyond the boundary layer: $T_{ki}\rightarrow 0$ as $r_{i}\rightarrow \infty$ . The initial condition for (2.31) is the static flow field: $T_{ki}(t=0)=0$ . Using these boundary and initial conditions, we solve the toroidal equation by means of the finite difference method. The integral related to $T_{ki}(r_{i},t)$ is calculated by composite numerical quadratures.

The expansion of the stress boundary condition (2.36) in terms of Legendre polynomials provides additional equations at the specific order of mode $k$ . The zeroth order ( $k=0$ ) represents the radial oscillation:

The left-hand side is the famous Rayleigh–Plesset equation plus the effects of the radial oscillation of the nearby bubble and the translational motions. The right-hand side of this equation is due to the membrane stress. Equation (2.47) describes the nonlinear pulsation in the radial mode.

We now discuss the rest of normal stress condition (2.36). The prefactor before $P_{1}(\cos \unicode[STIX]{x1D703}_{i})$ in (2.36) represents the translational motion. We express the translation equation in the form of inertia of added mass and forces:

where

and

The right-hand side of (2.48) is the forces exerting on an encapsulated bubble, projecting onto the translational direction. They are written as the sum of several parts (2.49)–(2.54), and their respective influences will be investigated in the next section. Equation (2.49) gives the secondary Bjerknes force, which is generated by the radial pulsation of the nearby bubble, which is indicated in $q_{0j}$ of (2.49). $F_{h1i}$ , $F_{h2i}$ and $F_{h3i}$ are higher-order perturbations. $F_{h1i}$ has a form similar to that of the classical Bjerknes force, representing the effect of the translational motion of the nearby bubble; $F_{h2i}$ comes from the coupling of its translational motion and the radial and translational motions of the nearby bubble; $F_{h3i}$ is due to the coupling of its translational motion and the second-order shape oscillation of the nearby bubble. $F_{mi}$ is the stress along the undeformed surface, denoted by the amplitude of $b_{1i}$ . $F_{\unicode[STIX]{x1D707}i}$ is the viscous drag. It includes the viscous friction on the bubble surface and the vorticity distribution due to the presence of the encapsulating membrane. The first term on the right-hand side of (2.53) is the Levich drag (Levich Reference Levich1949), which is the viscous drag of a spherical gas bubble with a free-slip surface in a potential flow. The subsequent terms are the contributions of viscous dissipation generated by the membrane. In this paper, we approximate the flow field as a potential flow with viscous correction. This approximation is valid when the two boundary layers around the bubbles do not overlap. For water, typical driving frequency $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 1~\text{MHz}$ , typical bubble radius $R_{0i}=2\sim 10~\unicode[STIX]{x03BC}\text{m}$ and typical distance between bubbles $L_{0}=3R_{01}$ , the thickness of the boundary layer calculated by (2.30) is less than $1/10$ th of the bubble radius $R_{0i}$ . This indicates that the vorticity generated by the two bubbles hardly interacts, and the Womersley number $(Wo=R_{i}/\unicode[STIX]{x1D712}_{i})\sim O(10)$ is high. Magnaudet & Eames (Reference Magnaudet and Eames2000) reported that when the boundary layer thickness is much smaller than the radius, viscous drag force can be approximated by the Levich drag even if the Reynolds number related to translation is small (also see Magnaudet & Legendre (Reference Magnaudet and Legendre1998)). Our main objective is to investigate the behaviour of two bubbles under ultrasound irradiation, for which the translational motion is oscillatory, and the short-time response of an accelerated flow dominates. In this situation, vorticity generated at the bubble surface during oscillation changes its sign within a short period, so that the boundary layer is confined.

The history force, generated from the unsteady diffusion of vorticity, is calculated by temporally solving the spatial distribution of the toroidal field $T_{1i}(r_{i},t)$ , given by (2.31), with (2.46) as a boundary condition. Combined with (2.46), the last two terms in (2.48) are

This expression indicates that the drag force comes from the toroidal field related to vorticity transport. It is the force excluding the added mass inertia and the effect of the second bubble. The last term in (2.56) is related to the fore–aft asymmetry of the vorticity distribution, which comes from the viscous correction part of the velocity field. This term vanishes when the integrand is evaluated at the bubble surface. Therefore, the dominant viscous drag is

This expression of drag force is suitable for both gas bubbles and encapsulated bubbles, depending on the boundary condition (2.46). Equation (2.46) is the tangential balance of membrane stress with the viscous friction due to the fluid. For a gas bubble, one simply drops the membrane stress.

The toroidal field equation (2.31), which is solved for the viscous correction, is the linearized Navier–Stokes equation for the vorticity transport (Prosperetti Reference Prosperetti1977). It is equivalent to solving the unsteady Stokes equation, whose solution is obtained analytically in appendix D. For the special case of an isovolumetric sphere with a no-slip or free-slip interface in an oscillatory motion, the analytical solution is well-known. Following the derivation in Landau & Lifshitz (Reference Landau and Lifshitz1987), we obtain the analytical solution for the encapsulated bubble as an extension of the classical solution (appendix C). By substituting the solution (D 3) into $T_{1i}$ in (2.57), we obtain a consistent expression for the drag force in (C 10). Using the derivation in appendix C, $F_{Di}$ becomes the Levich force for a gas bubble with a free-slip surface:

For a rigid sphere with a no-slip interface, $F_{Di}$ includes the history force and the quasisteady drag force:

For an encapsulated bubble with the material points moving along the surface, $F_{Di}$ is obtained numerically and its value lies in the intermediate range between (2.58) and (2.59).

The dynamic equation of shape oscillation at the $k$ th order is given by:

where

The first three terms of (2.60) account for the shape oscillation of a single bubble. The remaining terms come from the influences of radial oscillation and translational motion of the nearby bubble, the coupling effects with its own translational motion, the influences of both the higher- and lower-order shape modes $a_{k+1,i}$ and $a_{k-1,i}$ , and the viscous effects. For $k=2$ , $a_{k-1,i}$ is related to the translation. Here we describe the translation in terms of the velocity $v_{i}$ (or $v_{ri}$ ), so that we have a special expression for $G_{k}$ when $k=2$ (2.61). The right-hand side is the membrane stress at the order of $P_{k}(\cos \unicode[STIX]{x1D703}_{i})$ .

The continuity in tangential velocity is given by (2.4)

which can be written explicitly as

when $k=1$ it becomes

To summarize, the unknowns for bubble $i:R_{i},v_{i},b_{1i},a_{ki},b_{ki}~(k\geqslant 2)$ are obtained by solving (2.47), (2.48), (2.65), (2.60) and (2.64) by using the fourth-order Runge–Kutta method.

3.1 Membrane effects

To compare with the experimental results by Garbin et al. (Reference Garbin, Dollet, Overvelde, Cojoc, Di Fabrizio, Van Wijngaarden, Prosperetti, De Jong, Lohse and Versluis2009), we consider two bubbles of initial radius $2~\unicode[STIX]{x03BC}\text{m}$ and separation distance $8~\unicode[STIX]{x03BC}\text{m}$ . An eight-cycle 2.25 MHz ultrasound pulse is applied with an amplitude adjusted to fit the radial variation in the experiment (figure 4). The relative translation is calculated and shows in figure 5 that the two bubbles with the same size move towards each other. Our calculated translations for encapsulated bubbles with elastic membranes (the solid line in figure 5) damp more slowly than those observed in the experiment by Garbin et al. (Reference Garbin, Dollet, Overvelde, Cojoc, Di Fabrizio, Van Wijngaarden, Prosperetti, De Jong, Lohse and Versluis2009). The reason is that the coated bubbles used in Garbin et al. (Reference Garbin, Dollet, Overvelde, Cojoc, Di Fabrizio, Van Wijngaarden, Prosperetti, De Jong, Lohse and Versluis2009) are no-slip at the bubble surface, while our model allows the material points to move along the surface according to a viscoelastic constitutive law. As explained in § 2.6, the drag of an encapsulated bubble with viscoelastic membrane lies between the values of no-slip and free-slip interfaces. If we increase the elastic modulus, the membrane becomes hard and the behaviour of the encapsulated bubble should approach that of a no-slip bubble. As shown by the dotted curve in figure 5, the variation of the distance between the bubbles with large elastic modulus ( $G_{si}\times 10^{3}$ ) approaches that for no-slip bubbles (circle symbols). If we decrease the elastic modulus, the membrane becomes soft and the drag reduces, leading to a larger translational distance, as shown by the line with square symbols in figure 5 for $G_{si}\times 10^{-3}$ . When we drop the membrane stress and replace it with surface tension, we obtain the result for gas bubbles with a free-slip interface (the dashed line in figure 5), whose damping in translation is much smaller than that of encapsulated bubbles. This confirms that the no-velocity-jump condition for an elastic interface increases the drag significantly. The no-velocity-jump condition induces vorticity generation and thus contributes to the history force in the total drag force. The magnitude of history force depends on the mobility of the material points at the interface, revealed by the tangential displacement $b_{ki}$ . Figure 6 shows the time development of $b_{1i}$ , denoting the extension or compression along the undeformed interface. The material points on the free-slip interface have a much larger mobility, related to a smaller history force in drag than that of encapsulated bubbles. The comparison among encapsulated bubbles with an elastic membrane (with a normal $G_{si}$ ), a hard membrane (with $G_{si}\times 10^{3}$ ) and a soft membrane (with $G_{si}\times 10^{-3}$ ) reveals that the membrane stress suppresses the mobility of the material points, and this suppression effect is enhanced as the modulus increases. The membrane with a sufficiently large elastic modulus approaches the condition of a no-slip interface. The substantial difference in the tangential movements of material points between encapsulated bubbles and gas bubbles lies in the variation of $b_{1i}$ after the ultrasound driving is removed at $t=3.56~\unicode[STIX]{x03BC}\text{s}$ . For encapsulated bubbles, the stress in the membrane acts as a restoring force to move the material points back to their equilibrium positions. In these cases, $b_{1i}$ returns to zero after $t=3.56~\unicode[STIX]{x03BC}\text{s}$ . By comparison, the surface tension in the interface of a gas bubble does not have the restoring effect. Therefore, $b_{1i}$ for the gas bubble remains unchanged and will not return to zero after the ultrasound driving stops. For the deformation, we take the shape mode at the lowest order $k=2$ as an example and investigate the time developments of $a_{2i}$ and $b_{2i}$ , as shown in figures 7 and 8, respectively. The amplitudes of $a_{2i}$ and $b_{2i}$ for the encapsulated bubble with an elastic membrane are smaller than those with free-slip interfaces. The membrane has a stabilizing effect on the interface, making the deformation smaller than the gas bubble. The stabilizing effect becomes more appreciable as the elastic modulus increases. For an extremely large elastic modulus, the model gives a rigid sphere that would not deform. If one sets the membrane parameters to zero and assigns the surface tension, the model is for gas bubbles. By choosing appropriate interface parameters, the model can be used for free-slip, elastic and no-slip interfaces.

Figure 6.

Time evolutions of $b_{1i}$ for encapsulated bubbles with elastic membranes (solid line), gas bubbles with free-slip interfaces (dashed line), bubbles with hard membranes (line with circle symbols) and bubbles with soft membranes (line with square symbols). Inset: detailed comparison between encapsulated bubbles with elastic and hard membranes.

Figure 7.

Time evolutions of $a_{2i}$ for encapsulated bubbles with elastic membranes (solid line), gas bubbles with free-slip interfaces (dashed line), bubbles with hard membranes (line with circle symbols) and bubbles with soft membranes (line with square symbols).

Figure 8.

Time evolutions of $b_{2i}$ for encapsulated bubbles with elastic membranes (solid line), gas bubbles with free-slip interfaces (dashed line), bubbles with hard membranes (line with circle symbols) and bubbles with soft membranes (line with square symbols).

3.2 Validation of the potential part

Under the potential theory, the case of two equal-size bubbles is equivalent to that of a single bubble near a wall which is located at the midpoint between the centres of two bubbles. As a validation of the potential part in our model, we compare our results with the experimental observation by Lauterborn & Bolle (Reference Lauterborn and Bolle1975). The initial radius is $R_{01}=R_{02}=2.6\times 10^{-3}~\text{m}$ . The initial distance is $3R_{01}$ , equivalent to the distance from the wall being $1.5R_{01}$ . This distance is much larger than the boundary layer thickness around the bubble, which is calculated by (2.30) to be $\unicode[STIX]{x1D712}/R_{01}<0.1\,\%$ . In addition, the experimental results by Lauterborn & Bolle (Reference Lauterborn and Bolle1975) agree well with the numerical results by Plesset & Chapman (Reference Plesset and Chapman1971), who solved the bubble–wall interaction based on potential flow using the boundary integral method. Plesset & Chapman (Reference Plesset and Chapman1971) assumed that the normal derivative vanishes at the solid wall, implying that their bubble–wall system is equivalent to two equal-size bubbles. The agreement between Lauterborn & Bolle (Reference Lauterborn and Bolle1975) and Plesset & Chapman (Reference Plesset and Chapman1971) indicates the viscous effect is unimportant for the present radius and bubble–wall distance, and it is the added mass inertia and the force due to the pressure obeying the Bernoulli equation that dominate.

The pressure outside the bubbles is set to be $p_{\infty }=1\times 10^{5}~\text{Pa}$ , and the pressure inside the bubbles is $p_{g01}=p_{g02}=2337~\text{Pa}$ initially. The solid lines in figure 9 show the bubble shapes when $\bar{t}=t/[R_{01}(\unicode[STIX]{x1D70C}/\unicode[STIX]{x0394}p)^{1/2}]=0$ , 0.725, 0.825, 0.961, 0.991, 1.016, from larger to smaller shapes. The dotted symbols denote the experimental results at the same times. We can see the results agree well initially. The bubble is elongated in the direction perpendicular to the wall, and moves towards the wall. The results differ slightly at the last two times. Our model is based on a linear approximation for the higher-order shape mode, hence less accurate when the deformation is not small. However, the tendency towards jet formation at the backside of the bubble during collapse is predicted as the third-order shape mode emerges. In general, the agreement with the experimental result is satisfactory.

Figure 10 shows the behaviour of two encapsulated bubbles of equal size. We choose the same parameters as those in figure 9, except that the surface tension is neglected and the membrane model is used instead. We calculate the natural frequency of the radial mode using (Liu et al. Reference Liu, Sugiyama, Takagi and Matsumoto2012)

Coincidentally, the natural frequency of this encapsulated bubble is $2\unicode[STIX]{x03C0}\times 1258~\text{Hz}$ , close to that of the gas bubble ( $2\unicode[STIX]{x03C0}\times 1255~\text{Hz}$ ). Comparing the shapes at $\bar{t}=1.016$ , the deformation of the encapsulated bubble is smaller than that of the gas bubble, implying the interface of the encapsulated bubble is more stable.

Figure 9.

Two gas bubbles with equal size. $R_{01}=R_{02}=2.6\times 10^{-3}~\text{m}$ . The initial distance between the bubbles is $3R_{01}$ ; $p_{\infty }=1\times 10^{5}~\text{Pa}$ ; $p_{g01}=p_{g02}=2337~\text{Pa}$ . The solid lines denote the bubble shapes at $\bar{t}=0$ , 0.725, 0.825, 0.961, 0.991, 1.016, from larger to smaller shapes. The dotted symbols are the experimental results by Lauterborn & Bolle (Reference Lauterborn and Bolle1975).

Figure 10.

Two encapsulated bubbles with equal size. $R_{01}=R_{02}=2.6\times 10^{-3}~\text{m}$ . The initial distance between the bubbles is $3R_{01}$ ; $p_{\infty }=1\times 10^{5}~\text{Pa}$ ; $p_{g01}=p_{g02}=2337~\text{Pa}$ . The solid lines denote the bubble shapes at $\bar{t}=0$ , 0.725, 0.825, 0.961, 0.991, 1.016, from larger to smaller shapes. The dotted lines with symbols are the shapes of a gas bubble reproduced from figure 9.

Figure 11.

The final distances after 20 driving cycles at different driving frequencies with driving amplitude $\unicode[STIX]{x1D700}=0.1$ . $R_{01}=10~\unicode[STIX]{x03BC}\text{m}$ , $R_{02}=5~\unicode[STIX]{x03BC}\text{m}$ . (a) Gas bubbles, (b) encapsulated bubbles.

Figure 12.

The variations of radius at different driving frequencies for the gas bubbles at $\unicode[STIX]{x1D700}=0.1$ . $R_{01}=10~\unicode[STIX]{x03BC}\text{m}$ , $R_{02}=5~\unicode[STIX]{x03BC}\text{m}$ . (a) $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 0.1~\text{MHz}$ ; (b) $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 0.3~\text{MHz}$ ; (c) $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 0.5~\text{MHz}$ ; (d) $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 0.7~\text{MHz}$ ; (e) $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 0.9~\text{MHz}$ .

3.3 Comparison with Bjerknes’ theory

According to Bjerknes’ theory, the two bubbles will move away from each other when the driving frequency is between the two natural frequencies $\unicode[STIX]{x1D714}_{01}<\unicode[STIX]{x1D714}_{d}<\unicode[STIX]{x1D714}_{02}$ ; and when $\unicode[STIX]{x1D714}_{d}<\unicode[STIX]{x1D714}_{01}$ or $\unicode[STIX]{x1D714}_{d}>\unicode[STIX]{x1D714}_{02}$ , the two bubbles will move towards each other. We choose two bubbles with $R_{01}=10~\unicode[STIX]{x03BC}\text{m}$ and $R_{02}=5~\unicode[STIX]{x03BC}\text{m}$ and set their initial distance as $3R_{01}$ . A driving pressure characterized by a dimensionless amplitude $\unicode[STIX]{x1D700}$ and a driving frequency $\unicode[STIX]{x1D714}_{d}$ is applied:

where the ambient pressure $p_{0}$ is set as $1\times 10^{5}~\text{Pa}$ . We plot the final distances between the two bubbles after 20 driving cycles by sweeping the driving frequency in figure 11. Here the driving amplitude is $\unicode[STIX]{x1D700}=0.1$ . Figure 11(a) is the result for gas bubbles, and figure 11(b) for encapsulated bubbles. It is noted that the zeroth-order natural frequencies of the gas bubbles are $\unicode[STIX]{x1D714}_{01}=2\unicode[STIX]{x03C0}\times 0.34~\text{MHz}$ and $\unicode[STIX]{x1D714}_{02}=2\unicode[STIX]{x03C0}\times 0.72~\text{MHz}$ , respectively. From figure 11(a), we can see the two bubbles are attracted to each other when the driving frequency is less than both natural frequencies. When the driving frequency approaches the smaller natural frequency, i.e. $\unicode[STIX]{x1D714}_{d}=\unicode[STIX]{x1D714}_{01}$ , the distance reaches a minimum. After this minimum, the relative translational motion becomes repulsive. As the driving frequency increases to the vicinity of the larger natural frequency, i.e. $\unicode[STIX]{x1D714}_{d}=\unicode[STIX]{x1D714}_{02}$ , the two bubbles repel each other to the maximal distance. After this maximum, the bubbles return to being attractive. This picture agrees well with Bjerknes’ theory. Figure 11(a) shows that the resonance in radial oscillation leads to an extreme translation. At the same time, the resonance at subharmonics causes a secondary extreme translation. In figure 11(a), there is a local minimum near $2\unicode[STIX]{x03C0}\times 0.17~\text{MHz}$ , which is half of the zeroth-order natural frequency $\unicode[STIX]{x1D714}_{01}$ . Half of another natural frequency ( $\unicode[STIX]{x1D714}_{02}$ ) is $2\unicode[STIX]{x03C0}\times 0.36~\text{MHz}$ , close to $\unicode[STIX]{x1D714}_{01}=2\unicode[STIX]{x03C0}\times 0.34~\text{MHz}$ , which explains why the minimum near $2\unicode[STIX]{x03C0}\times 0.34~\text{MHz}$ has a broad spread to around $2\unicode[STIX]{x03C0}\times 0.36~\text{MHz}$ .

Figure 11(b) shows the results for encapsulated bubbles. Under the same driving amplitude, the translational distances of encapsulated bubbles are much smaller than those of gas bubbles. Our calculated zeroth-order natural frequencies of the two encapsulated bubbles are $\unicode[STIX]{x1D714}_{01}=2\unicode[STIX]{x03C0}\times 0.51~\text{MHz}$ and $\unicode[STIX]{x1D714}_{02}=2\unicode[STIX]{x03C0}\times 1.28~\text{MHz}$ , respectively. The minimal distance occurs near $\unicode[STIX]{x1D714}_{01}=2\unicode[STIX]{x03C0}\times 0.51~\text{MHz}$ and the maximum occurs near $\unicode[STIX]{x1D714}_{02}=2\unicode[STIX]{x03C0}\times 1.28~\text{MHz}$ . The behaviour in translation follows the prediction by the Bjerknes theory. Different from the gas bubble case, no local minimum appears near subharmonic resonance for encapsulated bubbles, which is attributed to the damping effects due to the encapsulating membrane, making the behaviour of encapsulated bubbles more stable.

Bjerknes’ theory also states that the bubbles move towards each other when oscillating in phase ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}<\unicode[STIX]{x03C0}/2$ or $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}>3\unicode[STIX]{x03C0}/2$ ); otherwise, they repel each other when oscillating out of phase ( $\unicode[STIX]{x03C0}/2<\unicode[STIX]{x0394}\unicode[STIX]{x1D719}<3\unicode[STIX]{x03C0}/2$ ), where $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ is the difference in phase angle. To examine the effects of the phase difference, we choose a set of driving frequencies: smaller than $\unicode[STIX]{x1D714}_{01}$ , near $\unicode[STIX]{x1D714}_{01}$ , between $\unicode[STIX]{x1D714}_{01}$ and $\unicode[STIX]{x1D714}_{02}$ , near $\unicode[STIX]{x1D714}_{02}$ , and larger than $\unicode[STIX]{x1D714}_{02}$ , that is, $\unicode[STIX]{x1D714}_{d}/2\unicode[STIX]{x03C0}=0.1$ , 0.3, 0.5, 0.7, 0.9 MHz for gas bubbles and $\unicode[STIX]{x1D714}_{d}/2\unicode[STIX]{x03C0}=0.3$ , 0.5, 0.8, 1.3, 1.5 MHz for encapsulated bubbles. The time evolutions of radius are plotted in figures 12 and 13, respectively. We can see the in-phase oscillations in figures 12(a), 12(b), 13(a) and 13(b), corresponding to the attractive motion we have observed in figure 11. For the driving frequencies greater than both natural frequencies, $\unicode[STIX]{x1D714}_{d}/2\unicode[STIX]{x03C0}=0.9~\text{MHz}$ for gas bubbles and $\unicode[STIX]{x1D714}_{d}/2\unicode[STIX]{x03C0}=1.5~\text{MHz}$ for encapsulated bubbles, the bubbles exhibit slight translation even when the two bubbles oscillate out of phase. This is perhaps due to the fact that the bubbles do not have enough time to move from their initial positions and keep a stable separation under such high driving frequencies, relating to ‘bubble grape’, which is supposed to form at high frequency (Harkin et al. Reference Harkin, Kaper and Nadim2001). For other cases (figures 12 c, 12 d, 13 c and 13 d), the two bubbles oscillate out of phase and consequently move away from each other. These results show that the phase difference comes from the nonlinear oscillation subjected to different driving frequencies, which are responsible for the attraction or repulsion of the two bubbles.

Figure 13.

The variations of radius at different driving frequencies for the encapsulated bubbles at $\unicode[STIX]{x1D700}=0.1$ . $R_{01}=10~\unicode[STIX]{x03BC}\text{m}$ , $R_{02}=5~\unicode[STIX]{x03BC}\text{m}$ . (a) $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 0.3~\text{MHz}$ ; (b) $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 0.5~\text{MHz}$ ; (c) $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 0.8~\text{MHz}$ ; (d) $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 1.3~\text{MHz}$ ; (e) $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 1.5~\text{MHz}$ .

3.4 Contributions of different forces

In this section, we investigate the roles of the different forces in (2.49)–(2.54). We choose the driving frequency $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 0.3~\text{MHz}$ for gas bubbles and $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 0.5~\text{MHz}$ for encapsulated bubbles. Under these driving frequencies, the larger bubble experiences radial resonance and the two bubbles exhibit maximal attractive motion. The forces are averaged over a driving cycle:

so are the distances and the velocities. Figures 14(a) and 14(b) show the time evolutions of the averaged forces exerting on bubble 1 (the larger bubble) for gas bubbles and encapsulated bubbles, respectively; figures 14(c) and 14(d) are for bubble 2 (the smaller bubble). Figures 14(e) and 14(f) show the respective variations of distances for gas and encapsulated bubbles, and figures 14(g) and 14(h) show the variations of velocities.

Figure 14.

The time evolutions of forces acting on the bubbles, the distance between the two bubbles and the velocity of each bubble. (a,c,e,g) Gas bubbles at $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 0.3~\text{MHz}$ . (b,d,f,h) Encapsulated bubbles at $\unicode[STIX]{x1D714}_{d}=2\unicode[STIX]{x03C0}\times 0.5~\text{MHz}$ . (a,b) Forces on the larger bubbles (bubble 1). (c,d) Forces on the smaller bubbles (bubble 2). (e,f) Distance variations. (g,h) Variation of velocity of each bubble.

As shown in figures 14(a) and 14(c), and figures 14(b) and 14(d), the time-averaged Bjerknes forces follow the envelopes of radial variation in figures 12(b) and 13(b), respectively. This is natural since the Bjerknes force is due to the radial oscillation of the nearby bubble. In general, $F_{B1}$ is positive, and $F_{B2}$ is negative. Referred to the coordinates we established in figure 1, the two bubbles move towards each other and lead to a decreasing distance (figures 14 e and 14 f). The forces with opposite sign of the Bjerknes force behave as a drag. For gas bubbles, the drag mainly comes from the viscous forces $F_{\unicode[STIX]{x1D707}i}$ , in which the quasisteady drag force $12\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}_{l}R_{i}v_{ri}$ dominates. For encapsulated bubbles, on the other hand, the drag is mainly from the membrane forces $F_{mi}$ , contributing to the history force through the tangential dynamic boundary condition (2.46).

$F_{h1i}$ for bubble 1 is a drag, since its sign is opposite to the Bjerknes force, resulting from the translation of bubble 2 (the smaller bubble), which has a larger (absolute) velocity than the larger bubble (bubble 1) (figures 14 g and 14 h). It is understandable, since the magnitudes of forces of the two bubbles are comparable (figures 14 a and 14 c, figures 14 b and 14 d), so that the smaller bubble experiences a larger acceleration and velocity. This larger velocity of bubble 2 contributes to the force $F_{h1i}$ on bubble 1, which is related to the Bjerknes-like force due to the translation of the nearby bubble. The other higher-order forces $F_{h2i}$ and $F_{h3i}$ are small, suggesting that the coupling effect of the translation and the radial oscillation of the nearby bubble is not important, and the bubbles do not exhibit appreciable deformation in the present cases.

The presence of the membrane does not necessarily damp radial oscillations. The amplitude of radial oscillation depends on its frequency. At resonance, the oscillatory amplitude of the encapsulated bubble is comparable to that of the gas bubble, as shown in figures 12(b) and 13(b). Determined by the radial oscillation of the nearby bubble, the magnitudes of the Bjerknes forces of encapsulated bubbles are thereby of the same order as those of gas bubbles, as revealed by comparing figures 14(a,c) and figures 14(b,d). However, the no-velocity-jump condition on membranes contributes significantly to the viscous drag, which has been explained in § 2.6, as it increases the viscous dissipation. The expression for the viscous drag (2.57) includes the history force, which was shown to be important for encapsulated bubbles in the previous sections and by Garbin et al. (Reference Garbin, Dollet, Overvelde, Cojoc, Di Fabrizio, Van Wijngaarden, Prosperetti, De Jong, Lohse and Versluis2009). Therefore, the translational distances of encapsulated bubbles are smaller than those of gas bubbles even when they experience comparative radial oscillations. In addition, due to the smaller translations of encapsulated bubbles compared with gas bubbles, the boundary layers, which are more important for encapsulated bubbles with no-velocity-jump condition than for gas bubbles with free-slip condition, do not overlap. This ensures that our model of assumption of potential flow with viscous correction remains valid.