IRROTATIONAL FLOW DUE TO FORCED OSCILLATIONS OF A BUBBLE

Abstract

The behaviour of an axisymmetric bubble in a pure liquid forced by an acoustic pressure field is analysed. The bubble is assumed to have a sharp deformable interface, which is subject both to surface tension and to Rayleigh viscosity damping. Two modelling regimes are considered. The first is a linearized solution, based on the assumption of small axisymmetric deformations to an otherwise spherical bubble. The second involves a semi-numerical solution of the fully nonlinear problem, using a novel spectral method of high accuracy. For large-amplitude nonspherical bubble oscillations, the fully nonlinear solutions show that a complicated resonance structure is possible and that curvature singularities may occur at the interface, even in the presence of surface tension. Rayleigh viscosity at the interface prevents singularity formation, but eventually causes the bubble to become purely spherical unless shape-mode resonances occur. An extended analysis is also presented for purely spherical bubbles, which allows for a more detailed study of the effects of resonance and the Rayleigh viscosity at the bubble surface.

1. Introduction

The stability, morphology and behaviour of spherical and nonspherical bubbles have been of immense interest for the past century. Bubbles occur in natural environments, in oceans and in the interstellar medium surrounding supernovae explosions [27], as well as in a vast range of technical applications such as in refrigeration equipment [1], ultrasonics and medicine. Early studies of bubbles were conducted by Besant and Rayleigh whose research into the collapse of spherical cavities led to the ubiquitous Rayleigh equation for a bubble of gas in a liquid [37]. Later, Plesset [32], Noltingk and Neppiras [30], Gilmore [11] and others generalized this model by adding viscosity, surface tension and an acoustic field. Assuming that the bubble remains spherical, the inclusion of viscosity only affects the boundary condition, not the governing equation itself, since purely radial fluid flow remains irrotational even when viscosity is present [34]. To account for energy loss from the radiation of acoustic energy, Keller and Kolodner [19] allowed the fluid surrounding the bubble to be weakly compressible by solving the wave equation rather than Laplace’s equation. This leads to the damped radial oscillations. Keller and Miksis [20] later added viscosity, surface tension and an acoustic field to this compressible model, and found periodic solutions even for high forcing amplitudes.

Complex models of single-bubble sonoluminescence can include heat and mass exchange through the bubble surface, gas diffusion, changes of state, chemical reactions and shocks [29]. This topic is thoroughly reviewed by Brenner et al. [4]. The radial oscillations of a forced spherical bubble are highly nonlinear; superharmonic resonances are common, and subharmonics and ultraharmonics can occur once a pressure threshold has been reached [12, 20]. Period doubling and chaos can occur, and multiple solution states can coexist for a given set of parameters leading to hysteresis [9]. These behaviours are often studied using resonance curves and bifurcation diagrams [41].

Surface tension and viscosity often work to make bubbles more spherical. However, if the bubbles are sufficiently large, or reside in a strong gravitational field, or are subject to disturbances created by other bubbles, rigid surfaces or fluid flow, the assumption of sphericity is unlikely to hold true [33, 35, 48]. Jet formation, bubble break-up or pinch-off and microstreaming can result from nonspherical disturbances [6, 18, 46]. Nonspherical effects are important in the study of sonoluminescence as the acceleration of the bubble wall during bubble collapse amplifies nonspherical shape deformations and can lead to bubble break-up [22, 39]. This is a Rayleigh–Taylor style instability, and it can also occur during afterbounces following the initial collapse [22].

Nonlinear interactions between shape and volume oscillation modes play an important role in these nonspherical bubble behaviours. The transfer of energy from radial oscillations to shape oscillations can lead to bubble break up [26]. Acoustic waves can cause nonlinear volume oscillations large enough to spark shape deformations and erratic bubble motion [39]. A translation of the bubble caused by travelling acoustic waves can also lead to shape deformation and volume oscillations. Forcing a bubble near its natural volume resonance can cause large volume oscillations that induce shape oscillations [26]. These shape oscillations can then affect the natural volume resonance, meaning that the bubble is no longer being forced at its natural resonance. This reduces the amplitude of volume oscillations, and the bubble returns to its initial state and volume resonance [9]. In this way, volume and shape oscillations can appear to compete. Anisotropic forcing pressures and near, but not exact, 2:1 resonance between volume and shape modes can also cause this continuous energy transfer between volume and shape modes [8]. This effect can also occur without external forcing in the absence of damping. If the 2:1 resonance is exact, an initial shape disturbance will be transferred completely to volume oscillations [10]. However, if there is a sufficient degree of detuning away from exact resonance, or a background extensional flow, energy can flow continuously between the modes [10, 47]. This resonance means that shape distortions can cause volume oscillations, and thus sound production at twice the frequency of the shape oscillations [25].

The upwards movement of a bubble due to gravity can also be a source of asymmetry. Large bubbles of radius R can become oblate and rise in erratic paths [45]. Very large bubbles tend to form spherical caps, sometimes with ragged skirts, and the speed at which they rise is proportional to $\sqrt {gR}$ [2]. To investigate bubble behaviour without gravity, low or zero gravity environments can be created in free-fall towers [42] or on parabolic flights [48]. Zero-gravity environments are relevant to many space systems and applications. For smaller bubbles, gravity has less effect, and surface tension and viscosity are better able to maintain sphericity of the bubble. Bubbles with radii less than approximately 0.1 mm generally rise in straight lines and their speed is proportional to $gR^2/\nu $ , where $\nu $ is the kinematic viscosity in the surrounding fluid [21]. For very small bubbles, the effects of gravity are minimal, and gravity is often ignored, as is the case in this paper.

One powerful technique for analysing fluid flow is the use of a potential. This use is often limited to the analysis of inviscid fluids; however, many of the desirable viscous effects, such as damping and dissipation, do not prevent the use of a potential flow as long as the flow remains irrotational. Further discussion of potential flows of viscous fluids may be found in Joseph [16]. Poritsky’s introduction of viscosity without rotation to purely spherical bubbles (see [16]), Shaw’s use of a Rayleigh dissipation function for axisymmetric bubbles [39], and Joseph and Wang’s use of a viscous pressure correction to account for the shear stress at the interface of a bubble [17] are all techniques used for this purpose in the study of bubbles.

In this paper, we consider an incompressible and irrotational liquid but model some viscous effects using an artificial Rayleigh viscosity at the surface of the bubble. Rayleigh introduced the concept of an artificial viscosity when he considered the standing waves formed around a disturbance in uniformly flowing water [36]. This idealized viscosity is incorporated in the dynamic boundary condition, and acts proportionally to the relative velocity to resist the movement of particles. It mimics the dissipative effects of real kinetic viscosity without introducing rotation. This technique is still used to analyse standing waves generated by moving disturbances [31]. It is also used in other problems such as for resonant sloshing in baffled tanks, where it provides necessary damping effects which stop divergence at resonant frequencies [5]. It is these damping and dissipative effects that are useful in our current study.

The remainder of this paper is structured as follows. In Section 2, we describe our nonspherical bubble model; we start with a simple linear case, and then formulate a numerical nonlinear model which is solved using spectral methods. The effects of surface tension and Rayleigh viscosity at the interface are discussed. This work builds upon an earlier paper by the authors, which considered both an inviscid model and a Boussinesq model [7]. The addition of the Rayleigh viscosity aids in the prevention of curvature singularities which previously occurred in the nonlinear inviscid model. This nonspherical analysis allows for interactions between shape and volume oscillations, and for resonances in shape mode oscillations. Section 3 continues with a nonlinear spherical extension. Results are given in Section 4, where the behaviour of the models is compared, and the effects of the Rayleigh viscosity, the forcing frequency and amplitude of the acoustic pressure field are explored. The paper closes with a few conclusions and final remarks.

2. Nonspherical forced oscillations

We consider an axisymmetric nonspherical bubble of gas in a liquid. The gas is assumed to be ideal, isentropic and diatomic, the liquid incompressible, and the flow irrotational. An artificial Rayleigh viscosity provides an approximate viscosity at the bubble interface, accounting for the dissipative effects of viscosity but not the swirling rotational effects. This simple Rayleigh viscosity still allows for Laplace’s equation to be used as the governing equation in place of the full Navier–Stokes equations. Analysis takes place in a spherical coordinate system where $(r,\phi ,\theta )$ are the distance from the origin, and the polar and azimuthal angles, respectively. The bubble is assumed to be axisymmetric about the z-axis, so that $ \partial /\partial \theta \equiv 0$ and all variables are independent of the azimuthal angle $\theta $ . There are two variables of interest, the radius of the bubble given by $R(\phi ,t)$ and the velocity potential $\Phi (r,\phi ,t)$ . Natural scaling variables can be deduced from a simplified spherical model. In the absence of surface tension and viscosity, the bubble has an equilibrium radius $R_0$ and a natural frequency $\omega _{\mathrm {eq}}$ given by

$$ \begin{align*} R_0 = \bigg[\frac{C}{p_\infty}\bigg]^{{1}/{3\gamma}} \quad \text{and} \quad \omega_{\mathrm{eq}}^2 = \frac{3\gamma p_\infty}{\rho R_0^2}, \end{align*} $$

where $C = K(3M_0/4\pi )^\gamma $ , in which K is a constant, $M_0$ is the mass of the gas inside the bubble, $\gamma $ is its ratio of specific heats, $\rho $ is the density of the liquid and $p_\infty $ is the pressure far away from the bubble. The value $\omega _{\mathrm {eq}}$ is the Minnaert resonant frequency [28].

All lengths and times are made dimensionless by scaling them on the scales $R_0$ and $1/\omega _{\mathrm {eq}}$ , respectively, while the velocity potential $\Phi $ is rendered dimensionless by reference to the quantity $\omega _{\mathrm {eq}}R_0^2$ . (For further details of this scaling, see [7].) In dimensionless variables, the system is described by Laplace’s equation,

(2.1) $$ \begin{align} \nabla^2\Phi = \frac{1}{r^2}\dfrac{\partial {}}{\partial {r}}\bigg(r^2\dfrac{\partial {\Phi}}{\partial {r}}\bigg) + \frac{1}{r^2 \sin\phi}\dfrac{\partial {}}{\partial {\phi}}\bigg(\sin\phi\dfrac{\partial {\Phi}}{\partial {\phi}}\bigg) = 0, \end{align} $$

and two boundary conditions on $r = R(\phi ,t)$ . As the bubble boundary is a material surface, there is a kinematic boundary condition which demands that

(2.2) $$ \begin{align} \dfrac{\partial {\Phi}}{\partial {r}} = \dfrac{\partial {R}}{\partial {t}} + \frac{1}{R^2}\dfrac{\partial {\Phi}}{\partial {\phi}}\dfrac{\partial {R}}{\partial {\phi}} \quad \text{on } r = R(\phi,t). \end{align} $$

There is also a dynamic boundary condition which requires that the pressure within the bubble at the surface equals that in the liquid at the surface plus effects from the surface tension and Rayleigh viscosity. In dimensionless terms, this is given by

(2.3) $$ \begin{align} p = p_B - \kappa \kappa_S + 3\gamma\mu \Phi \quad \text{on } r = R(\phi,t), \end{align} $$

where $\kappa $ and $\mu $ are the dimensionless surface tension and Rayleigh viscosity coefficients, and $\kappa _S$ is the axisymmetric curvature given by

(2.4) $$ \begin{align} \kappa_S = \bigg[\frac{R^2 + 2R^2_\phi - RR_{\phi \phi}}{(R^2 + R^2_\phi)^{3/2}}\bigg]. \end{align} $$

As the gas in the bubble is isentropic, its pressure $p_B$ and density $\rho _B$ (in dimensional variables) are related by $p_B = K_B\rho _B^\gamma $ , where $\rho _B$ can be found by dividing the mass of gas by the volume. Hence, the dimensionless pressure in the bubble is given by

$$ \begin{align*} p_B = \frac{2^\gamma}{\big[\!\int^\pi_0 R^3 \sin\phi \,d\phi\big]^\gamma}. \end{align*} $$

The pressure in the liquid is supplied by a dimensionless unsteady Bernoulli equation

$$ \begin{align*} p = p_\infty(t) - 3\gamma\bigg[\dfrac{\partial {\Phi}}{\partial {t}} + \frac{1}{2}||\nabla{\Phi}||^2\bigg] \end{align*} $$

and so, after rearranging, the dynamic condition becomes

(2.5) $$ \begin{align} \dfrac{\partial {\Phi}}{\partial {t}} + \frac{1}{2}\bigg[\bigg(\dfrac{\partial {\Phi}}{\partial {r}}\bigg)^2 + \frac{1}{R^2}\bigg(\dfrac{\partial {\Phi}}{\partial {\phi}}\bigg)^2\bigg] + \mu \Phi ={}& \frac{p_\infty(t)}{3\gamma} + \frac{\kappa}{3\gamma} \bigg[\frac{R^2 + 2R^2_\phi - RR_{\phi \phi}}{(R^2 + R^2_\phi)^{3/2}}\bigg] \notag\\ &- \frac{2^\gamma}{3\gamma\big[\!\int^\pi_0 R^3 \sin\phi \,d\phi\big]^\gamma}. \end{align} $$

Our bubble system is thus fully described by the dimensionless Laplace equation (2.1), which is solved subject to the kinematic requirement (2.2) and the dynamic constraint (2.5). We consider two solutions to this problem. We begin with a linear model which will provide a closed-form solution to aid our conceptual understanding of the physical behaviour of the bubble. A fully nonlinear model solved with semi-numerical methods will then capture the nonlinear dynamics missing from the linear version and a comparison of the two models will be useful both for considering the accuracy of the models and for gauging the impact of nonlinear effects.

Before continuing to the models themselves, it is worth considering the circumstances where our theory is likely to hold relevance. For a spherical bubble, to derive the classic Rayleigh–Plesset equation, the dimensional pressure at the bubble surface is given by

(2.6) $$ \begin{align} p(R) = p_B -\frac{2\sigma}{R} - \frac{4\mu_L\dot{R}}{R}, \end{align} $$

where $\sigma $ is the surface tension and $\mu _L$ is the liquid viscosity [34]. Voinov and Golovin use Lagrange equations for a system of spherical bubbles in a liquid of small viscosity and find the same viscosity term [43]. Our Rayleigh viscosity behaves similarly, since (if the bubble is spherical) the $\mu \Phi $ term in our pressure equation (2.3) behaves rather like $-\mu R\dot {R}$ . By comparing (2.3) and (2.6) in dimensional terms, we find that dimensionless surface tension and Rayleigh viscosity coefficients are given by

$$ \begin{align*} \kappa = \frac{2\sigma}{p_\infty R_0} \quad \text{and} \quad \mu = \frac{4\mu_L}{\rho\omega_{\mathrm{eq}} R_0^2}. \end{align*} $$

Simulations were run with $\kappa = 0.01$ , and $\mu $ values generally ranging from 0.01 to 0.1. The $\kappa $ , $\mu $ and $R_0$ values for bubbles of air in several real liquids are given in Table 1. Bubbles of these 10–100 $\mu $ m sizes are consistent with cavitation bubbles in experimental results [38]. Increasing the surface tension coefficient $\kappa $ corresponds to reducing the size of the bubble as surface tension has a greater effect on smaller bubbles. Simulations lasting 200 dimensionless time units would, in real time, last less than a millisecond. Given the small radii of the bubbles and the short timescales involved, it is unlikely that the addition of gravity would have any significant effect.

Table 1 Properties of various liquids [15, 48] and associated $\kappa $ and $\mu $ values.

( $^a$ Silicone oil can have a large range of viscosities and thus allows for a range of valid $\mu $ values [48].)

2.1. Linear model

We suppose that an axisymmetric perturbation is made to the spherical equilibrium. We let

$$ \begin{align*} &R(\phi,t) = R_{\mathrm{eq}} + \varepsilon R_1(\phi,t) + \mathcal{O}(\varepsilon^2),\\ &\Phi(r,\phi,t) = \varepsilon \Phi_1(r,\phi,t) + \mathcal{O}(\varepsilon^2),\\ &p_\infty(t) = 1 + \varepsilon p_1(t) + \mathcal{O}(\varepsilon^2), \end{align*} $$

and consider an acoustic field given by $p_1(t) = p_a \sin (\tau t)$ .

At $\mathcal {O}(\varepsilon )$ , the linearized Laplace’s equation and kinematic boundary condition give

(2.7) $$\begin{align}\begin{aligned} &\nabla^2\Phi_1 = 0,\\ &\dfrac{\partial {\Phi_1}}{\partial {r}} = \dfrac{\partial {R_1}}{\partial {t}} \quad \text{on } r = R_{\mathrm{eq}}. \end{aligned} \end{align}$$

The dynamic boundary condition gives, at $\mathcal {O}(1)$ and at $\mathcal {O}(\varepsilon )$ respectively,

(2.8a) $$ \begin{align} R_{\mathrm{eq}}^{3\gamma} + \kappa R_{\mathrm{eq}}^{3\gamma -1} &= 1, \end{align} $$

(2.8b) $$ \begin{align} \dfrac{\partial {\Phi_1}}{\partial {t}} + \mu \Phi_1 &= \frac{p_a \sin(\tau t)}{3\gamma} + \frac{1}{2R_{\mathrm{eq}}^{3\gamma + 1}}\int^\pi_0 R_1 \sin\phi \,d\phi - \frac{\kappa (R_1 + R_{1,\phi \phi})}{3\gamma R_{\mathrm{eq}}^2} \quad \text{on } r = R_{\mathrm{eq}}. \end{align} $$

We can find the equilibrium radius from the $\mathcal {O}(1)$ dynamic condition (2.8a). This equation cannot be solved in closed form, but given that $\kappa \ll 1$ , we find an approximate solution by expanding in powers of $\kappa $ ,

$$ \begin{align*} R_{\mathrm{eq}} \approx 1 - \frac{\kappa}{3\gamma} + O(\kappa^2). \end{align*} $$

Hence, we see that as $\kappa $ increases, the equilibrium radius of the bubble is reduced. Of course, a more accurate result could be found by solving for $R_{\mathrm {eq}}$ numerically.

We are now ready to satisfy Laplace’s equation by taking a series of Legendre polynomials for $\Phi _1(r,\phi ,t)$ and $R_1(\phi ,t)$ :

$$ \begin{align*} R_1(\phi,t) &= \sum_{n=0}^N a_n(t)P_n(\cos\phi), \\ \Phi_1(r,\phi,t) &= \sum_{n=0}^N\frac{b_n(t)}{r^{n+1}}P_n(\cos\phi), \end{align*} $$

where the Fourier coefficients, $a_n(t)$ and $b_n(t)$ , are unknown functions of time which must be found from the two boundary conditions. We choose to truncate these would-be infinite sums at some arbitrary number N. To find the Fourier coefficients, the series forms are substituted into the boundary conditions, which are then multiplied by $P_l(\cos \phi ) \sin \phi $ and integrated from $0$ to $2\pi $ with respect to $\phi $ .

Applying standard orthogonality conditions to $\mathcal {O}(\varepsilon )$ of the kinematic condition (2.7) gives a relationship between the $a_n(t)$ and $b_n(t)$ coefficients so that

(2.9) $$ \begin{align} b_l(t) = - \frac{R_{\mathrm{eq}}^{l+2}a_l'(t)}{l+1}. \end{align} $$

The dynamic constraint (2.8b) leads to

(2.10) $$ \begin{align} &\sum_{n=0}^N\frac{b_n'(t) + \mu b_n(t)}{R_{\mathrm{eq}}^{n+1}}\int^\pi_0 P_n(\cos\phi)P_l(\cos\phi)\sin\phi\ d\phi \notag \\&\quad= \bigg[\frac{p_a \sin(\tau t)}{3\gamma} + \frac{a_0(t)}{R_{\mathrm{eq}}^{3\gamma + 1}}\bigg]\int^\pi_0 P_l(\cos\phi)\sin\phi\ d\phi \notag \\& \qquad - \frac{\kappa}{3\gamma R_{\mathrm{eq}}^2}\sum_{n=0}^N [(1-n(n+1)]a_n(t)\int^\pi_0 P_n(\cos\phi)P_l(\cos\phi)\sin\phi\ d\phi \notag \\&\qquad - \frac{\kappa}{3\gamma R_{\mathrm{eq}}^2}\sum_{n=0}^N a_n(t)\int^\pi_0 P'_n(\cos\phi)P_l(\cos\phi)\cos\phi\sin\phi\ d\phi. \end{align} $$

The first three integrals are well-known Legendre orthogonality integrals, but the final integral is given by

$$ \begin{align*} \int^\pi_0 P'_n(\cos\phi)P_l(\cos\phi)\cos\phi\sin\phi\ d\phi = \begin{cases} \dfrac{2l}{2l+1} & \text{if } n = l,\\ 2 & \text{if } n + l \text{ is even and } n> l, \\ 0 & \text{otherwise,} \end{cases} \end{align*} $$

so that the final sum term in (2.10) yields

$$ \begin{align*} & \sum_{n=0}^N a_n(t)\int^\pi_0 P'_n(\cos\phi)P_l(\cos\phi)\cos\phi\sin\phi\ d\phi \\ &\quad = \sum_{n=0}^{l-1} a_n(t)[0] + a_l(t)\bigg[\frac{2l}{2l+1}\bigg] + \sum_{i=1}^{(N-l)/2} a_{l+2i}(t)[2] + \sum_{i=0}^{(N-l-1)/2} a_{l+2i+1}(t)[0] \end{align*} $$

and, after rearranging, (2.10) becomes

$$ \begin{align*} b_l'(t) + \mu b_l(t) ={}& R_{\mathrm{eq}}\bigg[\frac{p_a \sin(\tau t)}{3\gamma} + \frac{a_0(t)}{R_{\mathrm{eq}}^{3\gamma + 1}}\bigg]\delta_{l0} \notag \\ & - \frac{\kappa R_{\mathrm{eq}}^{l-1}}{3\gamma}\bigg[(1-l^2)a_l(t) + (2l+1)\sum_{i=1}^{(N-l)/2}a_{l+2i}(t)\bigg], \end{align*} $$

where $\delta _{l0} = 1$ if $l=0$ , and is zero otherwise. Then, using relation (2.9),

(2.11) $$ \begin{align} a_l"(t) + \mu a_l'(t) + \kappa\alpha_l a_l(t) + \frac{a_0(t)\delta_{l0}}{R_{\mathrm{eq}}^{3\gamma + 2}} = -\frac{p_a \sin(\tau t)\delta_{l0}}{3\gamma} + \kappa \beta_l \sum_{i = 1}^{(N-l)/2} a_{l+ 2i}(t) , \end{align} $$

where

(2.12) $$ \begin{align} \alpha_l = \frac{(l+1)(l^2 - 1)}{(3\gamma R_{\mathrm{eq}}^3)} \quad \text{and}\quad \beta_l = \frac{(l+1)(2l+1)}{(3\gamma R_{\mathrm{eq}}^3)}. \end{align} $$

In a fully nonlinear model, all the modes are coupled enabling complex nonlinear interactions to occur. That is not the case in this linear model, although surface tension does allow for some limited interaction. Consider the system of equations (2.11). The final equation in this set (with $l=N$ ) is homogeneous, but $a_N(t)$ appears in the right-hand side of the equation with $l=N-2$ . Moreover, both $a_N(t)$ and $a_{N-2}(t)$ appear in the right-hand side of the equation for $l=N-4$ . This pattern continues and the equation with $l=0$ includes every even Fourier coefficient in its right-hand side, while that with $l=1$ includes every odd Fourier coefficient. The system of equations admits the solution

$$ \begin{align*} a_0(t) &= e^{-\mu t/2} \bigg[\hat{A}_0\cos\bigg(\frac{\omega_0t}{2}\bigg) + \hat{B}_0\sin\bigg(\frac{\omega_0t}{2}\bigg) \notag \\ & \quad + \sum_{i=1}^{N/2} \frac{\kappa \beta_0 T_{0,i}}{\bar{\alpha}_0 - \kappa\alpha_{2i}} \bigg\{\hat{A}_{2i}\cos\bigg(\frac{\omega_{2i}t}{2}\bigg) + \hat{B}_{2i}\sin\bigg(\frac{\omega_{2i}t}{2}\bigg) \bigg\}\bigg] \notag\\ &\quad + \frac{p_a[\mu\tau\cos(\tau t) - (\bar{\alpha}_0 - \tau^2)\sin(\tau t)]}{3\gamma[(\bar{\alpha}_0 -\tau^2)^2 + (\mu\tau)^2]}, \\ a_1(t) &= \hat{A}_1 + \hat{B}_1 e^{-\mu t} \notag \\ &\quad - e^{-\mu t/2} \bigg[\sum_{i=1}^{(N-1)/2} \frac{\beta_1 T_{1,i}}{\alpha_{1+2i}} \bigg\{\hat{A}_{1+2i}\cos\bigg(\frac{\omega_{1+2i}t}{2}\bigg) + \hat{B}_{1+2i}\sin\bigg(\frac{\omega_{1+2i}t}{2}\bigg) \bigg\}\bigg], \end{align*} $$

and, for $l = 2,\ldots ,N$ ,

$$ \begin{align*} a_l(t) &= e^{-\mu t/2} \bigg[\hat{A}_l\cos\bigg(\frac{\omega_lt}{2}\bigg) + \hat{B}_l\sin\bigg(\frac{\omega_lt}{2}\bigg) \notag\\ & \quad + \sum_{i=1}^{N/2} \frac{\beta_l T_{l,i}}{\alpha_l - \alpha_{l+2i}} \bigg\{\hat{A}_{l+2i}\cos\bigg(\frac{\omega_{l+2i}t}{2}\bigg) + \hat{B}_{l+2i}\sin\bigg(\frac{\omega_{l+2i}t}{2}\bigg) \bigg\}\bigg], \end{align*} $$

where

(2.13) $$\begin{align}\begin{aligned} T_{l,i} = \prod_{j=1}^{i-1}\bigg(1 + \frac{\beta_{l+2j}}{\alpha_{l+2j} - \alpha_{l+2i}}\bigg), \quad \bar{\alpha}_0 = \kappa\alpha_0 + \frac{1}{R_{\mathrm{eq}}^{3\gamma + 2}}\\ \omega_0 = \sqrt{4\bar{\alpha}_0 - \mu^2} \quad \text{and} \quad \omega_l = \sqrt{4\kappa \alpha_{l} - \mu^2}. \end{aligned} \end{align}$$

The remaining $b_n(t)$ Fourier coefficients can be found using relation (2.9) and are listed in Appendix A. If the initial conditions, $a_l(0)$ and $b_l(0)$ , are known, then the $\hat {A}$ and $\hat {B}$ coefficients can be found in the order $A_N$ to $A_0$ excluding $A_1$ , then $B_N$ to $B_0$ , and finally $A_1$ . These constants are also noted in Appendix A.

In general, each part of the solution is given by trigonometric functions multiplied by a negative exponential, so the bubble is stable as perturbations will oscillate with decaying amplitude. This is owing to the dissipative effect of the artificial Rayleigh viscosity. The homogeneous part of $a_1(t)$ , which is neither a volume nor a shape mode, but instead a translation up or down the z-axis, does not correspond to a decaying oscillation, but instead decays to a constant (because $\alpha _1 = 0$ ). The inclusion of surface tension means that not only will a perturbation to a certain mode cause oscillations in that mode, but a perturbation to an even mode l will produce oscillations in every even mode $n \leq l$ , and likewise oscillations in every odd mode $n \leq l$ if the perturbation is odd.

As $t \rightarrow \infty $ , the solution tends to

(2.14) $$ \begin{align} &a_0(t) = \frac{p_a[\mu\tau\cos(\tau t) - (\bar{\alpha}_0 - \tau^2)\sin(\tau t)]}{3\gamma[(\bar{\alpha}_0 -\tau^2)^2 + (\mu\tau)^2]},\nonumber\\ &a_1(t) = \hat{A}_1,\\ &a_l(t) = 0.\nonumber \end{align} $$

Here, only the spherical forced response and a translation along the z-axis remain, since all nonspherical behaviour has been dissipated by the Rayleigh viscosity. Hence, it can be seen that the bubble is highly stable.

Our solution can be contrasted with an excessively simplified model which excludes surface tension, Rayleigh viscosity and the acoustic field. The solution for this model is just $a_0(t) = A_0\cos t + B_0\sin t$ and $a_n(t) = A_n + B_n t$ . In this case, spherical perturbations away from the equilibrium radius are neutrally stable as they cause spherical oscillations of constant amplitude, but nonspherical perturbations to the velocity potential (which cause $B_n$ to be nonzero) are unstable as they grow linearly with time, ultimately violating the underlying linear assumption. Rayleigh viscosity and surface tension clearly provide important stabilizing effects. At the resonant frequency, $\tau = \bar {\alpha }_0$ , (2.14) shows that $a_0(t) = p_a\cos (\bar {\alpha }_0 t)/(3\gamma \mu \bar {\alpha }_0)$ in the large time limit. This is a resonance peak with a maximum value which increases as the Rayleigh viscosity decreases. If there were no Rayleigh viscosity, $\mu = 0$ , there would be oscillations of infinite amplitude.

Herein lies a major attraction of the linear model: as the solution is written in closed form, the impact of physical effects, such as surface tension, Rayleigh viscosity and the acoustic forcing field, can be conceptually understood, and the behaviour at large time can be discerned at a glance. Surface tension causes nonspherical disturbances to oscillate rather than grow linearly with time and enables interactions between modes. The Rayleigh viscosity ensures that nonspherical disturbances dissipate while a forcing field causes enduring spherical oscillations.

2.2. Nonlinear model

Now that we have a linear model which provides a fundamental understanding of the physical behaviour of the bubble, we can turn our attention to a nonlinear model which will capture important effects necessarily absent from our linear model. We again begin by satisfying Laplace’s equation in the liquid surrounding the bubble with a Fourier–Legendre series. We then use semi-analytical methods to impose the full nonlinear boundary conditions. Suppose that the radius of the bubble and the velocity potential outside the bubble can be expressed as

(2.15) $$\begin{align} \begin{aligned} R(\phi,t) &= \sum_{n = 0}^N \bar{A}_n(t) P_n(\cos\phi),\\ \Phi(r,\phi,t) &=\sum_{n = 0}^N \frac{\bar{B}_n(t) P_n(\cos\phi)}{r^{n+1}} \quad \text{for } r> R(\phi,t). \end{aligned} \end{align}$$

These representations become increasingly accurate as N grows. Truncating the series forms in the linear model generally has no effect, as the sum naturally truncates itself at the highest initially perturbed mode. Surface tension allows perturbed modes to activate lower modes, but there is no mechanism for driving higher modes. Within the nonlinear model, this is no longer the case; higher modes can be excited and truncating the sum can have a significant effect. Simulations were run with $N = 81$ , and the validity of this choice is discussed in Section 4.

Equations for the functions $\bar {A}_n(t)$ and $\bar {B}_n(t)$ are derived by substituting the series solutions into the kinematic boundary condition (2.2) and the dynamic boundary condition (2.5). Using well-known orthogonality properties of Legendre polynomials, the boundary conditions imply that

(2.16) $$ \begin{align} \bar{A}_l'(t) = -\sum_{n = 0}^N \bar{B}_n(t)[I^{(1)}_{ln} + I^{(2)}_{ln}] \end{align} $$

and

(2.17) $$ \begin{align} \sum_{n = 0}^N \bar{B}_n'(t) I^{(3)}_{ln} = - \mu\sum_{n = 0}^N \bar{B}_n(t) I^{(3)}_{ln} - \frac{1}{2}I^{(4)}_{l} + \frac{2}{3\gamma}\bigg[p_\infty(t) - \frac{2^\gamma}{{I^{(5)}}^\gamma}\bigg] \delta_{l0} + \frac{\kappa}{3\gamma} I^{(6)}_{l} \end{align} $$

for $0\leq l\leq N$ . The system (2.16)–(2.17) constitutes a set of $2N+2$ differential equations for the evolution of the time-dependent Fourier-series coefficients, $A_n(t)$ and $B_n(t)$ . It was solved numerically with MATLAB’s ode45 which uses an explicit fourth/fifth-order Runge–Kutta method. Within these equations, the quantities $I^{(1)}_{ln}$ to $I_l^{(6)}$ are integrals which are listed in Appendix B, while $\delta _{00} = 1$ with $\delta _{l0} = 0$ for all other l. The various integrals were evaluated using Gauss–Legendre quadrature with $5N$ points whose weights and nodes were determined using an algorithm developed by von Winckel [44].

3. Extended results for spherical bubbles

We can expect the nonlinear nonspherical model to be unstable for high forcing pressure amplitudes and small Rayleigh viscosities, as well as being slow to run numerically. The linear model has shown that as time advances, dissipation from the Rayleigh viscosity causes nonspherical bubbles to become spherical. To extend our results, and to investigate high forcing pressure amplitudes and low Rayleigh viscosities in particular, we consider a purely spherical model.

When considering a purely spherical bubble, the system given in Section 2 is simplified considerably. We now require all $\phi $ derivatives to be zero, so $ {\partial /\partial \phi \equiv \partial /\partial \theta \equiv 0}$ . The radius, $R(t)$ , is now a function of t alone, and the velocity potential $\Phi (r,t)$ satisfies

$$ \begin{align*} \nabla^2\Phi = \frac{1}{r^2}\dfrac{\partial {}}{\partial {r}}\bigg(r^2\dfrac{\partial {\Phi}}{\partial {r}}\bigg) = 0. \end{align*} $$

The boundary conditions now apply on $r = R(t)$ ; the kinematic boundary condition is simply

$$ \begin{align*} \dfrac{\partial {\Phi}}{\partial {r}} = \dfrac{\partial {R}}{\partial {t}} \end{align*} $$

and the dynamic boundary condition is

$$ \begin{align*} \dfrac{\partial {\Phi}}{\partial {t}} + \frac{1}{2}\bigg(\dfrac{\partial {\Phi}}{\partial {r}}\bigg)^2 + \mu \Phi = \frac{p_\infty(t)}{3\gamma} - \frac{1}{3\gamma R(t)^{3\gamma}}, \end{align*} $$

where again, $p_\infty (t)$ is the dimensionless far field pressure and $\mu $ is an artificial Rayleigh viscosity. In the nonspherical model, surface tension ensures that initial disturbances cause oscillations rather than linear growth and causes initial disturbances to excite other lower modes of oscillation. In the spherical model, surface tension is not able or required to do these things and so, for simplicity, it has been excluded from this model. This means that the equilibrium radius is given by $R_{\mathrm {eq}} = 1$ .

A spherical linear model can be created using the same process as for the nonspherical linear model in Section 2.1. This spherical linear model is a special case of the nonspherical version and is given by

$$ \begin{align*} R(t) &= 1 + \varepsilon e^{-\mu t/2} \bigg[A\cos\bigg(\frac{t}{2}\sqrt{4 - \mu^2}\bigg) + B\sin\bigg(\frac{t}{2}\sqrt{4 - \mu^2}\bigg)\bigg] \notag\\ &\quad + \frac{\varepsilon p_a[\mu\tau\cos(\tau t) - (1 - \tau^2)\sin(\tau t)]}{3\gamma[(1 -\tau^2)^2 + (\mu\tau)^2]} + \mathcal{O}(\varepsilon^2). \end{align*} $$

As $t \rightarrow \infty $ , the first oscillatory part of this solution decays to zero, leaving only the oscillations at the forcing frequency $\tau $ . Once again, the Rayleigh viscosity prevents divergence at the natural frequency which is now at $\tau = 1$ . At this natural frequency, $R(t) = 1 + \varepsilon p_a\cos (t)/(3\gamma \mu )$ in the large time limit. So for $R(t)> 0$ , we require $\varepsilon p_a < 3\gamma \mu $ . If the forcing pressure amplitude increases beyond this, then this linear assumption breaks down as it predicts unphysical negative bubble radii.

3.1. Nonlinear spherical model

As for the nonspherical case, we can build a numerical model to implement the full nonlinear boundary conditions. We satisfy Laplace’s equation by taking

$$ \begin{align*} \Phi(r,t) = \frac{\bar{V}(t)}{r} \end{align*} $$

to be the velocity potential while the radius is simply given by $R(t)$ . The kinematic condition then gives

(3.1) $$ \begin{align} R'(t) = -\frac{\bar{V}(t)}{R(t)^2}, \end{align} $$

and the dynamic condition gives

(3.2) $$ \begin{align} \bar{V}'(t) = - \frac{\bar{V}(t)^2}{2R(t)^3} - \mu \bar{V}(t) + \frac{R(t)}{3\gamma}[ 1 + \varepsilon p_a \sin(\tau t) - R(t)^{-3\gamma}]. \end{align} $$

The system (3.1)–(3.2) was again solved numerically with MATLAB’s ode45 to find the evolution of the radius $R(t)$ and the velocity potential $\Phi (r,t) = \bar {V}(t)/r$ . As well as only being a system of two ordinary differential equations (ODEs), this requires no evaluation of integrals and no matrix inversion, and is thus extremely fast. Further, without nonspherical effects, the system is highly stable and curvature singularities do not occur. The speed and stability of this model enabled large numbers of simulations to be run (including in parallel) with varying initial conditions and parameter values ( $\mu $ , $p_a$ and $\tau $ ) including values which would cause the simulations to fail in the nonspherical model.

4. Results

4.1. Nonspherical bubbles

We begin this presentation of results by considering bubbles that are initially spherical, but which are perturbed by a small disturbance to the velocity using a single Legendre mode. Multi-modal initial velocity disturbances and initial radial disturbances are also possible in both the linear and nonlinear models, but for simplicity, we consider only single-mode disturbances here. For all nonspherical simulations presented, the Fourier sum truncates at $N = 81$ .

Before commenting further on the results, attention should be drawn to the convergence of the Fourier series in the nonlinear model. Calculating the curvature along the interface is a very sensitive test of the convergence of the Fourier series and the accuracy of the results. This is because the curvature in (2.4) involves the second derivative of the interface profile. Figure 1 presents the curvature $\kappa _S$ computed with total numbers $N=21$ , $61$ and $81$ of Fourier modes, for a bubble with a mode-5 initial disturbance. (The same bubble will be considered again in more detail in Figure 3.) Although there are modest errors in the curvature profile obtained with $N=21$ modes, the two sets of results for $N=61$ and $N=81$ are in very close agreement; the inset on the right side of Figure 1 shows that convergence for the curvature has been achieved to at least three significant figures. This confirms the high accuracy of our results and, indeed, the convergence of the series (2.15) for the interface $r=R(\phi ,t)$ is substantially better than this as no second derivatives are required.

In general, early in a simulation, a bubble has a strongly nonspherical shape due to the initial nonspherical velocity disturbance. It exhibits volume oscillations caused by the initial perturbation (occurring at the natural frequency and at other frequencies from higher mode interactions caused by surface tension) and has forced oscillations at the frequency $\tau $ . Then, as time advances, the artificial Rayleigh viscosity dissipates all but the spherical response to the forcing frequency; a larger viscosity, $\mu $ , means that this dissipation occurs faster. However, if the bubble is forced at the natural frequency, then resonance occurs, and the bubble becomes highly unstable and the nonlinear simulation eventually fails due to the formation of curvature singularities on the bubble surface. In this case, the linear model gives physically unrealistic results since it predicts very large oscillations with amplitude proportional to $1/\mu $ .

Figure 2 shows the evolution with time of the Fourier coefficients from both the linear and nonlinear nonspherical models for a bubble with an initial mode-4 disturbance and a forcing frequency of $\tau = 2$ . The figure shows excellent agreement between the two models. It can be clearly seen that the nonspherical behaviour dissipates with time and that the natural spherical oscillations rapidly decay, leaving only the oscillations at the forcing frequency. This bubble was initially perturbed in the fourth mode, so the two models show that modes zero and two have also been excited. The nonlinear model also shows a small disturbance to the sixth mode (and increasingly smaller ones for even higher even modes). This is a result of the nonlinear coupling which is absent in the linear model. This activation of higher modes is why the nonlinear model can be prone to curvature singularities (in simulations with low, or no Rayleigh viscosity) which do not occur in the linear model.

Figure 1 The curvature of the bubble interface from three simulations with $N=21$ , 61 and 81 Fourier modes, respectively. In each simulation, the bubble had an initial disturbance of $b_5(0) = -0.1$ , a forcing of $\varepsilon p_1(t) = 0.1\sin (1.5t)$ , a Rayleigh viscosity of $\mu = 0.1$ and was allowed to evolve to $t=100$ .

Figure 2 Fourier coefficients from the linear (dashed) and nonlinear (solid) models of a bubble with initial disturbance $b_4(0) = -0.1$ , Rayleigh viscosity $\mu = 0.05$ and a forcing of $\varepsilon p_1(t) = 0.1\sin (2t)$ .

In Figure 3, the bubble is subjected to an initial mode-5 disturbance and is forced at $\tau = 1.5$ . As the initial disturbance was odd, there is an enduring mode-1 vertical translation. Furthermore, while the linear model predicts that the mode-5 disturbance causes further disturbances to modes-3 and 1, the nonlinear model predicts that even modes will be excited along with the odd modes. The fact that even disturbances excite even modes only, but odd disturbances are capable of generating both even and odd modes is widely established. It is inherently nonlinear behaviour and has been found in both experimental and theoretical work (see, for example, Shaw [40] and Guédra and Inserra [13]). Notably, the sixth mode (which is predicted to be zero by the linear model) displays resonant behaviour in the nonlinear model. In spherical models, while the second subharmonic (when the forcing frequency is around twice the natural frequency) is the most important subharmonic, the largest resonance is caused by the natural frequency. The natural frequency is thus the critical frequency. However, for nonspherical models, Hsieh [14] has shown that the critical frequency for a shape mode occurs at the second subharmonic. The natural frequency of each nonspherical mode is given by $\omega _l/2$ (see (2.12) and (2.13)). For $\kappa = 0.01$ and $\mu = 0.1$ , the natural frequency of the sixth mode is $\omega _6/2 = 0.766$ . So when the bubble is forced at approximately twice that value, at $\tau = 1.5$ , the second subharmonic resonance occurs in the sixth mode. Other second subharmonics can also be achieved, such as a mode-4 resonance with $\tau = 2(\omega _4/2) = 0.842$ or a mode-8 resonance with $\tau = 2(\omega _8/2) = 2.332$ . These nonspherical resonances ultimately grow so large as to cause the simulations to fail. Despite the occurrence of the resonance in the nonlinear model, the two models still show good agreement, especially in the zeroth mode. Note that as the natural frequency of oscillations of shape modes is independent of the azimuthal mode number, axisymmetric models are more widely applicable than might be expected [17].

Figure 3 Fourier coefficients from the linear (dashed) and nonlinear (solid) models of a bubble with initial disturbance $b_5(0) = -0.1$ , Rayleigh viscosity $\mu = 0.1$ and forcing of $\varepsilon p_1(t) = 0.1\sin (1.5t)$ .

For Figure 4, the nonspherical nonlinear model was run from $t = 0$ to $t = 200$ , for increasing forcing pressures or frequencies. Over time, the Rayleigh viscosity causes bubbles to become increasingly spherical and so, in general, if a simulation reached $t = 200$ , the bubble would be completely spherical at that time. However, for large forcing pressures and some frequencies, the simulation failed before $t = 200$ . (The simulation was deemed to have failed if $|b_n(t)|> 0.15$ , as this occurs when the bubble is experiencing physically unrealistic, rapidly changing, high-mode surface oscillations which are indicative of curvature singularities.) Figure 4(a) shows the failure time against forcing pressure $\varepsilon p_a$ for two different forcing frequencies for bubbles with an initial mode-4 disturbance. For $\tau = 2$ , failure before $t=200$ began at the forcing pressure $\varepsilon p_a = 0.34$ . As $\tau = 0.9$ is close to the spherical resonant frequency, bubbles forced at this frequency are far less stable and fail before $t=200$ for forcing amplitudes greater than $\varepsilon p_a = 0.045$ .

Figure 4 The failure time of the nonlinear nonspherical bubble model for a range of forcing frequencies $\tau $ and forcing pressure amplitudes $\varepsilon p_a$ .

In Figure 4(b), the failure time is given over a range of forcing frequencies for a bubble with an initial mode-6 disturbance and a forcing amplitude of $\varepsilon p_a = 0.2$ . The bubble simulation fails before $t=200$ for $\tau $ values between 0.7 and 1.7. The fourth mode has a natural frequency $\omega _4/2 = 0.421$ , and its second subharmonic resonance causes the simulation to fail for $0.71 < \tau < 0.83$ , and likewise its third subharmonic causes failure for $1.15 < \tau < 1.28$ . Similarly, the mode-6 second subharmonic resonance causes failure for $1.29 < \tau <1.75$ . Near the natural frequency of the spherical mode ( $\tau = 1$ ), growth in the spherical oscillations causes the instability, although the fourth mode remains significant. At half the natural frequency, $\tau = 0.5$ , the simulation also fails early since nonlinear mode coupling excites the primary resonance at $\tau =1$ .

4.2. Spherical bubbles

The spherical simulations are not subject to curvature singularities and so are far more stable than their nonspherical counterparts. This enables simulations to be run with higher forcing pressures, lower forcing frequencies and lower Rayleigh viscosities. Of particular interest are regions where the nonlinear model predicts the existence of multiple solutions caused by subharmonic and superharmonic structures. The spherical simulations are extremely fast. Running in parallel, the 200 nonspherical simulations for Figure 4(b) took approximately 12 hours (50 simulations each for four cores on a standard laptop), whereas one core can run 80,000 spherical simulations in two hours.

Analysis is undertaken by comparing the predicted amplitudes for the two models. The large time limit of the linear model has the form

$$ \begin{align*} R(t) = 1 + \varepsilon\frac{p_a[\mu\tau\cos(\tau t) - (1 - \tau^2)\sin(\tau t)]}{3\gamma[(1 -\tau^2)^2 + (\mu\tau)^2]}. \end{align*} $$

From this form, the (peak-to-trough) amplitude can be calculated as

$$ \begin{align*} \text{Amp} = \frac{2\varepsilon p_a}{3\gamma\sqrt{(\mu\tau)^2 + (1 - \tau^2)^2}}. \end{align*} $$

To find the amplitude from the fully nonlinear model, first, the simulation was allowed to run until it had converged to a stable solution. The bubble was considered to have converged if the $R(t)$ and $\bar {V}(t)$ values at the end of each period fell within a margin of $0.3\%$ for 10 consecutive periods. Fundamental frequency solutions and superharmonics were found using the forcing period. If this was unsuccessful, an integer multiple of the forcing period was used, which enabled subharmonics and ultraharmonics to be found. When the solution had converged, the period was noted, as was the number of maxima per period and their values. The absolute amplitude was then calculated from the last period by subtracting the minimum $R(t)$ from the maximum $R(t)$ .

Figure 5 The linear (dashed) and nonlinear (solid) predictions of the amplitude of bubble oscillations against forcing frequency $\tau $ for two different forcing amplitudes from simulations with $\mu = 0.05$ . The legend gives the period as a multiple of the forcing period and the number of distinct maxima for the nonlinear solution.

For Figures 5–9, a series of simulations was run with varying values of the forcing frequency ( $\tau $ ), the forcing amplitude ( $\varepsilon p_a$ ) or the Rayleigh viscosity ( $\mu $ ). The two initial conditions are given by the radius $R(0)$ and the velocity potential $\Phi (r,0) = \bar {V}(0)/r$ . The first simulation started with a low forcing pressure or a forcing frequency away from any resonances to ensure that the nonlinear model would give a very similar result to the linear model when started with the same initial conditions. The simulation was then rerun several times with varying initial conditions to find other possible solutions. The parameter of interest was then iterated, and the process was restarted using the final $R(t)$ and $\bar {V}(t)$ values of each converged solution from the previous iterative step as initial conditions, as well as using the linear guess and other variations. This process allowed for multiple solutions to be found, including subharmonics and superharmonics. Note that for these spherical models, surface tension is not included.

Figure 6 The nonlinear predictions of the amplitude of bubble oscillations against forcing frequency $\tau $ for two different forcing pressure amplitudes from simulations with $\mu = 0.01$ . The legend gives the period as a multiple of the forcing period and the number of distinct maxima. The horizontal axis is on a logarithmic scale.

Figure 7 The linear (dashed) and nonlinear (solid) predictions of the amplitude of bubble oscillations against the Rayleigh viscosity parameter $\mu $ for various forcing frequencies for simulations with $\varepsilon p_a = 0.3$ .

Figure 8 The nonlinear predictions for the maxima of the stable bubble oscillations against the forcing pressure. The legend gives the period as a multiple of the forcing period and the number of distinct maxima.

Figure 9 The nonlinear predictions for the maxima of the stable bubble oscillations against the forcing pressure from simulations with $\tau = 0.25$ and $\mu = 0.05$ . The lower plot shows the $0.8 < \varepsilon p_a < 1$ region in greater detail. The legend groups the period (as a multiple of the forcing period) by period-doubling cascades.

Figure 5 shows resonance plots of the amplitude of oscillations against the forcing frequency for the two forcing pressures $\varepsilon p_a = 0.3$ and $\varepsilon p_a = 0.03$ . When $\varepsilon p_a = 0.03$ , the two models give very similar resonance spikes at the natural frequency $\tau = 1$ . The nonlinear model also predicts a second, far smaller, superharmonic resonance spike at $\tau = 1/2$ . When $\varepsilon p_a = 0.3$ , the form of the linear solution is unchanged, with the height of the resonance spike merely increasing by the expected factor of 10, but the nonlinear resonance spike has been deflected to the left. The left side of the nonlinear resonance spike is missing. Presumably, this part of the solution forms the unstable portion of a hysteresis region, and the nonlinear model cannot converge to it, but instead converges to the higher-amplitude solution on the right of the resonance spike, or lower-amplitude solution on the linear curve. For this larger $\varepsilon p_a$ value, the superharmonic resonance at $\tau = 1/2$ is larger, and superharmonic resonances also occur at $\tau = 1/3$ and $1/4$ , as well as a subharmonic resonance at $\tau = 2$ .

The high resonance spikes from the linear model can lead to unrealistic negative $R(t)$ values. The linear model predicts smooth sinusoidal oscillations, so if the peak-to-trough amplitude increases above 2 (as it does in Figure 5(b) near $\tau = 1$ ), then the radius becomes negative. Importantly, the nonlinear model never predicts negative radii. As the minima approach zero, they become spiked instead of rounded. This is indicative of the extremely sudden bubble collapse that occurs in cavitation bubbles.

Figure 6 shows two more plots of amplitude against forcing frequency, both with small Rayleigh viscosity $\mu = 0.01$ . In Figure 6(a), the forcing pressure amplitude $\varepsilon p_a = 0.2$ causes a large left-leaning resonance spike at $\tau = 1$ , superharmonics at $\tau = 1/4$ , $1/3$ and $1/2$ , and subharmonics at $\tau = 2$ and 3. For the larger forcing pressure amplitude $\varepsilon p_a = 0.5$ in Figure 6(b), more superharmonic resonances are visible, and they decrease in maximum amplitude as $\tau $ is distanced from the resonance at $\tau = 1$ . As $\tau $ increases, the resonant spikes are increasingly bent to the left. For $\tau \gtrsim 1/5$ , this causes the unstable left sides of the resonance spikes to be absent. For lower $\tau $ values, the resonant spikes are essentially upright and so very little is missing.

For the natural resonant spike at $\tau = 1$ and low-amplitude solutions that lie on the linear solution curve, the bubble oscillates with the forcing period and has one maximum per period. For the superharmonic resonances around forcing frequencies $\tau _n = 1/n$ , the bubble still oscillates with the forcing period, but there are n maxima per period. For the subharmonics, the bubble oscillates with an integer multiple of the forcing period, but with only one maximum. Hence, for the superharmonic resonance at $\tau = 1/2$ , there are two maxima per forcing period, and for the $\tau = 2$ subharmonic, the bubble has a period of $2(2\pi /\tau )$ with one maximum. For $\varepsilon p_a = 0.5$ , an ultraharmonic is present at $\tau = 2/3$ , and the bubble oscillates with a period twice the forcing period and has three maxima per period. Resonance curves qualitatively similar to those in Figure 6 are widely available in papers such as those of Lauterborn and Kurz [24] and Gong et al. [12].

It is also possible to plot the curves of the amplitude against the Rayleigh viscosity. This is done for various forcing frequencies in Figure 7. In general, away from resonances, such as at $\tau = 0.6$ or 0.7, the linear and nonlinear predictions are very similar and the amplitude does not change significantly with $\mu $ . At $\tau = 0.5$ , the nonlinear model reacts to the second superharmonic and predicts high amplitude solutions with two maxima per forcing period for low Rayleigh viscosities (plotted in light blue). As the Rayleigh viscosity increases, it dampens the resonance, and at $\tau \approx 0.3$ , the nonlinear solution returns to having one maximum per period (plotted in dark blue) like the linear model which does not detect harmonic resonances. Both models are sensitive to the natural resonance at $\tau = 1$ , however, for low Rayleigh viscosities, at $\tau = 0.9$ , the linear prediction is lower than the nonlinear prediction, and at $\tau = 1.1$ , the reverse is true. This behaviour is caused by the nonlinear resonance spike curving to the left while the linear spike remains straight. This can be seen in Figure 5(b). As the Rayleigh viscosity $\mu $ increases, the strength of the resonance fades, and the linear and nonlinear predictions become extremely similar. In all other simulations, a small Rayleigh viscosity was chosen with $\mu $ set to 0.1, 0.05 or 0.01.

In a conventional bifurcation plot, the bubble radius $R(t)$ is noted at the end of each forcing period. Subharmonics, which oscillate with a multiple of the forcing period, are easily identified with this method; however, superharmonics, which have multiple maxima but oscillate with the forcing frequency, and ultraharmonics are not. By instead plotting the maxima of the stable oscillations, all harmonic resonances can be identified. Further discussion of this point and a comparison of the methods can be found from Sojahrood et al. [41]. In Figure 8, two such plots are given for different $\tau $ and $\mu $ values. Only one stable solution is present for the simulations with $\tau = 0.3$ and $\mu = 0.1$ . It starts with a period one oscillation, initially with one maximum but then with two and three maxima from the $\tau = 1/2$ and 1/3 superharmonics. The period then doubles and doubles again. The first solution for the simulations with $\tau = 0.5$ and $\mu = 0.05$ behaves similarly, although the period-4 maxima begin to converge after the initial divergence. This is evidence of so-called “period-bubbling” in which a closed structure is formed in the bifurcation diagram when solutions undergo period halving after period doubling [3]. Similar behaviours have been reported by Sojahrood et al. [41] and Lauterborn and Kurz [24]. A second high-amplitude, period-one solution appears at $\varepsilon p_a = 0.335$ . The existence of a second high-amplitude solution is common and occurs earlier for lower Rayleigh viscosities. For $\tau = 0.6$ , it appears at $\varepsilon p_a \approx 0.5$ and 0.25 for $\mu = 0.1$ and 0.05, respectively. A time-dependent plot of the radius of a very similar solution can be seen in Figure 12. There is also a short-lived third solution which begins with period-2 and three maxima (the $\tau = 2/3$ ultraharmonic), and doubles twice in quick succession to a period-8 solution with 12 maxima.

Figure 9 shows the period-doubling route to chaos for simulations with $\tau = 0.25$ and $\mu = 0.05$ . Until $\varepsilon p_a = 0.8$ , the results are broadly similar to those in Figure 8(a) which has a similar forcing frequency and twice the Rayleigh viscosity. Before the period doubling begins, the lower solution has four maxima from the $\tau = 1/4$ superharmonic, and the higher solution which begins at $\varepsilon p_a \approx 0.75$ has three maxima. When the period doubles, the number of maxima likewise doubles (although the number of maxima is only noted in the legend for the P1 solutions). For the lower solution, a period-doubling cascade begins at $\varepsilon p_a = 0.815$ , and doublings occur with increasingly close intervals at $\varepsilon p_a = 0.857$ , 0.867 and 0.870. The convergence check only allows for solutions up to period-20, so the final solution found before chaos is period-16. Beyond that, the region is generally blank, except for some period-3 solutions which double to P6, then P12, and period-5 solutions which double to P10. Period-3 solutions are indicative of chaos. At $\varepsilon p_a = 0.957$ , a period-8 solution reappears and the period successively halves to period-1 at $\varepsilon p_a = 0.987$ forming a Feigenbaum remerging tree [3]. This behaviour of period-doubling routes to chaos, accompanied by superharmonics and ultraharmonics, then followed by period halving has been found experimentally by Lauterborn and Cramer [23].

Figure 10 Contour plots of the maximum amplitude of oscillations for varying initial conditions at four different forcing frequencies for a bubble with $\mu = 0.05$ and $\varepsilon p_a = 0.5$ .

Figure 11 On the left is a time history of the three solutions that exist for the bubble with forcing frequency $\tau = 0.42$ , $\varepsilon p_a = 0.5$ and $\mu = 0.05$ from Figure 10. On the right is a phase plane diagram of the three solutions.

Figure 12 On the left is a time history of the two solutions that exist for the bubble with forcing frequency $\tau = 0.62$ , $\varepsilon p_a = 0.5$ and $\mu = 0.05$ from Figure 10. On the right is a phase plane diagram of the two solutions.

The spherical results thus far have examined the effect of varying the parameters $\mu $ , $\tau $ and $\varepsilon p_a$ . Figures 10–12 relate instead to varying the initial values, $R(0)$ and $\Phi (r,0) = \bar {V}(0)/r$ , of the bubble in the nonlinear model. Simulations were run for a grid of initial conditions with $0.1 \leq R(0) \leq 5$ and $-5 \leq \bar {V}(0) \leq 5$ with increments of 0.025. Each simulation was run until the bubble oscillations had converged to a stable solution, and then the amplitude was calculated. This was done for a range of forcing frequencies. Almost 80,000 simulations were run per forcing frequency, with each set taking approximately two hours of computing time to complete. This was performed on a standard four-core laptop running four frequencies in parallel.

For simulations with Rayleigh viscosity $\mu = 0.05$ and forcing pressure amplitude $\varepsilon p_a = 0.5$ , the nonlinear model predicts the existence of multiple solutions for the forcing frequencies $\tau = 0.34$ to $0.796$ . A contour plot of the amplitudes from the nonlinear model at four frequencies is given in Figure 10. The yellow (then mustard, green and aqua) colour represents a region of a high-amplitude solution with a single bubble collapse per forcing period. As $\tau $ increases, this region expands and spirals in towards the origin until, at $\tau = 0.79$ , it covers all but a thin spiralling line. Beyond $\tau = 0.8$ , it is the only remaining solution. Its amplitude reduces from approximately 4 at $\tau = 0.42$ to approximately 1.5 at $\tau = 0.8$ . The blue represents a region of lower-amplitude solutions that have the same period, but generally have two maxima per period due to the $\tau = 1/2$ superharmonic.

The case with forcing frequency $\tau = 0.42$ is examined in more detail in Figure 11, which shows both a plot of $R(t)$ against t and a phase plane diagram of $\bar {V}(t)$ and $R(t)$ . The highest-amplitude solution, plotted in yellow, has cavitation bubble-like collapses where $R(t)$ is extremely close to zero. In the phase plane, this solution looks like a single large loop. The mid-amplitude solution (light blue) has two maxima of similar values and two nearly overlapping loops in the phase plane. This solution has two collapses per forcing period. The lowest-amplitude smooth solution (dark blue) does not exhibit bubble collapse, as the minima remain smooth rather than pointed and do not approach zero. This solution has one large maximum followed by a smaller one and can be seen in the phase-plane as a small loop inside a larger loop. For $0.34 < \tau < 0.42$ , the smooth solution (dark blue) is the only low-amplitude solution present. Then, for $0.42 \leq \tau \leq 0.44$ , both low-amplitude solutions are present, with the mid-amplitude solution with collapses (light blue) appearing as a growing swirl inside the low-amplitude region. For $\tau> 0.44$ , the light-blue solution with its double collapse is the only remaining low-amplitude solution. As $\tau $ increases and becomes distanced from the $\tau = 1/2$ superharmonic, the strength of the second collapse is reduced and the first minimum becomes less pointed, until, at $\tau = 0.62$ , as seen in Figure 12, the second minimum can be seen only as a small dip between two similarly valued maxima. In the phase plane, the second minimum remains as a small kink on the right of the limit cycle. Beyond $\tau = 0.7$ , the second minimum has ceased to exist.

5. Conclusion

We have developed linear and nonlinear models with the aim of improving our understanding of the forced oscillations of spherical and nonspherical bubbles. The substantial benefit of the simpler linear model was that due to its closed-form solution, it allowed the physical effects of surface tension, Rayleigh viscosity and the acoustic forcing field to be clearly understood, and the behaviour at large time to be determined without costly simulations. In general, the Rayleigh viscosity caused nonspherical disturbances and oscillations at frequencies other than the forcing frequency to decay, and surface tension caused nonspherical disturbances to oscillate. It also provided limited mode coupling in the linear model to mimic some of the full nonlinear coupling in the full model. Hence, the models had excellent agreement and both predicted complex nonspherical shapes. The more realistic nonlinear model captured shape-mode resonances, which were the only enduring nonspherical behaviour. The nonlinear model sometimes experienced eventual failure when nonspherical disturbances grew too large. This tended to occur when the bubble was forced at these resonant forcing frequencies, or with high forcing pressures and low Rayleigh viscosities. Our numerical results suggest that failure is caused by curvature singularities that develop within a finite time on the bubble interface. Eliminating such singular behaviour is possible when fluid viscosity and finite-width interfacial regions are introduced into the nonlinear model (see [7]).

The inherently more stable spherical model allowed for consideration of these less stable parameter ranges. The nonlinear model exhibited complicated resonance behaviour including subharmonics, superharmonics and ultraharmonics, each with curved resonance spikes. Period-doubling and halving cascades were also evident. These nonlinear behaviours were stronger for high forcing pressures and low Rayleigh viscosities. Resonant behaviour in the linear model was limited to a straight resonance spike at the natural frequency.

A full description of the resonance behaviour within a general nonspherical bubble and the development of a weakly nonlinear model to bridge the gap between the linear and nonlinear spherical models remain topics for further study.