Birmingham Microbubble dynamics in a viscous compressible liquid subject to ultrasound

When a microbubble is subject to ultrasound, non-spherical oscillation or surface modes can be generated after many acoustic cycles. This phenomenon has wide applications, including ultrasonic cleaning, sonochemistry, and biomedical ultrasonics. Yet, the nonlinear development of the bubble shape modes over dozens of cycles is not well understood. Here, we describe a grid-free and robust model to simulate the phenomenon. A viscous pressure correction is introduced to compensate the non-zero tangential stress at the free surface in the potential ﬂow model, based on conservation of energy. Consequently, the phenomenon is modeled using the boundary integral method, in which the compressible and viscous effects are incorporated into the model through the boundary conditions. The computations have been carried out for axisymmetric cases; however, the numerical model can be extended for three-dimensional cases in a straightforward manner. The numerical results are shown to be in good agreement for many cycles with some independent viscous and compressible theories for axisymmetric bubbles and experiments for microbubbles undergoing shape oscillation subject to ultrasound. The development of the shape oscillation of a bubble after a dozen cycles, the formation of a reentry jet and its penetration through the bubble, and the topological transformation of the bubble are simulated and analyzed in terms of the amplitude and frequency of the ultrasound. The computations and physical analysis are carried out for the development of shape modes due to a resonant volume oscillation, strong pressure wave, or the matching of the acoustic wave frequency with the shape mode frequency.

It is observed in experiments that bubbles may be activated into repeated stable shape oscillations in an acoustic field (Asaki and Marston, 1995;Versluis et al., 2010).The experimental results are limited by the resolution in space and time, especially for microbubbles whose size and period are at O(10 À6 ) m at O(10 À6 ) s. Theoretical studies were carried on shape oscillations of bubbles at small amplitude using perturbation methods via spherical harmonics, predicting the natural frequency of shape modes and the stability threshold (Plesset, 1954;Prosperetti, 1977;Shaw, 2009Shaw, , 2017;;Doinikov, 2004;Gu edra et al., 2017;Gu edra and Inserra, 2018).We aim to implement a numerical model to simulate and analyze the development of shape modes of bubbles at large amplitude.
Viscous effects are important for microbubbles, since the associated Reynolds numbers can be O(10) (Lauterborn and Kurz, 2010).The radical oscillation and shape modes of microbubbles and the jetting velocity are damped by viscous effects (Boulton-Stone and Blake, 1993;Smith and Wang, 2017).The compressible effects of liquid are essential, which are associated with the acoustic radiation at the minimum volumes of bubbles (Prosperetti and Lezzi, 1986;Lezzi and Prosperetti, 1987;Smith and Wang, 2018), when significant energy is radiated to the far field (Wang, 2016).This phenomenon is associated with three length scales: the thickness of the viscous boundary layer at the bubble surface, the bubble radius, and the wavelength of acoustic waves.The boundary layer thickness is smaller than the bubble radius, both changing by an order of magnitude with time, and the wavelength is in turn much larger than the bubble radius.
Bubble dynamics in a viscous liquid was simulated based on the Navier-Stokes equations using the finite volume method (FVM) (Popinet and Zaleski, 2002;Minsier et al., 2009;Hua and Lou, 2007) or finite element method (Chen et al., 2016).The compressible effects were modeled by Tiwari et al. (2013) using a diffuse interface model, and Han et al. (2015) and Lechner et al. (2017) using the FVM.However, domain numerical simulations of the three-scaled problem for dozens of cycles of oscillation or more have proven to be computationally demanding, even if feasible on supercomputers in the future.Cleve et al. (2018) evidenced the possibility of inducing steady-state shape oscillations over hundreds or thousands of oscillation cycles.Consequently, any theoretical development that can reduce the computational complexity is desirable and this creates the opportunity for a relatively simple computational analysis of a wide range of models.Our objective is to describe a new model for non-spherical microbubble dynamics in a compressible viscous flow.
Free surface flows at high Reynolds numbers are often approximated by potential flow theory.Miksis et al. (1982) included the normal viscous force at the surface of a rising bubble.Lundgren and Mansour (1988) developed the boundary layer theory for a pulsating drop.The method was later developed for bubble dynamics, by Boulton-Stone and Blake (1993) for bursting bubbles near an interface, and by Tsiglifis andPelekasis (2005, 2007) for bubbles that deform subject to an initial elongation, overpressure or acoustic disturbance, where highly deformed shapes and the details of the final stages of jet formation were captured with viscous dissipation.This rational theory is only for axisymmetric cases and tedious to be implemented.
Alternatively, Joseph and Wang (2004) introduced an auxiliary function, the viscous correction to the pressure due to potential flows, to address the shear stress not vanishing at a gas-liquid interface.We will derive a formula for the pressure correction in terms of the velocity and the shear stress of the potential flow at the interface based on conservation of energy.This model can be readily developed for threedimensional cases.
The liquid flow associated with bubbles is irrotational in the bulk volume of the liquid (Boulton-Stone and Blake, 1993).The compressible effects associated with the acoustic radiation are modeled using weakly compressible theory of Wang andBlake (2010, 2011).In the theory, the flow far away from the bubble is shown to satisfy the linear wave equation to second order in terms of the Mach number and it is obtained analytically.The flow near the bubble is shown to satisfy Laplace's equation to second order too.Wang (2013Wang ( , 2014) ) showed the computational results based on the weakly compressible theory agreed well with the experiments for underwater explosion bubbles (Hung and Hwangfu, 2010) and laser generated bubbles near a rigid boundary (Philipp and Lauterborn, 1998).
The remainder of the paper is organized as follows.The physical and mathematical model is described in Sec.II based on the weakly compressible theory, the viscous potential flow theory and the boundary integral method (BIM).In Sec.III, the formula for the viscous correction pressure is derived using conservation of energy at the interface between the liquid and gas.In Sec.IV, our numerical model is validated by comparing with Shaw's nonlinear asymptotic theory for oscillation of a bubble in a compressible and viscous liquid, Tsiglifis andPelekasis's modeling (2005, 2007) based on the viscous boundary layer and boundary integral method, and the experiments (Versluis et al., 2010) for shape oscillation of a bubble subject to ultrasound.In Sec.V, we perform a parametric analysis for microbubble dynamics in a compressible viscous liquid in terms of the amplitude and frequency of acoustic waves.The summary and conclusions are presented in Sec.VI.

II. PHYSICAL AND MATHEMATICAL MODEL A. Physical model
Consider a gas bubble suspended in a compressible and slightly viscous fluid of infinite extent, subject to an acoustic wave.The reference length, density, and pressure are chosen as the equilibrium bubble radius R 0 , the density of the undisturbed liquid q 1 , and Dp ¼ p 1 À p v , respectively.Here, p 1 is the hydrostatic pressure and p v the vapor pressure.The reference velocity is thus obtained as U ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Dp=q 1 p , and the Reynolds number Re for the liquid flow is Re ¼ R 0 ffiffiffiffiffiffiffiffiffiffiffiffi q 1 Dp p =l L .The Reynolds number is often O(10) or larger, the flow is potential in the bulk volume of the liquid except for a thin viscous boundary layer at the bubble surface, which can thus be described approximately by the viscous potential flow theory.
With the above considerations, we describe the flow in the bulk volume of the liquid as inviscid and compressible.A Cartesiancoordinate system is chosen, with the origin at the center of the initial spherical bubble, and the z-axis is along the direction of the acoustic wave, as shown in Fig. 1.The flow in the bulk volume of the liquid is governed by the continuity equation, and the Euler equation, where t is the time, u the velocity of the flow, p L the liquid pressure, and g the gravity.
The highest speed of the liquid flow induced by bubble dynamics is usually associated with the velocity of the bubble jet, which is often lower than 200 mÁs À1 at normal ambient pressure, as observed in experiments (Philipp and Lauterborn, 1998;Lindau and Lauterborn, 2003;Brujan and Matsumoto, 2012;Yang et al., 2013;Zhang et al., 2015).As the speed c of sound in water is about 1500 mÁs À1 , the flow induced by the bubble dynamics is assumed to be associated with a low Mach number e, e ¼ U c ( 1: (3) The wavelength k of incident waves or acoustic radiation is usually large compared to the bubble scale R 0 .The wavelength of ultrasound is k ¼ c/f, where f is the frequency of ultrasound, thus k is from 15 mm-150 lm for the typical range of ultrasounds of 0.1-10 MHz.The wavelength of ultrasound is usually much larger than the size of microbubbles.The wavelength k of acoustic radiation is also much larger than the bubble size, since

B. Weakly compressible theory
As the bulk volume of the liquid has two length scales, the wavelength k and the bubble radius R 0 , it is divided into two asymptotic regions: the inner region near the bubble where (x, y, z) ¼ O(R 0 ), and the outer region far away from the bubble where (x, y, z) ¼ O(k), are illustrated in Fig. 1.
Using the method of matched asymptotic expansions to the governing Eqs. ( 1) and ( 2), the outer solution was shown to satisfy the linear wave equation to second order in terms of the Mach number e, where u is the velocity potential, u ¼ ru.Strictly speaking the order of the errors in (5) is Oðe 2 Þ and this holds for the subsequent equations.
An analytical solution for the outer region was obtained as follows (Wang and Blake, 2010;Wang, 2013), where V(t) is the transient bubble volume at time t, r ¼ jrj, r ¼ (x, y, z), and b, k, and x are the amplitude, wave number, and angular frequency of the acoustic wave, respectively.The first term of the solution (6) is the incident acoustic wave and the second term is associated with the acoustic radiation due to the volume oscillation of the bubble.The solution (6) satisfies the wave Eq. ( 5), the matching conditions with the inner expansion, and the boundary condition of an incident acoustic wave at infinity, The corresponding boundary condition at infinity in terms of the pressure is where p a is the pressure amplitude of the wave.
The inner solution to the second order for Eqs.(1) and (2) satisfies Laplace's equation and the kinematic boundary condition on the bubble surface, S, as follows (Wang andBlake, 2010, 2011): To the second order, the inner flow is incompressible, which can be interpreted as follows.The wavelength k, of the incident wave or acoustic radiation, is much larger than the scale R 0 of the inner region, as such the variation of physical quantities associated with the acoustic wave over the inner region is small.The time period for an acoustic wave traveling across the inner region is of O(R 0 /c), which is much smaller than the period of bubble oscillation of O(R 0 /U).
The far-field boundary condition of the inner solution is obtained by matching with the outer solution as follows (Wang and Blake, 2010;Wang, 2016): where the first two terms are associated with the incident wave and the last term is associated with the acoustic radiation due to volume oscillation of the bubble.The initial condition on the boundary is given as where n is the unit normal at the bubble surface pointing to the gas side and R t0 is the initial rate of change of the bubble radius.

C. Viscous potential flow model
Boulton-Stone and Blake (1993) noticed that a thin viscous boundary layer exists at the bubble surface if the associated Reynolds number Re is large, whose thickness is of O R 0 = ffiffiffiffiffi Re p À Á .We will approximate the viscous effects using the viscous potential flow theory.
In the viscous potential flow theory, the normal stress balance at the bubble surface S is given as follows: where p vc is the viscous pressure correction to be discussed later, s L n the normal viscous stress, and p G the pressure of the bubble contents.We FIG. 1. Illustration of a bubble suspended in a viscous compressible liquid subject to an acoustic wave.The liquid flow domain consists of three asymptotic regions: the outer region far away from the bubble with its length scale of the wavelength k, the inner region with its length scale of the equilibrium bubble radius R 0 , and a thin viscous boundary layer with its thickness being R 0 = ffiffiffiffiffi Re p .
assume here that the viscous stress of the bubble gas is negligible and the pressure is uniform inside the bubble, since the density and viscosity of gases are two to three orders of magnitude smaller than liquids.
The normal viscous stress s L n at the bubble surface due to the potential flow is (e.g., Miksis et al., 1982) We assume that the expansion and compression of the bubble gas is adiabatic and thus the partial pressure of the bubble gas follows the adiabatic law.According to Dalton's law, the internal pressure of the bubble is where V 0 is the initial volume of the bubble, and p g0 is the initial pressure of the non-condensable bubble gas.We do not consider the thermal effects associated with this phenomenon (Szeri et al., 2003;Fuster and Montel, 2015).The effects of viscoelasticity were included in potential flow theory of microbubbles by Lind and Phillips (2010, 2012, 2013).
The tangential stress of the liquid flow at the bubble surface should approximately vanish as a result of the relatively low viscosity of the gas inside the bubble.However, the shear stress of a potential liquid flow is non-zero at the bubble surface.As shown in Sec.III, the viscous pressure correction p vc to the pressure due to potential flows can be introduced to address this discrepancy based on conservation of energy.It is given in terms of the normal velocity u n , tangential velocity u s , and shear stress s L s due to the potential flow at the bubble surface as follows: Using the Bernoulli equation, the dynamic boundary condition (10) at the bubble surface can be written as where g is the acceleration of gravity.
For the axisymmetric case, p vc ¼ Àu s s L s =u n .As the above equation is invalid for u n ¼ 0, p vc is thus set to be zero as j u n /Uj 0.01 in our calculations.Noticing s L s ¼ l L ð@u n =@s þ @u s =@nÞ and the flow is irrotational, i.e., @u n =@s À @u s =@n ¼ 0, the shear stress s L s is (Boulton-Stone and Blake, 1993) where j 1 is the curvature of the intersection curve of the bubble surface with the azimuthal plane, and s is the arc length parameter of the curve.Denoting the intersection curve as: r c (s) and z c (s), we have where the overdot denotes a derivative with respect to s.
Substituting ( 12)-( 15) into ( 16) yields Examining the initial and boundary value problem of ( 9), one can see that the compressible effects and the viscous effects appear in the dynamics boundary condition at the bubble surface (17).As the basic equation is Laplace's equation, this problem can be modeled using the boundary integral method (BIM).The details on the numerical model using the BIM for this problem can be found in (Wang et al., 1996a(Wang et al., , 1996b;;Curtiss et al., 2013).The computations have been carried out for axisymmetric cases; however, the numerical model can be developed for three-dimensional cases in a straightforward manner, by updating the potential using ( 14) instead of ( 17).The effects of buoyancy is usually negligible in for microbubbles applications, since the associated buoyancy parameter A non-spherical bubble collapse often leads to the formation of a high-speed liquid jet.The jet subsequently impacts the opposite bubble surface and penetrates the bubble, and the liquid domain is then transformed from a singly connected to a doubly connected domain.The solution to a potential problem in a doubly connected domain is nonunique.The doubly connected domain can be made singly connected by, for example, using a branch cut by Best (1993) or a vortex sheet by Zhang et al. (1993) and Zhang and Duncan (1994).Wang et al. (1996b) developed a vortex ring model from these earlier ideas to model the topological transition of a singly connected bubble to a subsequent toroidal bubble.In the vortex ring model, a vortex ring is put inside the toroidal bubble after jet impact.The circulation of the vortex ring is equal to the jump of the potential u across the contact point at the time of jet impact, where u N and u S are the potentials at the impact point.Here, we assume that jet impact occurs at a single point.The potential u is then decomposed as follows: where u vr is the potential of the vortex ring, which can be obtained from the Biot-Savart law (Wang et al., 1996b;Liu et al., 2016).With the potential jump being accounted for by the vortex ring using ( 19), the remnant potential / is continuous in the flow field and can be simulated using the BIM model.

III. VISCOUS PRESSURE CORRECTION
In the classical potential flow theory, the fluids are assumed to be inviscid in order to guarantee that the flow is irrotational.However, there are many situations where the viscosity of fluids is not small but the vorticity is small and not important, which may be well approximated by a potential flow model.The typical cases studied using the viscous potential flow theory are two-phase gas-liquid flows associated with an interface.
A key issue is that the shear stress due to the irrotational flow does not vanish at the interface, which contradicts the physical boundary condition.This results in the violation of the conservation of energy, as to be shown later on.A viscous pressure correction to the irrotational pressure was introduced to resolve this discrepancy, initially regarding the drag on a rising spherical bubble at a constant velocity by Moore (1959Moore ( , 1963)).Kang andLeal (1988a, 1988b) showed that the viscous pressure correction at the surface of a rising spherical bubble can be expressed in terms of the spherical harmonic series.The viscous potential flow theory has been applied for many cases, including a rising spherical cap bubble (Davies and Taylor, 1950;Joseph, 2003), the decay and stability of free surface waves (Wang and Joseph, 2006;Wang et al., 2005), and the Kelvin-Hemholtz instability (Joseph, 2006;Padrino et al., 2011).

A. Conservation of energy at a free surface
Consider a two-phase Newtonian flow of liquid and gas.The mechanical energy of a compressible Newtonian flow of the liquid side follows: where q and u are the fluid density and velocity, respectively, f is the body force per unit mass, e is the rate of the strain tensor, and r s is the stress tensor, where 1 is the unit tensor, p the pressure, and l the viscosity the flow.Equation ( 20) can be obtained from the Navier-Stokes equations and continuity equation for a compressible flow.An arbitrary material control region X around an area S I on the interface and slightly into the liquid side is prescribed, as shown in Fig. 2. The boundary R of the region X consists of S I , S L , and a narrow surface S J joining S I and S L .
We apply the integral form of the energy Eq. ( 20) to the material control region X, where T is the stress We replace the volume integrals in ( 22) by an average value times the volume of integration region as follows, based on the continuum assumption: We limit the region toward the interface so that the volume X and the narrow joining surface S J vanish, and only the surfaces S I and S L remain, S L ! S I , n L ! -n, where n and n L are the unit outward normal at the surfaces S I and S L , respectively.Providing the velocities and their partial derivatives are continuous, (23) becomes where T L is the stress at S L due to the liquid flow, and T I is the stress at the interface due to the gas flow.Since the integration surface S I is arbitrary, the integrand must be zero, We make the following approximations: where p L , s L n , and s L s , and p G , s G n , and s G s are the pressure, normal viscous stress, and shear stress of the liquid and gas at the interface, respectively.As the viscosity of gases is much smaller than liquids, the normal viscous stress and shear stresses s G n and s G s are negligible.Substituting (26) into (25) yields This equation is satisfied at the free surface in the viscous modeling with the balance of the normal stress and vanishing of the shear stress, but it is not satisfied in the potential flow theory, where s at the interface is usually non-zero.Accordingly, the energy is not conserved at the interface in potential flow theory.

B. Viscous pressure correction
Now we introduce the viscous potential flow theory.It is assumed that the flow is irrotational in the flow domain.It satisfies Eq. ( 27) at the gas-liquid free surface to satisfy the conservation of energy.Equation ( 27) can be rewritten as Introducing a viscous correction pressure p vc as follows: Equation ( 28) becomes or FIG. 2. Arbitrary material control region control region X on the liquid side around an area S I on the interface between liquid and gas.

Physics of Fluids
Adding a viscous correction pressure p vc of (29) in the balance of normal stress (31) at the interface leads to the satisfaction of (28) as well as ( 27) and ( 22), the energy conservation at the interface.The normal stress Àp L þ s L n due to the liquid potential flow outside the viscous boundary layer is corrected by adding p vc to include the viscous effects of the boundary layer and to conserve the energy at the interface.With the inclusion of surface tension, (31) becomes The integral form of ( 29) on the interface I can be written as ð This relation was introduced by Joseph and Wang (2004).Relation ( 33) is a foundation for and widely used in the viscous potential flow theory (c.f.Joseph and Wang, 2004;Wang and Joseph, 2006;Wang et al., 2005;Joseph, 2006;Padrino et al., 2011).It has led to improvements to the potential flow theory for many problems, including the rising of a spherical bubble/drop, the decay of free gravity waves on water and the Kelvin-Helmholtz instability.Motivated by these great successes, we have developed their hypothesis locally for potential flow theory.The viscous pressure correction (29) can be used for three-dimensional free surface flows, but it does not apply when the viscous boundary layer separates into the liquid bulk.
A remarkable example of the viscous pressure correction is for a spherical gas bubble of radius R rising in a viscous liquid at high Reynolds number.Levich (1949) obtained the drag on the bubble at 12pRlU, where U is the rising velocity of the bubble, by calculating the dissipation of the irrotational flow around the bubble.Moore (1959) calculated the drag directly by integrating the pressure and viscous normal stress of the potential flow and neglecting the viscous shear stress, obtaining the value 8pRlU.This approach results in one third of the relative error, because the energy of the flow is not conserved at the free surface in potential flow theory.The total energy is conserves as well, since the energy is conserved in the inner flow domain via Bernoulli's equation.Joseph and Wang (2004) obtained the correct value for the drag, 12pRlU, using the potential flow theory with the viscous pressure correction.The viscous pressure correction guarantees the conservation of the energy at the free surface as well as globally.

IV. VALIDATIONS
In this section, comprehensive validations will be carried out for the viscous compressible potential flow theory, because it is new.It will first be compared with the Keller-Miksis equation for the oscillation of spherical bubbles.It will be then compared with Shaw's nonlinear asymptotic theory (2017) for nonspherical bubbles in a compressible and viscous liquid as well as Tsiglifis and Pelekasis's model ( 2005) based on the viscous boundary layer and boundary integral method.It will be further compared with experiments for the dynamics of non-spherical bubbles.We also check the convergence of the results in terms of the mesh size.
We first compare the results obtained from the viscous compressible BIM (VCBIM) and the Keller-Miksis equation (KME).The case considered is a bubble having an equilibrium radius of 26 lm suspended in water subject to an acoustic wave with the pressure amplitude of 20 kPa and frequency of 130 kHz. Figure 3 compares the time histories of the bubble radius obtained from the VCBIM and KME, respectively.The two results have excellent agreement for the first nine cycles of oscillation.
Next, the VCBIM is compared with other numerical results.Shaw (2017) developed a nonlinear asymptotic model to study the shape mode oscillation of parametrically forced bubbles including viscous and compressible effects, which accounts for nonlinear shape mode interactions to third order.Many cycles of oscillation were considered, showing the importance of the inclusion of viscous and compressible effects.The case chosen for comparison considers a microbubble with an equilibrium radius of 144 lm subject to an acoustic wave with a pressure amplitude of 13 kPa and a frequency of 10 kHz.The remaining parameters are the same as in Fig. 3.As is seen in Fig. 4, the VCBIM obtains excellent agreement with the results of the asymptotic model.Non-spherical oscillation is accurately predicted at each time, despite the numerous cycles of oscillation leading up to this time.The bubble displays mixed shape modes, with mode 3 being predominant at t ¼ 10.07 [Fig.4(a)] and t ¼ 10.17 [Fig.4(c)] and mode 6 predominant being at t ¼ 10.11 [Fig.4(b)] and t ¼ 10.22 [Fig.4(d)].This demonstrates the accuracy and robustness of the VCBIM to model bubble dynamics for dozens of cycles of oscillation with important viscous and compressible effects.This also shows that Shaw's nonlinear asymptotic model is accurate for non-spherical oscillations of bubbles at large amplitude.
Additionally, Tsiglifis and Pelekasis (2005) examined the weakly viscous oscillation of elongated bubbles using the viscous boundary layer theory coupled with the boundary integral method.They revealed that small initial elongations would return to a spherical shape for any Ohnesorge number, Oh ¼ l=ðqRrÞ 1=2 , while for larger elongations there is a threshold value of Oh À1 above which the bubble breaks up.The case chosen for comparison considers a microbubble with an equilibrium radius of 5.8 lm and an initial elongation parameter of S ¼ 0.6.The elongation is defined as S ¼ a/R eq0 with a and R eq0 FIG. 3. Comparison of the results obtained from the VCBIM and KME for the radius history R(t) for a bubble having an equilibrium radius of 26 lm subject to an acoustic wave with the amplitude 20 kPa and frequency 130 kHz for nine cycles.The other parameters are j ¼ 1.4, r ¼ 0.073 N m À1 , p 1 ¼ 100 kPa, p v ¼ 3 kPa, and q L ¼ 1000 kg m À3 .denoting the length of the smaller semi-axis and the radius of a bubble with the same initial volume, V 0 , as the elongated bubble.The inverse Ohnesorge number is given by Oh À1 ¼ 1000, with the rest of the parameters as in Fig. 3. Figure 5 shows the comparison between the results of Tsiglifis and Pelekasis (2005) with those of the VCBIM during jet formation.Very good agreement is achieved between the bubble shapes, with good agreement between the jet profiles.We next compare the numerical results of the VCBIM with experiments.Versluis et al. (2010) carried out a series of experiments for the shape oscillation of a microbubble driven by an ultrasonic wave.They revealed shape oscillations for various mode numbers n ¼ 2-6 of microbubbles, with the ultrahigh-speed imaging.They found that the mode number n is dependent on the bubble radius but is independent of the pressure amplitude.The experimental case chosen for the comparison is for amplitude of 120 kPa and frequency of 130 kHz. Figure 6 shows the comparison of the bubble shapes using the experimental images and the computational results, at successive maximum and minimum bubble volumes.The wave propagates from the left to the right side for this case as well as the subsequent cases.

Physics of Fluids
Figure 6(a) shows the bubble shapes during the first 5 cycles of oscillation.The bubble is in equilibrium before the arrival of the acoustic wave (t ¼ 0 ls).It starts to collapse as the acoustic wave defined by ( 8) is initially positive at the bubble location.During the first 5 cycles, the bubble oscillates in a spherical shape with the amplitude increasing with each cycle.Figure 6(b) shows the bubble shape from the sixth to ninth cycles of oscillations.At the sixth maximum volume (t ¼ 48.25 ls), the bubble is approximately spherical, with its right side slightly over protruding.However, the bubble quickly becomes non-spherical during the following collapse, when its cross section takes a square shape at the sixth minimum volume, with rounded corners and a weak jet forming on the right side.The jet further develops as the bubble expands to the seventh maximum volume.The surface mode n ¼ 4 becomes obvious at and after the seventh minimum volume.The right jet remains as the bubble oscillates during the eighth cycle.At the end of the eighth minimum volume, a left jet starts to develop and becomes obvious at the ninth maximum volume, when two opposing jets occur along the axis of symmetry.
The computational results agree well with the experimental images for all nine cycles of oscillation, in terms of the bubble shapes and the time sequence.All the above features were reproduced by the computation.However, the jet is not visible in the experimental images due to the opaqueness of the bubble.case shown in Fig. 6.The numerical results for m ¼ 61 and 71 have excellent agreement for the first 9 cycles of oscillation.The remaining calculations in this paper were undertaken for m ¼ 61.

V. NUMERICAL RESULTS AND DISCUSSIONS
In this section, we perform a parametric analysis of the oscillation of a microbubble suspended in an infinite liquid subject to acoustic waves of various amplitudes and frequencies.Particular attention is paid to the bubble oscillations when the driving frequency of acoustic waves is set to be the natural frequencies of either spherical oscillation or shape oscillation.

A. Effects of the pressure amplitude of ultrasound p a
To study the effects of the pressure amplitude of acoustic waves, we consider a bubble with the equilibrium radius R 0 ¼ 30 lm, subject to an acoustic wave having the frequency f ¼ 130 kHz and the various pressure amplitudes p a ¼ 40, 47, and 50 kPa.The remaining parameters are the same as in Fig. 6.For the lower pressure amplitude p a ¼ 40 kPa the bubble remains approximately spherical as far as the 23th cycle of oscillation, with the spherical bubble shapes at the 23th minimum and 23th maximum volumes being shown in Fig. 8(a).
For the intermediate pressure amplitude p a ¼ 47 kPa, the bubble remains approximately spherical for the first 13 cycles of oscillation.As shown in Fig. 8(b), the bubble is still approximately spherical at the 14th minimum and 14th maximum volumes but its left side, facing the wave propagation direction, is slightly flattened.The right side is, in turn, flattened at the 15th minimum and maximum volumes, while its left side returns to a spherical shape.A jet starts forming at the minimum volume.The turnover of jetting at the minimum volume and the flattened left and right sides repeats and enhances during subsequent cycles.Surface mode n ¼ 3 becomes obvious gradually from 14th oscillation to 19th oscillation.During the 18th and 19th cycles of oscillation, shape mode 3 has become well developed.The cross section of the bubble roughly takes the form of an equilateral triangle with one of its sides alternating between either facing or backing onto the wave direction during successive cycles of oscillation.The jet also alternates between the left and right sides during successive cycles, facing the acoustic wave and along the direction of the acoustic wave.
For the larger pressure amplitude p a ¼ 50 kPa the bubble is approximately spherical for the first 8 cycles of oscillation.Its left side is flattened and the right side elongated at the ninth minimum and ninth maximum volumes.The jet forms at the right side at the tenth minimum volume and at the left side at the 11th minimum volume.Subsequently, the jet alternates between the left and right sides.The surface mode n ¼ 3 of the bubble becomes obvious in cycle 12, and continues developing.During cycles 13 and 14, jets occur both at the minimum and maximum values.
The bubble initially oscillates in a purely spherical mode and surface modes can be generated after several acoustic cycles if the acoustic pressure is above a critical threshold.When the pressure amplitude is slightly higher than the threshold, the shape mode develops gradually for many cycles of oscillation.The parametric instability for a bubble due to an acoustic wave is a cumulative effect, requiring many oscillation cycles to build up.As the pressure amplitude increases, the development of shape modes starts earlier and grows faster, and the subsequent shape mode has larger amplitude.This is because the Bjerknes force F B acting on the bubble due to the pressure wave p ac (r, t) is proportional to the pressure amplitude.The Bjerknes force F B is given by where r c is the geometrical center of the bubble and z c is its z-coordinate.Here, we used the pressure wave given in (8).

B. Effects of the frequency of ultrasound
To consider the effects of the driving frequency of ultrasound, we repeat the case shown in Fig. 8(b) for the acoustic pressure amplitude p a ¼ 47 kPa and the equilibrium radius R 0 ¼ 30 lm with different frequencies, the remaining parameters being kept the same.The natural frequency for spherical oscillation for a bubble is given as (Brennen, 1995) Using R 0 ¼ 30 lm, j ¼ 1.4, r ¼ 0.073 NÁm À1 , p 1 ¼ 100 kPa, p v ¼ 3 kPa, and q L ¼ 1000 kgÁm À3 , we have f 0 ¼ 110 kHz.We want to compare the three cases for f ¼ 85 kHz < f 0 , f ¼110 kHz ¼ f 0 , and f ¼ 130 kHz > f 0 .The last case for the driving frequency being larger than the natural frequency was considered in Fig. 8(b).For f ¼ 85 kHz, the driving frequency is smaller than the natural frequency f 0 of the bubble.The bubble remains a spherical shape till the 24th cycle, as shown in Fig. 9(a).
For f ¼ 110 kHz, the wave frequency is equal to the natural frequency f 0 of the bubble, the bubble undergoes resonant oscillation.As shown in Fig. 9(b), the bubble reaches a much larger maximum volume during the first and second cycles as compared to f ¼ 85 or 130 kHz [shown in Figs.9(a) and 8(b), respectively].Consequently, a much larger Bjerknes force acts on the bubble, since it is proportional to the bubble volume as shown in (34).During the third cycle of oscillation, a reentry liquid jet forms at the bubble side facing the acoustic wave at the minimum volume and it retakes a spherical shape at the subsequent maximum volume.This repeats during the fourth cycle of oscillation, when the jet develops further.Subsequently, the nonspherical oscillation develops quickly.The bubble shape at the fifth maximum volume (t Ã ¼ 15.73) is elongated along the wave direction.Two opposing jets develop and impact each other at the sixth minimum volume.The bubble then rejoins before reaching the sixth maximum volume.In subsequent time, the singly connected bubble continues to oscillate.
The bubble remains spherical as f ¼ 85 kHz < f 0 [Fig.9(a)], becomes non-spherical after 15 cycles of oscillation for f ¼ 130 kHz > f 0 [Fig. 8(b)], but undergoes obvious non-spherical oscillation even from the third cycle if the driving frequency is equal to the natural frequency of the bubble for f ¼ 110 kHz ¼ f 0 , when the bubble undergoes resonant oscillation.

C. Driving frequency equal to resonance frequency of shape mode
The natural frequency f n of shape modes n of bubbles is given as (Lamb, 1932) where n >¼2.For R 0 ¼ 39 lm, r ¼ 0.073 N m À1 and j ¼ 1.4, the natural frequencies for shape modes of the bubble for n ¼ 3, 4, and 6 are f 3 ¼ 70.6 kHz, f 4 ¼ 145 kHz, and f 6 ¼ 193 kHz, respectively.We are now considering the bubble suspended in a liquid subject to an acoustic wave at the natural frequencies of shape modes, when the shape instability of the bubble is prone to appear.We first consider the case for the driving frequency being set at the natural frequency of mode 3: f ¼ f 3 ¼ 70.6 kHz and the amplitude p a ¼ 40 kPa.As shown in Fig. 10, the bubble oscillates for nine cycles in a spherical volumetric mode.The bubble starts developing surface mode 3 during the tenth cycle of oscillation, the right part of the bubble is flattened at the tenth minimum volume and the left part becomes flattened at the tenth maximum volume (t Ã ¼ 36.97).This is FIG.9. Dynamics of a bubble with an equilibrium radius R 0 ¼ 30 lm driven by an acoustic wave with the pressure amplitude p a ¼ 47 kPa and frequencies: (a) f ¼ 85 and (b) f ¼ 110 kHz, respectively.The dimensionless time is shown on the upper right corner in each frame.The remaining parameters are the same as in Fig. 6.
reversed during the subsequent period, when the left side becomes flattened at the 11th minimum volume, where a jet forms, and the right part becomes flattened at the 11th maximum volume.Subsequently, the above features repeat and mode 3 develops and becomes obvious in the 19th cycle of oscillation.
Figure 11 shows the development of shape mode 4 of the bubble as the driving frequency is equal to the corresponding natural frequency f 4 ¼ 145 kHz and the pressure amplitude p a ¼ 75 kPa.The bubble becomes non-spherical at the third minimum volume, when both left and right sides become flattened with weak jets, but is approximately spherical at the third maximum volume.During the fourth cycle of oscillation, the bubble is elongated along the wave direction at both minimum and maximum volumes.At the fifth minimum and maximum volumes, shape mode 4 becomes obvious, when the bubble cross section takes a square shape with rounded corners and its sides parallel or perpendicular to the wave direction.At the sixth minimum volume, the bubble has a square cross section, with its diagonals now being parallel or perpendicular to the wave direction.A larger expansion of the bubble takes place in this cycle and the bubble regains a spherical shape at the sixth maximum volume.From seventh to eighth cycle of oscillation, the surface mode grows in magnitude.During the ninth cycle, two opposing jets develop along the wave direction and impact each other at the ninth minimum volume.Subsequently, the bubble becomes toroidal but it still displays the character of mode 4.
Figure 12 shows that the bubble displays shape mode 6 as the driving frequency is equal to the corresponding natural frequency f ¼ f 6 ¼ 194 kHz and the pressure amplitude p a ¼ 235 kPa.The surface mode n ¼ 6 develops at the 13th minimum volume and it has a square shape at the 13th maximum volume.For the next four cycles, the bubble displays the surface mode n ¼ 4.During 19th and 20th oscillation, the bubble displays surface mode n ¼ 6.As predicted by Doinikov (2004), even modes can interact with each other.Energy can be transferred back and forth between the two modes.

D. Viscous and compressible effects
To analyze the viscous and compressible effects, we compare the computational results with and without these effects.The importance of viscous effects is demonstrated in Fig. 13, where the bubble shapes of the viscous and inviscid models are compared during bubble collapse.This is the case considered in Fig. 5 for the oscillation of an elongated bubble.Including viscous effects leads to the formation of a wider jet, at both the opening of the jet and the tip.This is as expected when considering viscous effects and agrees with the results of Tsiglifis and Pelekasis (2005).
Four variations of the case described in Fig. 6 are considered, including (a) both viscous and compressible effects, (b) only viscous effects, (c) only compressible effects, and (d) neither viscous nor compressible effects.The corresponding bubble shapes are shown in Fig. 14 at the ninth minimum bubble volume, where the jet is clearly weakened due to the viscous effects.However, there are no significant dissipative compressible effects for this case, since there is no significant collapse here.This is consistent with the observation of Calvisi et al. (2007).
The ninth minimum occurs at different times in Fig. 14 depending on whether or not viscous effects are included.Compressible effects do not modify the time for the ninth minimum.The frequency FIG. 10.Development of shape mode 3 of a bubble with an equilibrium radius R 0 ¼ 39 lm subject to an acoustic wave with the pressure amplitude p a ¼ 40 kPa, and the driving frequency being equal to the natural frequency of mode 3 of the bubble, f ¼ f 3 ¼ 70.6 kHz.The dimensionless time is shown on the upper right corner in each frame.The remaining parameters are the same as in Fig. 6. and the phase shift for the bubble oscillations are modified by viscosity (Smith and Wang, 2017) but not by compressibility (Smith and Wang, 2018).

VI. SUMMARY AND CONCLUSIONS
In potential flow theory for free surface flows with large Reynolds number, the key issue is that the shear stress should approximately vanish at a gas-liquid interface, but it does not.This results in the violation of the conservation of energy.We have derived a formula for the viscous pressure correction in terms of the velocity and the shear stress of the potential flow at a free surface, to enforce the conservation of energy.The formula approximates directly the local viscous effects at an interface and is applicable for three-dimensional free surface flows.We have developed an accurate, grid free in the flow domain and robust numerical model for a bubble oscillation over many cycles involving topological transformation, with the viscous pressure correction as well as the weakly compressible theory and vortex ring model.
The numerical results are shown to be in good agreement with Shaw's nonlinear asymptotic theory (2017) for bubbles in viscous and compressible liquids, the Tsiglifis and Pelekasis model (2005) based on the viscous boundary layer theory and boundary integral method for the oscillation of elongated bubbles, as well as experiments for microbubbles undergoing shape oscillation when subject to ultrasound.
The computations have been carried out for axisymmetric cases however the numerical model can be developed for three dimensional straightforwardly.We observed the following features for microbubble dynamics subject to a pressure wave.
A bubble subject to an acoustic wave first oscillates in a spherical shape due to surface tension.Shape modes develop gradually if the wave amplitude is beyond a critical threshold, this being a supercritical bifurcation.The parametric instability for a bubble due to an acoustic wave is a cumulative effect, requiring many oscillation cycles to build up.As the wave amplitude increases, the shape modes start earlier and develop faster, because the Bjerknes force is proportional to the wave amplitude.
If the wave frequency is near to the natural frequency of spherical oscillation, the bubble undergoes resonant spherical oscillation at large amplitude and shape mode 3 develops quickly.The large Bjerknes FIG.11.Development of shape mode 4 of a bubble with an equilibrium radius R 0 ¼ 39 lm subject to an acoustic wave with the pressure amplitude p a ¼ 75 kPa, and the driving frequency being equal to the natural frequency of mode 4 of the bubble, f ¼ f 4 ¼ 145 kHz.The dimensionless time is shown on the upper right corner in each frame.The remaining parameters are the same as in Fig. 6.
force along the direction of acoustic wave propagation, which is proportional to the bubble volume, is responsible for the growth of shape mode 3.
If the driving frequency is equal to the natural frequencies of different shape modes of a bubble, the bubble is excited into the corresponding shape modes.The orientation of the shape mode and the jetting are determined by the direction of propagation of the driving acoustic wave.For shape mode 3, the bubble takes a triangular shape with the sides facing and backing onto the wave direction being flattened in turn during each successive cycle of oscillation, where reentry liquid jets often form parallel to the wave direction.For shape mode 4, the bubble takes a square shape with its sides and diagonals being parallel to the wave direction during each successive cycle of oscillation.If the driving frequency is equal to the natural frequency of shape mode 6, the bubble alternates between shape modes 4 and 6.This is associated with the exchange of energy between even shape modes, it being due to the nonlinearities in the kinematic and normal stress conditions at the free surface.
The standard inviscid model and BIM capture the essential shape of the bubble, but prediction of the amplitude of the bubble oscillation requires the conservation of energy including viscosity.Viscous effects FIG.12. Development of shape mode 6 of a bubble with an equilibrium radius R 0 ¼ 39 lm subject to an acoustic wave, with the driving frequency equal to the natural frequency of mode 6 of the bubble, f ¼ f 6 ¼ 194 kHz and the pressure amplitude p a ¼ 235 kPa.The dimensionless time is shown on the upper right corner in each frame.The remaining parameters are the same as in Fig. 6.FIG. 13.Comparison between the VCBIM (solid line) and Tsiglifis and Pelekasis (2005) (dotted line) during bubble collapse at t ¼ 1.3453.The parameters are the same as in Fig. 5, except Oh À1 is infinite in Tsiglifis and Pelekasis (2005).
result in a reduction of the amplitude and a modification of the phase shift for strongly nonlinear non-spherical bubble oscillations in comparison with the inviscid fluid flow.The reduction in the amplitude also results in changes to the oscillation frequency which is dependent on amplitude in the finite-amplitude case.Furthermore, during the collapse, viscosity dampens the jet formation, which leads to flatter bubble jets.

Figure 7
FIG. 4.Comparison of the results obtained from the VCBIM (red dotted line) with those ofShaw (2017) where a bubble with equilibrium radius 144 lm is subject to an acoustic wave with a pressure amplitude of 13 kPa and a frequency of 10 kPa.The selected times are at (a) t ¼ 10.07;(b) t ¼ 10.11; (c) t ¼ 10.17; and (d) t ¼ 10.22 ms.The other parameters are j ¼ 1.4, r ¼ 0.073 N m À1 , p 1 ¼ 100 kPa, p v ¼ 3 kPa, and q L ¼ 1000 kg m À3 .

FIG. 7 .
FIG. 7.Convergence test in terms of the number m of mesh elements for the time histories of the equivalent bubble radius R eq using the CVBIM for the case shown in Fig.6.