A tractable mathematical model for rectified diffusion

Abstract

A core unresolved challenge in rectified diffusion is the lack of predictive models, which allow accurate modelling of the growth of bubbles over millions of cycles of oscillation. The spherically symmetric problem involves two asymptotic regions in space for the convection and diffusion of the dissolved gas in the liquid, an inner region adjacent to the bubble and an outer region. It also involves three time scales: the shortest time scale associated with the bubble oscillation, the intermediate and the longest time scales for the mass transport across the bubble interface and the outer region, respectively. The problem is analysed systematically using matched asymptotic expansions in space and multi-scales in time, the large parameter being the Péclet number for mass transfer. A tractable mathematical model is formulated for rectified diffusion. It comprises the Rayleigh–Plesset equation for the spherical oscillations of the bubble on the shortest time scale, an ordinary differential equation for the mass of gas in the bubble on the intermediate time scale and a diffusion equation on the longest time scale and the long length scale. Predictions of this simplified model have excellent agreement with experimental results for a single spherical bubble in the bulk over millions of cycles of oscillation. A parametric study is undertaken with several new phenomena/features observed for acoustic frequencies in the range 12 to 62 kHz, acoustic field $0.12$ to $0.22$ bar and the initial gas concentration in the liquid being $0.95$ to $1.0$ of the saturation concentration.

1. Introduction

A complete explanation for rectified diffusion of spherical bubbles in the bulk does not yet exist. Rectified diffusion cannot be traced with computational fluid dynamics packages because this requires solving the spherically symmetric partial differential equations involved for billions of time steps in order to simulate millions of cycles of oscillation. It is impossible to record the detailed process of the bubble growth experimentally. As a compromise, the conditions or thresholds for rectified diffusion occurrence have been provided approximately by experiments for specific applications and specific set-ups, which do not provide systematic understanding. The number of applications of rectified diffusion grows daily and the gap in understanding holds back many important areas of technology. A clear understanding of the mechanisms providing these interesting cavitation effects is required to optimize such applications and thereby to make them safer, more precise, more efficacious and more affordable.

A bubble exchanges gas molecules across its interface with the surrounding liquid. Due to surface tension, the gas pressure inside a bubble is higher than the partial gas pressure and thus a bubble dissolves. Blake (1949) first established that, if a gas bubble undergoes large-amplitude oscillations, then it may increase in mass. Essentially, more gas diffuses into the bubble than out of the bubble across the bubble interface. The presence of ultrasonic pressure oscillations may lead to bubble growth depending on whether the amplitude of the pressure oscillations is above a threshold value (Plesset & Prosperetti 1977; Young 1989; Leighton 1994; Lauterborn & Kurz 2010; Brennen 2013). Fyrillas & Szeri (1994) consider the case where the slow time scale of diffusion can be separated from the fast time scale of the bubble oscillation, which is still considered to be the state-of-the-art work in this field (Lauterborn & Kurz 2010). Both the oscillatory and smooth problems are then treated by singular perturbation methods. Their results for bubble growth rates are quite seriously under-predicted for larger bubbles. The existing literature is limited to threshold conditions in which there is no transport of gas across the bubble interface or to small-amplitude oscillations.

Fyrillas & Szeri (1994) describe the theoretical difficulty as hinging on the determination of the net flux of dissolved gas at the bubble interface; in particular, the influence of convection associated with bubble volume oscillations. Based on the Lagrangian change of variables developed by Plesset & Zwick (1952), Smith & Wang (2022) recently found an exact solution for the gas flux at the bubble interface including the influence of convection, but only for the initial rapid bubble growth over hundreds of cycles of oscillations. Smith & Wang (2022) performed matched asymptotic expansions in space. In this article, our approach is to extend the result in Smith & Wang (2022) to allow for the simulation of millions of cycles of oscillation. We will establish the model below using both matched asymptotic expansions in space and multi-scales in time.

Recent experimental research by Ashokkumar and co-workers has improved the database of observations for the growth of a single bubble by rectified diffusion in the bulk (Lee, Kentish & Ashokkumar 2005; Leong et al. 2010, 2011). Prior to this research, the bubble radius was measured indirectly through the rise velocity of the bubble with the sound field switched off (Eller 1969; Gould 1974; Crum 1980, 1984). Lee et al. (2005) and Leong et al. (2010, 2011) used a strobe technique which enabled the bubble radius within an acoustic cycle to be imaged directly. We utilize both sets of experimental results to validate our tractable mathematical model for the rectified diffusion in the bulk of air–water systems.

The analysis of the basic model for rectified diffusion (and more sophisticated models such as the inclusion of thermal effects) has been made possible by the latest developments in singular perturbation theory. In the last few years, these state-of-the-art techniques have led to new results in the strongly nonlinear analysis of the oscillations of spherical bubbles damped by viscosity (Smith & Wang 2017), the oscillations of spherical bubbles damped by acoustic radiation (Smith & Wang 2018), the axisymmetric oscillations of liquid drops damped by viscosity (Smith 2010), the finite-amplitude travelling waves in two-dimensional plane Poiseuille flow (Smith & Wissink 2015), the finite-amplitude travelling and standing waves in two-dimensional Kolmogorov flow (Smith & Wissink 2018) and steady states and quasi-steady states in vortex dynamics (Smith & Wang 2021). In rectified diffusion, the time scales in the mass transport boundary layer adjacent to the bubble may only be systematically analysed with these techniques.

The remainder of the paper is organized as follows. The basic mathematical model for the bubble growth is described in § 2. In § 3, the basic model is analysed using strongly nonlinear analysis in time and matched asymptotic expansions in space in which the large parameter is the Péclet number for mass transfer. A self-contained summary of the resulting simplified model has been included in § 3.6 for the convenience of those readers who would prefer to bypass the asymptotic analysis. In § 4, the numerical method is discussed and numerical solutions are compared against experimental observations. In § 5, a parametric analysis is then carried out with the above theory for the applied pressure frequency, the applied pressure amplitude and the initial gas concentration in the liquid. Finally, in § 6, this study is summarized and the key outcomes are identified.

2. Physical and mathematical model

2.1. Physical model

Bubbles of approximately $100$ nm diameter are always present in most practical liquids. In almost all of the applications of sonochemistry, the bubble grows from above the threshold for rectified diffusion to the microscale (see figure 1). Bubbles take a spherical shape for the majority of this long process of rectified diffusion. This is due to two main reasons: (a) the surface tension effects are strong because they are inversely proportional to the short bubble radius; and (b) the Bjerknes force is weak because it is proportional to the small bubble volume (Wang & Blake 2010, 2011). As the spherical bubble grows larger, the Bjerknes force increases and the surface tension decreases. When the balance between these two competing aspects shifts, the bubble will become unstable to shape modes (the shape mode threshold in figure 1). Further growth only exacerbates the imbalance, the shape modes grow until the bubble eventually fragments (Brennen 2002). The smaller daughter bubbles created during the fragmentation are again spherical and the long process of rectified diffusion recommences. The spherically symmetric model is appropriate for this long process of rectified diffusion before the growth of shape modes and after the fragmentation.

Figure 1. A schematic of the typical regimes involved in rectified diffusion as a function of the bubble volume (Young 1989). The spherically symmetric model is valid for bubble volumes less than the shape mode threshold. The shape mode threshold is typically on the microscale and it follows from (4.1).

Louisnard & Gomez (2003) considered an alternative to the scenario described in figure 1. In nearly saturated liquids, subject to low acoustic frequencies and high acoustic amplitudes, the transient cavitation threshold occurs at lower bubble volumes than the shape mode threshold. Furthermore, the transient cavitation threshold may merge with the rectified diffusion threshold, which means that all of the spherical bubble growth is explosive growth.

Computational studies for non-spherical bubbles have been mainly undertaken for a few cycles of oscillation (Wang 1998, 2014; Tiwari, Pentano & Freund 2015; Lechner et al. 2017; Zeng et al. 2018), limited by computational resources, accumulated numerical errors and/or robustness of numerical modelling. Wang et al. (2022) recently simulated non-spherical bubble oscillations for dozens of cycles of oscillation. These studies are far away from simulating the millions of cycles of oscillation that are required for the modelling of rectified diffusion for non-spherical bubbles.

The physical mechanisms underlying bubble growth (or rectified diffusion) in the absence of surfactants are threefold (Fyrillas & Szeri 1994). Firstly, the gas pressure inside the bubble varies during each cycle of the bubble oscillation which changes the concentration of gas in the liquid at the bubble wall via Henry's law. Dissolved gas will diffuse into the bubble during the part of the oscillation cycle when the partial pressure of gas in the bubble is lower. Similarly, gas will diffuse into the liquid during the remainder of the cycle because the partial pressure of gas is higher. History effects make this mechanism more complicated than this simple description (Peñas-López et al. 2016), but it serves as a guide. Secondly, the rate of diffusion is proportional to the surface area. Therefore, more gas will enter the bubble when the surface area is greater than its mean value than leave the bubble when the surface area is smaller than its mean value. This second mechanism is known as the area effect (Young 1989). Thirdly, the flux of gas in the liquid is proportional to the concentration gradient of the dissolved gas. As the bubble expands, convection of the liquid surrounding the bubble causes the concentration gradient of the gas in the liquid to increase. Conversely, as the bubble contracts, convection of the adjacent liquid causes the concentration gradient of the gas in the liquid to decrease. The net effect of this convection is to promote bubble growth which slowly accumulates over many oscillations. This third mechanism is known as the shell effect (Young 1989).

Compressibility of the liquid damps the oscillation of bubbles by the emission of acoustic radiation. It is significant for large-amplitude oscillations when the time derivatives of the bubble radius become very large, usually during the collapse (Wang 2016; Smith & Wang 2018). The published experimental results for rectified diffusion do not consider such cases.

The current model is not able to take into account the effect of surfactants. Henry's law is no longer applicable in the presence of surfactants. It must be replaced by a flux boundary condition in which the non-equilibrium partitioning of the gas across the interface drives the mass transport (Fyrillas & Szeri 1995).

Our research is focussed on rectified diffusion in the bulk, that is bubble growth subject to acoustic waves, under the following hypothesis. We assume that (i) bubbles are spherical, since their sizes during rectified diffusion are typically small compared with other physical scales including the distances between a bubble and its neighbouring bubbles or boundaries, and the wavelength of acoustic waves; (ii) in the bubble, the concentration of gas is uniform (Eller & Flynn 1965); (iii) the pressure of the gas in the bubble satisfies a polytropic relation; (iv) the concentration of gas in the liquid and the flow of liquid surrounding the bubble are assumed to be spherically symmetric; (v) the bubbles do not move relative to the host liquid; (vi) the liquid is incompressible and Newtonian; (vii) the viscous heating of the liquid is assumed to take place on a longer time scale than those considered herein; and (viii) phase change at the bubble interface is neglected. The localized inaccuracy of the polytropic bubble model near resonance as discussed in Appendix C of Fyrillas & Szeri (1994) should have a limited influence on our results over many millions of bubble oscillations.

2.2. Mathematical model

The Rayleigh–Plesset equation for a gas bubble in an incompressible Newtonian liquid is

(2.1)\begin{equation} \bar{R} \frac{{{\rm d}}^{2} \bar{R}}{{{\rm d}} \bar{t}^{2}} + \frac{3}{2} \left( \frac{{{\rm d}} {\bar{R}}}{{{\rm d}} {\bar{t}}}\right)^{2} =\frac{\bar{p}_l}{\rho}, \end{equation}

in which $\bar {p}_l$ is the pressure of the liquid at the bubble surface minus the external pressure in the liquid (Prosperetti 1984)

(2.2)\begin{equation} \bar{p}_l = \bar{p}_{g} - \frac{2 \sigma}{\bar{R}} -(\bar{p}_{\infty} - \bar{p}_v) - \frac{4 \mu}{\bar{R}} \frac{{{\rm d}} {\bar{R}}}{{{\rm d}} {\bar{t}}}-\bar{p}_a \sin (2 {\rm \pi}\bar{\nu}_a \bar{t}), \end{equation}

and $\bar {p}_{g}$ is the partial pressure of the gas inside the bubble

(2.3)\begin{equation} \bar{p}_{g} = \bar{p}_{g0} \left ( \frac{\bar{m}}{\bar{m}_i} \right )^{\kappa} \left( \frac{\bar{R}_{eqi}}{\bar{R}} \right)^{3 \kappa}, \end{equation}

where $\bar {R}(\bar {t})$ is the spherical bubble radius at time $\bar {t}$, $\bar {m} (\bar {t})$ is the mass of gas in the bubble, $\bar {R}_{eqi}$ the initial equilibrium bubble radius, $\bar {m}_i$ the initial mass of gas in the bubble, $\rho$ the liquid density, $\bar {p}_\infty$ the hydrostatic pressure of the liquid, $\bar {p}_v$ the vapour pressure of the liquid, $\bar {p}_{g0}$ the initial partial pressure of the bubble gas, $\kappa > 1$ the polytropic index, $\sigma$ the surface tension, $\mu$ the liquid viscosity and $\bar {p}_a$ ($\bar {\nu }_a$) is the amplitude (frequency) of the external pressure. The concentration of gas in the liquid, $\bar {c}(\bar {r}, \bar {t})$, is governed by a spherically symmetric convection–diffusion equation in a domain with a moving boundary of the form

(2.4)\begin{equation} \frac{\partial {\bar{c}}}{\partial {\bar{t}}} + \frac{{{\rm d}} {\bar{R}}}{{{\rm d}} {\bar{t}}} \left ( \frac{\bar{R}}{\bar{r}} \right )^{2} \frac{\partial {\bar{c}}}{\partial {\bar{r}}} = \frac{D}{\bar{r}^{2}} \frac{\partial}{\partial {\bar{r}}} \left( \bar{r}^{2} \frac{\partial {\bar{c}}}{\partial {\bar{r}}} \right)\quad \text{for} \ \bar{r} > \bar{R}(\bar{t}), \end{equation}

where $\bar {r}$ is the radial distance from the centre of the bubble and $D$ is the mass diffusivity (Eller & Flynn 1965). The far-field boundary condition is

(2.5)\begin{equation} \bar{c} \rightarrow \bar{C}_i \quad \text{as} \ \bar{r} \rightarrow \infty, \end{equation}

where $\bar {C}_i$ is the initial uniform concentration of the gas in the liquid. The boundary condition at the bubble interface is provided by Henry's law, which relates the concentration of the gas in the liquid at the bubble interface, $\bar {c} (\bar {R}(\bar {t}), \bar {t})$, to the partial pressure of the gas inside the bubble, $\bar {p}_{g}$ (Eller & Flynn 1965; Fyrillas & Szeri 1994; Brennen 2013). We have

(2.6)\begin{equation} \bar{c} (\bar{R}(\bar{t}), \bar{t}) = \frac{\bar{p}_{g}}{k_H} = \frac{1}{k_H} \bar{p}_{g0} \left ( \frac{\bar{m}}{\bar{m}_i} \right )^{\kappa} \left( \frac{\bar{R}_{eqi}}{\bar{R}} \right)^{3 \kappa}, \end{equation}

where $k_H$ is Henry's law constant. The mass of gas in the bubble must satisfy the conservation law

(2.7)\begin{equation} \frac{{{\rm d}} {\bar{m}}}{{{\rm d}} {\bar{t}}} = 4 {\rm \pi}\bar{R}^{2} D \frac{\partial {\bar{c}}}{\partial {\bar{r}}}(\bar{R}(\bar{t}),\bar{t}). \end{equation}

The initial conditions are given by

(2.8a–d)\begin{equation} \bar{R}(0) = \bar{R}_{eqi}, \quad\frac{{{\rm d}} {\bar{R}}}{{{\rm d}} {\bar{t}}}(0) = 0, \quad \bar{c} ( \bar{r}, 0 ) = \bar{C}_i, \quad \bar{m}(0) = \bar{m}_i, \end{equation}

for $\bar {r} > \bar {R}_{eqi}$. As we are only interested in solutions which are periodic on the time scale for bubble oscillation, we shall also impose periodicity on this time scale. The physical system (2.1)–(2.8a–d) has a global statement of conservation of gas in the form

(2.9)\begin{equation} 4 {\rm \pi}\int_{\bar{r}=\bar{R}(\bar{t})}^{\infty} ( \bar{c} - \bar{C}_i ) \bar{r}^{2} \; {{\rm d}} \bar{r} = \bar{m}_i - \bar{m}. \end{equation}

In summary, the basic model comprises the following four aspects (Epstein & Plesset 1950; Eller & Flynn 1965; Fyrillas & Szeri 1994). (i) The time dependence of the spherical bubble radius is modelled by the Rayleigh–Plesset equation, which is a nonlinear ordinary differential equation. Both viscous effects and surface tension will be considered. (ii) Henry's law will be used to model the concentration of gas at the bubble surface. (iii) The concentration of gas in the liquid is governed by a spherically symmetric convection–diffusion equation, which is a nonlinear partial differential equation in the moving liquid domain. (iv) The mass of gas in the bubble satisfies a conservation law given by the gas flux across the bubble surface, which is a second nonlinear ordinary differential equation. The mass of the gas couples into the other equations via the partial pressure of the gas, this being the only significant difference between this mathematical model and the model solved by Fyrillas & Szeri (1994). The initial concentration of gas in the liquid is uniform, which is also the boundary condition for the concentration of gas at infinity. This mathematical model simulates bubble growth, bubble dissolution and bubble oscillation.

Equations (2.1)–(2.9) are scaled using $\bar {R} = \bar {R}_{eqi} \hat {R}$, $\bar {c} = \bar {C}_i + \bar {C}_i \hat {c}$, $\bar {m} = \bar {m}_i \hat {m}$,

(2.10a,b)\begin{equation} \bar{t} = \frac{\bar{R}_{eqi}}{U} \hat{t},\quad \bar{r} = \bar{R}_{eqi} r, \end{equation}

where

(2.11a,b)\begin{equation} \Delta = \bar{p}_{\infty} - \bar{p}_v, \quad U = \sqrt{\frac{\Delta}{\rho}}, \end{equation}

in which $\Delta$ is the characteristic pressure of the liquid and $U$ is a reference velocity. The overbars denote dimensional dependent variables and the hats their dimensionless counterparts. The dimensionless Rayleigh–Plesset equation takes the form

(2.12)\begin{equation} \hat{R} \frac{{{\rm d}}^{2} \hat{R}}{{{\rm d}} \hat{t}^{2}} + \frac{3}{2} \left( \frac{{{\rm d}} {\hat{R}}}{{{\rm d}} {\hat{t}}} \right)^{2} = \frac{p_{g0} \hat{m}^{\kappa}}{\hat{R}^{3 \kappa}} - \frac{2}{We \hat{R}} - 1 - \frac{4}{Re \hat{R}} \frac{{{\rm d}} {\hat{R}}}{{{\rm d}} {\hat{t}}} - p_a \sin (2 {\rm \pi}\nu_a \hat{t}), \end{equation}

in which $Re = \rho U \bar {R}_{eqi}/ \mu$ is the Reynolds number, $We = \bar {R}_{eqi} \Delta / \sigma$ is the Weber number, $p_{g0} = \bar {p}_{g0} / \Delta$ is the dimensionless initial partial pressure of the bubble gases and $p_a = \bar {p}_a / \Delta$ ($\nu _a = \bar {\nu }_a \bar {R}_{eqi} / U$) is the dimensionless amplitude (frequency) of the applied external pressure. The convection–diffusion equation for the gas concentration in the liquid becomes

(2.13)\begin{equation} \frac{\partial {\hat{c}}}{\partial {\hat{t}}} + \frac{{{\rm d}} {\hat{R}}}{{{\rm d}} {\hat{t}}} \left ( \frac{\hat{R}}{r} \right )^{2} \frac{\partial {\hat{c}}}{\partial {r}} = \frac{1}{Pe}\frac{1}{r^{2}} \frac{\partial}{\partial {r}} \left ( r^{2} \frac{\partial {\hat{c}}}{\partial {r}} \right ) \quad \text{for} \ r > \hat{R}(\hat{t}), \end{equation}

where $Pe = U \bar {R}_{eqi} / D$ is the large Péclet number for mass transfer. The boundary conditions are given by

(2.14)\begin{equation} \hat{c} \rightarrow 0 \quad \text{as} \ r \rightarrow \infty, \end{equation}

and

(2.15)\begin{equation} \hat{c} (\hat{R}(\hat{t}), \hat{t}) ={-} 1 + d p_{g0} \frac{\hat{m}^{\kappa}}{\hat{R}^{3 \kappa}}, \end{equation}

in which $d = \Delta / k_H \bar {C}_i$. Conservation of the mass of gas in the bubble becomes

(2.16)\begin{equation} \frac{{{\rm d}} {\hat{m}}}{{{\rm d}} {\hat{t}}} = \frac{4 {\rm \pi}M}{Pe} \hat{R}^{2} \frac{\partial {\hat{c}}}{\partial {r}} (\hat{R}(\hat{t}),\hat{t}), \end{equation}

where

(2.17)\begin{equation} M = \frac{\bar{R}_{eqi}^{3} \bar{C}_i}{\bar{m}_i}. \end{equation}

The initial conditions are given by

(2.18a–d)\begin{equation} \hat{R}(0) = 1, \quad \frac{{{\rm d}} {\hat{R}}}{{{\rm d}} {\hat{t}}}(0) = 0, \quad \hat{c} ( r, 0 ) = 0, \quad \hat{m}(0) = 1, \end{equation}

for $r > 1$. Finally, we impose periodicity on the time scale for bubble oscillation. The dimensionless global conservation of gas becomes

(2.19)\begin{equation} 4 {\rm \pi}M \int_{r=\hat{R}(\hat{t})}^{\infty} \hat{c} r^{2} \, {{\rm d}}r = 1 - \hat{m}. \end{equation}

Equations (2.12)–(2.18a–d) constitute the dimensionless version of the mathematical model for our three unknowns $\hat {R}(\hat {t})$, $\hat {c}(r,\hat {t})$ and $\hat {m}(\hat {t})$ to be studied herein.

In principle, (2.12)–(2.18a–d) can be solved numerically. Firstly, the initial boundary value problem (2.13)–(2.15) and (2.18a–d) may be solved for the concentration of gas in the liquid. Secondly, the initial value problem (2.16) and (2.18a–d) may be solved for the mass of gas in the bubble. Lastly, the initial value problem (2.12) and (2.18a–d) may be solved for the bubble radius. However, the process takes millions of cycles of oscillation and the numerical simulations require billions of time steps. While it is practical to accurately integrate the two ordinary differential equations for bubble radius and bubble mass for billions of time steps, unfortunately, it is impractical to accurately integrate the initial boundary value problem for the concentration of gas in the liquid (2.13)–(2.15) with the appropriate initial conditions in (2.18a–d) for billions of time steps even with modern supercomputers. To resolve this challenge, we will develop a tractable model for rectified diffusion in § 3.

3. Asymptotic analysis

The problem is characterized by having several physical processes: the inertial bubble oscillation, the mass transport at the bubble surface and the convection and diffusion of the dissolved gas in the liquid. These processes act simultaneously and interactively with different scales.

The typical value of the Péclet number for mass transfer, $Pe$, is greater than one hundred thousand. Therefore, this problem is susceptible to perturbation theory. We divide the liquid domain into an inner region and an outer region, as illustrated in figure 2. The inner region is the concentration boundary layer adjacent to the bubble surface with a thickness ${{O}}(Pe^{-1/2})$, where $r=\hat {R}(\hat {t})+{{O}}(Pe^{-1/2})$. The outer region is outside the concentration boundary layer, where $r=\hat {R}(\hat {t})+{{O}}(1)$. We now study these two regions and the asymptotic matching between them (Bender & Orszag 1987).

Figure 2. A schematic of the spatial asymptotic structure. The bubble interface, $r=\hat {R}(\hat {t})$, oscillates on the shortest time scale and grows/dissolves due to the mass flux on the intermediate time scale. The concentration of gas in the inner region oscillates on the shortest time scale as governed by Henry's law and convects on the shortest time scale due to the bubble motion. The concentration of gas in the outer region convects on the shortest time scale due to the bubble motion, but diffusion of gas only takes place on the longest time scale.

The dimensionless diffusion time scale over a dimensionless length scale $\hat {L}$ is $Pe\, \hat {L}^{2}$. The dimensionless diffusion time scales are ${{O}}(1)$ in the inner region with a thickness ${{O}}(Pe^{-1/2})$ and ${{O}}(Pe)$ in the outer region with a length scale of ${{O}}(1)$. In addition, it is known from (2.16) that the time scale for the mass transfer across the bubble surface is ${{O}}(Pe^{1/2})$.

The analysis is organized as follows. In § 3.1, a Lagrangian change of independent variables is introduced and multiple scale analyses are applied in the inner region. Multiple scale analyses are also applied in the outer region in § 3.2. Matching between the inner and outer regions is undertaken in § 3.3. The gradient of the concentration of gas at the bubble interface is determined in § 3.4. In § 3.5, the far-field concentration of gas in the outer liquid region is evaluated in the form of a similarity solution.

3.1. Inner region

3.1.1. Leading-order problem

In the inner region, $r=\hat {R}(\hat {t})+{{O}}(Pe^{-1/2})$, we introduce the following Lagrangian change of independent variables ($x={{O}}(1$)) first used by Plesset & Zwick (1952):

(3.1)\begin{equation} x = \frac{Pe^{1/2}}{3} [r^{3} - \hat{R}(\hat{t})^{3}], \quad \tau = \hat{t}, \end{equation}

where the new dependent variables are the bubble radius $R(\tau ) = \hat {R}(\hat {t})$, the concentration of gas in the inner liquid region $c(x,\tau ) = \hat {c}(r,\hat {t})$ and the mass of gas in the bubble $m(\tau ) = \hat {m}(\hat {t})$. The dependent variables without bars or hats correspond to dimensionless variables in the Lagrangian frame of reference. The convection and diffusion equation (2.13) for the concentration of gas in the liquid may be rewritten as a diffusion equation

(3.2)\begin{equation} \frac{\partial {c}}{\partial {\tau}} = \frac{\partial}{\partial {x}} \left [ \left ( R^{3} + \frac{3 x}{Pe^{1/2}} \right)^{4/3} \frac{\partial {c}}{\partial {x}} \right ] \quad \text{for} \ x > 0. \end{equation}

The Rayleigh–Plesset equation (2.12) for the bubble radius becomes

(3.3)\begin{equation} R \frac{{{\rm d}}^{2} R}{{{\rm d}} \tau^{2}} + \frac{3}{2} \left( \frac{{{\rm d}} {R}}{{{\rm d}} {\tau}} \right)^{2} = \frac{p_{g0} m^{\kappa}}{R^{3 \kappa}} - \frac{2}{We R} - 1 - \frac{4}{Re R} \frac{{{\rm d}} {R}}{{{\rm d}} {\tau}} - p_a \sin (2 {\rm \pi}\nu_a \tau). \end{equation}

Equation (2.16) for conservation of mass of the gas in the bubble transforms to

(3.4)\begin{equation} \frac{{{\rm d}} {m}}{{{\rm d}} {\tau}} = \frac{4 {\rm \pi}M}{Pe^{1/2}} R^{4} \frac{\partial {c}}{\partial {x}} (0,\tau). \end{equation}

The growth of bubbles in a fluid is modelled using three time scales. The concentration of gas in the liquid, $c$, varies over a short time scale $t$ corresponding to the period of inertial bubble oscillation, an intermediate time scale $\varLambda$ associated with the time scale for mass transfer across the bubble surface and a long time scale $\lambda$ which follows from diffusion on the long length scale of the outer region. Accordingly, we introduce three time variables. The short time variable for the bubble oscillation is

(3.5a)\begin{equation} t = \tau; \end{equation}

the intermediate time variable of the mass transfer at the bubble surface is

(3.5b)\begin{equation} \varLambda = \frac{\tau}{Pe^{1/2}}; \end{equation}

and the long time variable for the diffusion of the gas in the outer liquid region is

(3.5c)\begin{equation} \lambda = \frac{\tau}{Pe}. \end{equation}

As real time $\tau$ increases, the short time scale $t$ increases at the same rate, while the intermediate time scale $\varLambda$ increases slowly and the long time scale $\lambda$ increases slowest. Thus, we have

(3.6)\begin{equation} \frac{\partial}{\partial {\tau}} = \frac{\partial}{\partial {t}}+ \frac{1}{Pe^{1/2}} \frac{\partial}{\partial {\varLambda}}+ \frac{1}{Pe} \frac{\partial}{\partial {\lambda}}. \end{equation}

Our choice of a short time scale will leave the period of the leading-order solution varying on the longer time scales. This is inconsistent with formal perturbation theory; however, there will be no changes to our subsequent results and it is convenient to retain our short time scale.

We introduce expansions of the first two orders, since the disturbances are ${{O}}(Pe^{-1/2})$ as shown in (3.2), (3.3) and (3.4):

(3.7a–c)\begin{equation} c \sim c_0 + \frac{c_1}{Pe^{1/2}}, \quad R \sim R_0 + \frac{R_1}{Pe^{1/2}}, \quad m \sim m_0 + \frac{m_1}{Pe^{1/2}}, \end{equation}

as $Pe \rightarrow \infty$. In the above, the subscript $0$ corresponds to the leading-order approximation and subscript $1$ to the first correction. Substituting into the governing equations (3.2)–(3.4) and comparing coefficients of $Pe$, we find a sequence of problems. At leading order, (3.2)–(3.4) become

(3.8)$$\begin{gather} \frac{\partial {c_0}}{\partial {t}} = R_0^{4} \frac{\partial^{2} {c_0}}{\partial {x}^{2}}, \end{gather}$$

(3.9)$$\begin{gather}R_0 \frac{\partial^{2} {R_0}}{\partial {t}^{2}} + \frac{3}{2} \left( \frac{\partial {R_0}}{\partial {t}} \right)^{2} = \frac{p_{g0} m_0^{\kappa}}{R_0^{3 \kappa}} - \frac{2}{We R_0} - 1 - \frac{4}{Re R_0} \frac{\partial {R_0}}{\partial {t}} - p_a \sin (2 {\rm \pi}\nu_a t), \end{gather}$$

(3.10)$$\begin{gather}\frac{\partial {m_0}}{\partial {t}} = 0. \end{gather}$$

Equation (3.10) shows that the mass of gas in the bubble does not change with the fast time but changes in $\varLambda$ and $\lambda$ through the accumulated effect of weak diffusion. It is readily integrated to yield $m_0=m_0(\varLambda, \lambda )$. Equation (3.8) is not straightforward to integrate and will be solved below in § 3.4. The boundary condition at the bubble interface (2.15) when $x=0$ is given by

(3.11)\begin{equation} c_0 (0, t, \varLambda,\lambda) = f(t, \varLambda,\lambda) \equiv{-} 1 + d p_{g0} \frac{m_0^{\kappa}}{R_0^{3 \kappa}}. \end{equation}

A second boundary condition will be determined by matching with the outer region in § 3.3. The initial conditions from (2.18a–d) are given by

(3.12a–d)\begin{equation} R_0(0,0,0) = 1, \quad \frac{\partial {R_0}}{\partial {t}}(0,0,0) = 0, \quad c_0 ( x, 0,0,0 ) = 0, \quad m_0(0,0) = 1 \end{equation}

for $x > 0$.

As $m_0=m_0(\varLambda, \lambda )$ does not vary over the short time scale $t$, the bubble oscillates repeatedly according to (3.9) over this time scale. As $R_0$ is a periodic function on the short time scale $t$, the concentration of the dissolved gas $c_0$ at a fixed point in the liquid domain oscillates repeatedly according to (3.8). The period of this oscillation $P$ is a function of the intermediate time scale $\varLambda$ and long time scale $\lambda$ for an initial transient, after which it settles down to the acoustic frequency.

3.1.2. First correction

In the leading-order analysis, we obtain the variation on the short time scale, but not the variation on the intermediate time scale. For this purpose, we will discuss the first-order correction. The equations at next order ${{O}}(Pe^{-1/2})$ for (3.2)–(3.4) are complicated. They are expressed in the concise matrix format. At next order, we obtain

(3.13) \begin{equation} L \boldsymbol{z} = \begin{pmatrix} \displaystyle -\dfrac{\partial {c_0}}{\partial {\varLambda}} + 4 R_0 \dfrac{\partial}{\partial {x}} \left \{ x \dfrac{\partial {c_0}}{\partial {x}} \right \}\\ \displaystyle -2 R_0 \dfrac{\partial^{2} R_0}{\partial t \partial \varLambda} - \dfrac{\partial {R_0}}{\partial {\varLambda}} \left \{ 3 \dfrac{\partial {R_0}}{\partial {t}} + \dfrac{4}{Re R_0} \right \}\\ \displaystyle -\dfrac{\partial {m_0}}{\partial {\varLambda}} + 4 {\rm \pi}M R_0^{4} \dfrac{\partial {c_0}}{\partial {x}} (0,t,\varLambda,\lambda) \end{pmatrix}, \end{equation}

in which the matrix differential operator $L$ and unknown vector function $z$ are

(3.14) \begin{equation} L=\left(\begin{array}{@{}ccc@{}} \displaystyle \dfrac{\partial}{\partial {t}} - R_0^{4} \dfrac{\partial^{2}}{\partial x^{2}} & \displaystyle -4 R_0^{3} \dfrac{\partial^{2} {c_0}}{\partial {x}^{2}} & \displaystyle 0\\ \displaystyle 0 & \displaystyle \bar{L} & \displaystyle -\dfrac{p_{g0} \kappa m_0^{\kappa-1}}{R_0^{3 \kappa}}\\ \displaystyle 0 & \displaystyle 0 & \displaystyle \dfrac{\partial}{\partial {t}} \end{array}\right),\quad \boldsymbol{z}=\left(\begin{array}{@{}c@{}} c_1\\ R_1\\ m_1 \end{array}\right). \end{equation}

In the operator $L$, $\bar {L}$ is defined as

(3.15)\begin{align} \bar{L} R_1 &= R_0 \frac{\partial^{2} {R_1}}{\partial {t}^{2}} + \frac{\partial {R_1}}{\partial {t}} \left [ 3 \frac{\partial {R_0}}{\partial {t}} + \frac{4}{Re R_0} \right ]\nonumber\\ &\quad + R_1 \left [ \frac{\partial^{2} {R_0}}{\partial {t}^{2}} + \frac{3 p_{g0} \kappa m_0^{\kappa}}{R_0^{3 \kappa + 1}} - \frac{2}{We R_0^{2}} - \frac{4}{Re R_0^{2}} \frac{\partial {R_0}}{\partial {t}} \right ]. \end{align}

As discussed in § 2, the periodicity conditions on the shortest time scale corresponds to

(3.16)$$\begin{gather} c_1 (x, t, \varLambda, \lambda) = c_1 (x, t+P, \varLambda, \lambda), \end{gather}$$

(3.17)$$\begin{gather}{[}R_1,m_1](t,\varLambda, \lambda) = [R_1,m_1](t+P,\varLambda, \lambda). \end{gather}$$

These periodicity conditions suppress the secular terms in (3.13) and are necessary to ensure the uniformity of the asymptotic expansions. There are also boundary conditions on $c_1$, but we do not specify these conditions as we will not require them in what follows.

We determine an adjoint matrix differential equation operator $L^{*}$ with unknown vector function $\boldsymbol {r}$:

(3.18) \begin{equation} L^{*}=\left(\begin{array}{@{}ccc@{}} \displaystyle -\dfrac{\partial}{\partial {t}} - R_0^{4} \dfrac{\partial^{2}}{\partial x^{2}} & \displaystyle 0 & \displaystyle 0 \\ \displaystyle -4 R_0^{3} \dfrac{\partial^{2} {c_0}}{\partial {x}^{2}} & \displaystyle \hat{L} & \displaystyle 0\\ \displaystyle 0 & \displaystyle -\dfrac{p_{g0} \kappa m_0^{\kappa-1}}{R_0^{3 \kappa}} & \displaystyle -\dfrac{\partial}{\partial {t}} \end{array}\right),\quad \boldsymbol{r}=\left(\begin{array}{@{}c@{}} \displaystyle \alpha\\ \displaystyle \beta\\ \displaystyle \gamma \end{array}\right), \end{equation}

such that the right-hand side of the following equation is in the form

(3.19) \begin{align} &\boldsymbol{r}^{T} L \boldsymbol{z} - \boldsymbol{z}^{T} L^{*} \boldsymbol{r}= \frac{\partial}{\partial {t}} \left[ \alpha c_1 + \beta R_0 \frac{\partial {R_1}}{\partial {t}} - R_1 \frac{\partial}{\partial {t}}(\beta R_0) \right . \nonumber\\ &\quad + \left. 3 \beta \frac{\partial {R_0}}{\partial {t}} R_1 +\frac{4 \beta R_1}{Re R_0} + \gamma m_1 \right] +\frac{\partial}{\partial {x}} \left[ - \alpha R_0^{4} \frac{\partial {c_1}}{\partial {x}} + c_1 R_0^{4} \frac{\partial {\alpha}}{\partial {x}} \right]. \end{align}

The left-hand side of (3.19) takes a standard form adopted when seeking conditions for which the first correction (3.13) has solutions. In the operator $L^{*}$, $\hat {L}$ is defined as

(3.20)\begin{align} \hat{L} \beta &= \frac{\partial^{2}}{\partial t^{2}} (\beta R_0) - 3 \frac{\partial}{\partial {t}} \left ( \beta \frac{\partial {R_0}}{\partial {t}} \right ) - \frac{4}{Re} \frac{\partial}{\partial {t}} \left ( \frac{\beta}{R_0} \right )\nonumber\\ &\quad + \beta \left [ \frac{\partial^{2} {R_0}}{\partial {t}^{2}} + \frac{3 p_{g0} \kappa m_0^{\kappa}}{R_0^{3 \kappa + 1}} - \frac{2}{We R_0^{2}} - \frac{4}{Re R_0^{2}} \frac{\partial {R_0}}{\partial {t}} \right ]. \end{align}

At this point of the analysis, it would usually be necessary to specify the periodicity and boundary conditions on $\boldsymbol {r}$. We do not require the periodicity and boundary conditions in this case, as the problem simplifies. In order to suppress secular terms on the right-hand side of the first correction (3.13), we require solutions of $L^{*} \boldsymbol {r} = \boldsymbol {0}$. The second equation in $L^{*} \boldsymbol {r} = \boldsymbol {0}$ is only satisfied when $\alpha =0$, because the coefficient of $\alpha$ is spatially dependent, whereas the remaining terms are independent of space. The only solution of $\hat {L} \beta =0$ that we have been able to find is $\beta =0$. Therefore,

(3.21)\begin{equation} \boldsymbol{r} = (0,0,1)^{\rm T} \end{equation}

is the one linearly independent solution of $L^{*} \boldsymbol {r} = \boldsymbol {0}$. Using this solution, (3.19) simplifies to

(3.22)\begin{equation} -\frac{\partial {m_0}}{\partial {\varLambda}} + 4 {\rm \pi}M R_0^{4} \frac{\partial {c_0}}{\partial {x}} (0,t, \varLambda, \lambda) = \frac{\partial {m_1}}{\partial {t}}. \end{equation}

In order to suppress the secular terms in (3.13), we use the periodicity conditions (3.17) to establish that

(3.23)\begin{equation} \left\langle \frac{\partial {m_1}}{\partial {t}} \right\rangle = 0, \end{equation}

in which

(3.24)\begin{equation} \langle \,\cdot\, \rangle = \frac{1}{P} \int^{t+P}_{t} \cdot\, {{\rm d}}t \end{equation}

denotes the average value over a single cycle of oscillation, where $P$ is the period. Therefore, we obtain

(3.25)\begin{equation} \frac{\partial {m_0}}{\partial {\varLambda}} = 4 {\rm \pi}M \left \langle R_0^{4} \frac{\partial {c_0}}{\partial {x}} (0,t,\varLambda,\lambda) \right \rangle. \end{equation}

The mass of gas in the bubble varies following (3.25). The bubble is being compelled to retain only a short-term memory of the concentration gradient by integrating on the order-one time scale. This equation will be used to update the bubble gas mass on the intermediate time scale $\varLambda$.

3.2. Outer region

In the outer expansion, $r=\hat {R}(\hat {t})+{{O}}(1)$, we introduce the following Lagrangian change of variables ($y={{O}}(1$)):

(3.26)\begin{equation} y = \frac{x}{Pe^{1/2}} = \frac{1}{3} [r^{3} - \hat{R}(\hat{t})^{3}],\quad \tau = \hat{t}. \end{equation}

The concentration of gas in the outer liquid region, $C$, is described using two time scales on this long length scale. We therefore introduce two time variables, $t=\tau$ corresponding to the period of inertial bubble oscillation and $\lambda =\tau /Pe$ associated with the diffusion time scale in the outer region, as follows:

(3.27)\begin{equation} \frac{\partial}{\partial {\tau}} = \frac{\partial}{\partial {t}}+ \frac{1}{Pe} \frac{\partial}{\partial {\lambda}}. \end{equation}

The convection–diffusion equation for mass transport (2.13) for $C(y, t, \lambda ) = \hat {c}(r, \hat {t})$ becomes

(3.28)\begin{equation} \frac{\partial {C}}{\partial {t}}+ \frac{1}{Pe} \frac{\partial {C}}{\partial {\lambda}} = \frac{1}{Pe} \frac{\partial}{\partial {y}} \left [( R^{3} + 3 y )^{4/3} \frac{\partial {C}}{\partial {y}} \right ] \quad \text{for} \ y > 0 \end{equation}

with far-field boundary condition

(3.29)\begin{equation} C \rightarrow 0 \quad \text{as} \ y \rightarrow \infty, \end{equation}

the initial condition

(3.30)\begin{equation} C = 0 \quad \text{on} \ t=\lambda=0 \end{equation}

and the periodicity condition

(3.31)\begin{equation} C (y, t, \lambda) = C (y, t+P, \lambda) \end{equation}

as discussed in § 2. A second boundary condition must be determined by matching with the inner region. Equation (3.28) shows that the diffusion in the outer region is of the order of ${{O}}(Pe^{-1})$ and it justifies the choice of the long time scale $\lambda$.

Since the disturbance is ${{O}}(Pe^{-1})$ as shown in (3.28), we introduce an expansion of the following form:

(3.32a,b)\begin{equation} C \sim C_0 + \frac{C_1}{Pe}, \quad R \sim R_0, \end{equation}

as $Pe \rightarrow \infty$. We will need to study (3.28) at both the leading order and the first correction, so two terms are required in the expansion for the concentration of gas in the outer liquid region $C$. Only the leading-order term in the bubble radius $R$ is required in the analysis which follows. At leading order in (3.28), we obtain

(3.33)\begin{equation} \frac{\partial {C_0}}{\partial {t}}=0, \end{equation}

which is readily integrated to yield $C_0 = C_0 (y, \lambda )$.

To determine $C_0 (y, \lambda )$, we need to consider the first-order correction. At next order in (3.28), we have

(3.34)\begin{equation} \frac{\partial {C_1}}{\partial {t}}+ \frac{\partial {C_0}}{\partial {\lambda}} = \frac{\partial}{\partial {y}} \left [ ( R_0^{3} + 3 y )^{4/3} \frac{\partial {C_0}}{\partial {y}} \right ] \quad \text{for} \ y > 0. \end{equation}

From (3.31), $C_1$ is periodic in the short time scale $t$. We thus have the periodicity condition

(3.35)\begin{equation} C_1 (y, t, \lambda) = C_1 (y, t+P, \lambda). \end{equation}

In order to suppress the secular terms in (3.34) on the shortest time scale, we use the periodicity condition (3.35) to establish that

(3.36)\begin{equation} \left \langle \frac{\partial {C_1}}{\partial {t}} \right \rangle = 0. \end{equation}

Therefore, we obtain the following secularity condition:

(3.37)\begin{equation} \frac{\partial {C_0}}{\partial {\lambda}} = \frac{\partial}{\partial {y}} \left [ \langle ( R_0^{3} + 3 y )^{4/3} \rangle \frac{\partial {C_0}}{\partial {y}} \right ] \quad \text{for} \ y > 0 ,\end{equation}

with the far-field boundary condition

(3.38)\begin{equation} C_0 \rightarrow 0 \quad \text{as} \ y \rightarrow \infty, \end{equation}

which follows from (3.29) and the initial condition

(3.39)\begin{equation} C_0 (y,0) = 0, \end{equation}

which follows from (3.30).

In the outer region, the concentration of gas in the liquid has been shown to vary on the long time scale $\lambda$. Fyrillas & Szeri (1994) obtained the same equation for their smooth problem on the same length and time scales. Fyrillas & Szeri (1994) then solved this equation when the time derivatives are zero, which corresponds to near the threshold of rectified diffusion. They subsequently found that their theoretical results are accurate in comparison with the experimental threshold of Eller (1969). The theoretical predictions of threshold (Fyrillas & Szeri 1994) have also been successful in applications related to stable single bubble sonoluminescence (Brenner, Hilgenfeldt & Lohse 2002).

The concentration $C_0 (y, \lambda )$ is governed by the diffusion equation (3.37) with a boundary condition at infinity (3.38) and initial condition at $\lambda =0$ (3.39). Another boundary condition is still needed at $y = 0$, which will be obtained by matching with the inner solution, to be discussed in § 3.3.

3.3. Matching

Firstly, we consider the outer problem for $y \ll 1$. In this limit, (3.37) takes the form

(3.40)\begin{equation} \frac{\partial}{\partial {y}} \left [ \langle ( R_0^{3} + 3 y )^{4/3} \rangle \frac{\partial {C_0}}{\partial {y}} \right ] = 0. \end{equation}

We solve this equation to find that the matching condition is

(3.41)\begin{equation} C_0 (y, \lambda) \sim g(\lambda) \hspace{3mm} \text{as} \ y \rightarrow 0, \end{equation}

where $g(\lambda )$ remains to be determined and $g(0)=0$ due to the initial condition (3.39). Therefore, according to Prandtl's limit matching technique, the matching condition for the inner problem is

(3.42)\begin{equation} c_0 (x, t, \varLambda,\lambda) \rightarrow g(\lambda) \quad \text{as} \ x \rightarrow \infty. \end{equation}

As such, the concentration in the transition region only depends on the long time scale.

In the absence of a general solution to the inner problem, we adopt the global statement for conservation of mass (2.19) in the form

(3.43)\begin{equation} 4 {\rm \pi}M \left \{ \frac{1}{Pe^{1/2}} \int_{x=0}^{\infty} c \, {{\rm d}}\kern0.7pt x + \int_{y=0}^{\infty} C \, {{\rm d}}y \right \} = 1-m, \end{equation}

where we assume that the integrals in the inner and outer regions converge to finite values. The liquid domain here is split in terms of $r$ from $\hat {R}$ to greater than $\hat {R}+{{O}}(Pe^{-1/2})$ for the inner region and from $\hat {R}+o(1)$ to the far field for the outer region. This corresponds to $x$ from $0$ to $\infty$ and $y$ from $0$ to $\infty$, respectively. At leading order, global conservation of mass becomes

(3.44)\begin{equation} 4 {\rm \pi}M \int_{y=0}^{\infty} C_0 \, {{\rm d}}y = 1-m_0. \end{equation}

Equation (3.44) allows us to evaluate the matching condition $g(\lambda )$ in the transition region between inner and outer regions.

3.4. Leading-order inner solution

Our approach is to find an exact solution to the diffusion equation (3.8) for the leading-order concentration of gas in the liquid. Using this exact solution, the flux of gas at the bubble surface may be evaluated. We consider the variable ${\mathcal {C}} (x, t, \varLambda,\lambda ) = c_0 (x, t, \varLambda,\lambda ) - g(\lambda )$ which satisfies the initial boundary value problem

(3.45)\begin{equation} \frac{\partial {{\mathcal{C}}}}{\partial {t}} = R_0^{4} \frac{\partial^{2} {{\mathcal{C}}}}{\partial {x}^{2}} \quad \text{for} \ x > 0, \end{equation}

subject to the boundary conditions

(3.46)\begin{equation} {\mathcal{C}} \rightarrow 0 \quad \text{as} \ x \rightarrow \infty, \quad {\mathcal{C}}(0,t, \varLambda,\lambda) = f(t, \varLambda,\lambda)-g(\lambda) ,\end{equation}

and the initial condition

(3.47)\begin{equation} {\mathcal{C}}(x,0, \varLambda,\lambda)=0 \quad \text{for} \ x>0. \end{equation}

The far-field boundary condition in (3.46) follows from the matching condition (3.42) and the boundary condition at the bubble interface follows from (3.11). The initial condition (3.47) has been chosen to be consistent with the matching condition (3.42).

The following is a modified version of the initial-stage analysis described by Smith & Wang (2022). The initial boundary value problem (3.45)–(3.47) is nonlinear, but the bubble radius in (3.45) only depends on time. As the spatial domain is semi-infinite, the spectrum of this problem is continuous. We take a Fourier sine transform of (3.45)–(3.47) in space to obtain the initial value problem

(3.48)\begin{equation} \frac{\partial {\tilde{c}}}{\partial {t}} + s^{2} R_0^{4} \tilde{c} = s R_0^{4} (f-g) \quad \text{with} \ \tilde{c}(s,0, \varLambda, \lambda)=0, \end{equation}

in which

(3.49)\begin{equation} \tilde{c}(s,t, \varLambda, \lambda) = \int_{x=0}^{\infty} {\mathcal{C}}(x,t, \varLambda, \lambda) \sin (s x) \,{{\rm d}}\kern0.7pt x. \end{equation}

This initial value problem is readily solved to obtain

(3.50)\begin{align} \tilde{c}(s,t, \varLambda, \lambda) = \int_{\xi=0}^{t} [f(\xi, \varLambda, \lambda)-g(\lambda)] R_0(\xi, \varLambda, \lambda)^{4} s \exp \left({-}s^{2} \int_{\xi}^{t} R_0(\eta, \varLambda, \lambda)^{4}\, {{\rm d}} \eta \right) {{\rm d}} \xi. \end{align}

Using the inverse Fourier sine transform, we have a solution for the concentration of gas in the liquid in the form

(3.51)\begin{align} {\mathcal{C}}(x,t, \varLambda, \lambda) &= \frac{2}{\rm \pi} \int_{s=0}^{\infty} \int_{\xi=0}^{t} [f(\xi, \varLambda, \lambda)-g(\lambda)] R_0(\xi, \varLambda, \lambda)^{4}\nonumber\\ &\quad \times \exp \left({-}s^{2} \int_{\xi}^{t} R_0(\eta, \varLambda, \lambda)^{4} \,{{\rm d}} \eta \right) {{\rm d}} \xi s \sin( s x ) \,{{\rm d}}s. \end{align}

We wish to change the order of the integration, but it is unclear whether the integral of the absolute value exists, which is a requirement to apply Fubini's theorem. Therefore, we apply integration by parts to the inner integral which yields

(3.52)\begin{align} {\mathcal{C}}(x,t, \varLambda, \lambda) &= \frac{2}{\rm \pi} \int_{s=0}^{\infty} \left \{ -\int_{\xi=0}^{t} \frac{\partial {f}}{\partial {\xi}}(\xi, \varLambda, \lambda) \exp \left({-}s^{2} \int_{\xi}^{t} R_0(\eta, \varLambda, \lambda)^{4} \,{{\rm d}} \eta \right) {{\rm d}} \xi\right. \nonumber\\ &\quad +f(t, \varLambda, \lambda)-g(\lambda) - [f(0, \varLambda, \lambda)-g(\lambda)] \nonumber\\ &\quad \times \left. \exp \left({-}s^{2} \int_{0}^{t} R_0(\eta, \varLambda, \lambda)^{4} \,{{\rm d}} \eta \right) \right \} \frac{\sin( s x )}{s} \,{{\rm d}}s. \end{align}

Fubini's theorem may now be applied to change the order of the integration. Using tables of integral transforms (Erdélyi et al. 1954), or otherwise, we evaluate the integrals in $s$ to find

(3.53)\begin{align} {\mathcal{C}}(x,t, \varLambda, \lambda) &= f(t, \varLambda, \lambda)-g(\lambda) - [f(0, \varLambda, \lambda)-g(\lambda)] \,\text{erf} \left ( \frac{x}{2} \left \{ \int_{0}^{t} R_0(\eta, \varLambda, \lambda)^{4} \,{{\rm d}} \eta \right \}^{{-}1/2} \right )\nonumber\\ &\quad -\int_{\xi=0}^{t} \frac{\partial {f}}{\partial {\xi}}(\xi, \varLambda, \lambda) \,\text{erf} \left( \frac{x}{2} \left \{ \int_{\xi}^{t} R_0(\eta, \varLambda, \lambda)^{4} \,{{\rm d}} \eta \right \}^{{-}1/2} \right) {{\rm d}} \xi, \end{align}

in which

(3.54)\begin{equation} \text{erf} (\xi) = \frac{2}{\sqrt{\rm \pi}} \int_{y=0}^{\xi} \exp({-}y^{2}) \,{{\rm d}}y. \end{equation}

It is straightforward to check that this result is an exact solution of the leading-order approximation by substituting this formula back into (3.45)–(3.47). We take the partial derivative of this analytical solution with respect to $x$ to obtain the concentration gradient at the bubble interface located at $x=0$

(3.55)\begin{align} \frac{\partial {c_0}}{\partial {x}}(0,t, \varLambda, \lambda) &={-}\frac{[f(0, \varLambda, \lambda)-g(\lambda)]}{\sqrt{\rm \pi}} \left \{ \int_{0}^{t} R_0 (\eta, \varLambda, \lambda)^{4} \,{{\rm d}} \eta \right \}^{{-}1/2} \nonumber\\ &\quad - \frac{1}{\sqrt{\rm \pi}} \int_{\xi=0}^{t} \frac{\partial {f}}{\partial {\xi}}(\xi, \varLambda, \lambda) \left \{ \int_{\xi}^{t} R_0(\eta, \varLambda, \lambda)^{4} \,{{\rm d}} \eta \right \}^{{-}1/2} {{\rm d}} \xi, \end{align}

in which

(3.56)$$\begin{gather} f(0, \varLambda, \lambda) ={-} 1 + d p_{g0} \frac{m_0(\varLambda, \lambda)^{\kappa}}{R_0(0, \varLambda, \lambda)^{3 \kappa}}, \end{gather}$$

(3.57)$$\begin{gather}\frac{\partial {f}}{\partial {\xi}} ={-} 3 \kappa d p_{g0} \frac{m_0^{\kappa}}{R_0^{3 \kappa + 1}} \frac{\partial {R_0}}{\partial {\xi}}. \end{gather}$$

Therefore, the flux of dissolved gas has been determined as an exact solution of the leading-order approximation. The first term on the right-hand side of (3.55) is the memory of the initial and matching conditions. The second is the history term (Peñas-López et al. 2016), which accounts for the preceding time history of the first derivative of the concentration at the bubble interface. We note that these memory integrals only span the shortest of our three time scales; in other words, they represent a short-term memory.

3.5. Leading-order far-field outer solution

Noticing the far field of the outer solution is at $y = \infty$, a boundary condition at a truncated far field $y = L \gg 1$ is needed to improve the accuracy and reduce the computational cost of solving (3.37) for the concentration. In this section, we will find the leading-order outer solution for $y \gg 1$. This is then exploited to derive the far-field boundary condition for the outer problem. When $y$ is sufficiently large, $R_0^{3}$ becomes negligible compared with $3y$ in (3.37). The resulting far-field partial differential equation is

(3.58)\begin{equation} \frac{\partial {C_0}}{\partial {\lambda}} = \frac{\partial}{\partial {y}} \left [ ( 3 y )^{4/3} \frac{\partial {C_0}}{\partial {y}} \right ] \quad \text{for} \ y \gg 1, \end{equation}

with boundary condition (3.38) and initial condition (3.39). The similarity solution to this problem is of the form

(3.59a,b)\begin{equation} C_0 = \phi(\eta), \quad \eta = \frac{y}{\lambda^{3/2}}. \end{equation}

The unknown function $\phi (\eta )$ satisfies the ordinary differential equation

(3.60)\begin{equation} \frac{{{\rm d}}^{2} {\phi}}{{{\rm d}} {\eta}^{2}} + \left [ \frac{4}{3 \eta} + \frac{1}{2 (3 \eta)^{1/3}} \right ] \frac{{{\rm d}} {\phi}}{{{\rm d}} {\eta}} = 0, \end{equation}

with boundary condition

(3.61)\begin{equation} \phi \rightarrow 0 \quad \text{as} \ \eta \rightarrow \infty. \end{equation}

This ordinary differential equation may be integrated to yield

(3.62)\begin{equation} \frac{{{\rm d}} {\phi}}{{{\rm d}} {\eta}} = \frac{A}{\eta^{4/3}} \exp \left ( -\frac{3^{2/3}}{4} \eta^{2/3} \right ), \end{equation}

in which $A$ is a constant of integration. We may integrate again and use the boundary condition (3.61) to find that

(3.63)\begin{equation} C_0 ={-} \int_{s=y/\lambda^{3/2}}^{\infty} \frac{A}{s^{4/3}} \exp \left ( -\frac{3^{2/3}}{4} s^{2/3} \right ) {{\rm d}}s. \end{equation}

The domain for the leading-order outer solution may now be truncated by implementing the far-field boundary condition in the form

(3.64)\begin{align} \frac{\partial {C_0}}{\partial {y}}(L, \lambda) ={-}\frac{C_0(L, \lambda) \lambda^{1/2}}{L^{4/3}} \exp \left ( -\frac{(3 L)^{2/3}}{4 \lambda} \right ) \left \{ \int_{s=L/\lambda^{3/2}}^{\infty} \frac{1}{s^{4/3}} \exp \left ( -\frac{3^{2/3}}{4} s^{2/3} \right ) {{\rm d}}s \right \}^{{-}1}, \end{align}

in which the value of $L \gg 1$ needs to be chosen such that $R_0^{3}$ is negligible in comparison with $3 L$.

The analytical solution (3.63) is used to determine the integrand in the global conservation of mass (3.44) for $y > L$, where the value of $A$ follows from the value of $C_0(L, \lambda )$.

3.6. Summary of the tractable model

For the convenience of readers, we now summarize the leading-order model to be solved numerically, in a self-contained manner. The bubble radius, $R_0 (t, \varLambda,\lambda )$, satisfies the Rayleigh–Plesset equation (3.9) on the shortest time scale $t$

(3.65)\begin{equation} R_0 \frac{\partial^{2} {R_0}}{\partial {t}^{2}} + \frac{3}{2} \left( \frac{\partial {R_0}}{\partial {t}} \right)^{2} = \frac{p_{g0} m_0^{\kappa}}{R_0^{3 \kappa}} - \frac{2}{We R_0} - 1 - \frac{4}{Re R_0} \frac{\partial {R_0}}{\partial {t}} - p_a \sin (2 {\rm \pi}\nu_a t), \end{equation}

in which $Re$ is the Reynolds number, $We$ is the Weber number, $p_{g0}$ is the dimensionless initial partial pressure of the bubble gases and $p_a$ ($\nu _a$) is the dimensionless amplitude (frequency) of the applied external pressure. The mass of gas in the bubble, $m_0 (\varLambda,\lambda )$, is governed by conservation of mass (3.25) on the intermediate time scale $\varLambda$

(3.66)\begin{equation} \frac{\partial {m_0}}{\partial {\varLambda}} = 4 {\rm \pi}M \left \langle R_0^{4} \frac{\partial {c_0}}{\partial {x}} (0,t,\varLambda,\lambda) \right \rangle, \end{equation}

in which the concentration gradient is determined from (3.55) as follows:

(3.67)\begin{align} \frac{\partial {c_0}}{\partial {x}}(0,t, \varLambda, \lambda) &={-}\frac{[f(0, \varLambda, \lambda)-g(\lambda)]}{\sqrt{\rm \pi}} \left \{ \int_{0}^{t} R_0 (\eta, \varLambda, \lambda)^{4} \,{{\rm d}} \eta \right \}^{{-}1/2} \nonumber\\ &\quad - \frac{1}{\sqrt{\rm \pi}} \int_{\xi=0}^{t} \frac{\partial {f}}{\partial {\xi}}(\xi, \varLambda, \lambda) \left \{ \int_{\xi}^{t} R_0(\eta, \varLambda, \lambda)^{4} \,{{\rm d}} \eta \right \}^{{-}1/2} {{\rm d}} \xi, \end{align}

where

(3.68)$$\begin{gather} \langle \,\cdot\, \rangle = \frac{1}{P} \int^{t+P}_{t} {\cdot} \,{{\rm d}}t, \end{gather}$$

(3.69)$$\begin{gather}f(0, \varLambda, \lambda) ={-} 1 + d p_{g0} \frac{m_0(\varLambda, \lambda)^{\kappa}}{R_0(0, \varLambda, \lambda)^{3 \kappa}}, \end{gather}$$

(3.70)$$\begin{gather}\frac{\partial {f}}{\partial {\xi}} ={-} 3 \kappa d p_{g0} \frac{m_0^{\kappa}}{R_0^{3 \kappa + 1}} \frac{\partial {R_0}}{\partial {\xi}}, \end{gather}$$

and the period of bubble oscillation is $P(\varLambda, \lambda )$. The concentration of gas in the outer liquid region, $C_0 (y, \lambda )$, is determined by the diffusion equation (3.37) on the longest time scale $\lambda$

(3.71)\begin{equation} \frac{\partial {C_0}}{\partial {\lambda}} = \frac{\partial}{\partial {y}} \left [ \langle ( R_0^{3} + 3 y )^{4/3} \rangle \frac{\partial {C_0}}{\partial {y}} \right ] \quad \text{for} \ 0< y < L, \end{equation}

with matching condition $C_0(0, \lambda ) = g(\lambda )$ and far-field boundary condition from (3.64)

(3.72)\begin{align} \frac{\partial {C_0}}{\partial {y}}(L, \lambda) ={-}\frac{C_0(L, \lambda) \lambda^{1/2}}{L^{4/3}} \exp \left ( -\frac{(3 L)^{2/3}}{4 \lambda} \right ) \left \{ \int_{s=L/\lambda^{3/2}}^{\infty} \frac{1}{s^{4/3}} \exp \left ( -\frac{3^{2/3}}{4} s^{2/3} \right ) {{\rm d}}s \right \}^{{-}1} \end{align}

to truncate the domain at $L \gg 1$. Equation (3.44),

(3.73)\begin{equation} 4 {\rm \pi}M \int_{y=0}^{\infty} C_0 \, {{\rm d}}y = 1-m_0, \end{equation}

determines $g(\lambda )$ in the matching condition $C_0(0, \lambda ) = g(\lambda )$, where $C_0$ is evaluated on the semi-infinite interval $[L, \infty )$ via the analytical solution (3.63)

(3.74)\begin{equation} C_0 ={-} \int_{s=y/\lambda^{3/2}}^{\infty} \frac{A}{s^{4/3}} \exp \left ( -\frac{3^{2/3}}{4} s^{2/3} \right ) {{\rm d}}s, \end{equation}

in which the value of $A$ is taken to be

(3.75)\begin{equation} A ={-}C_0(L, \lambda) \left \{ \int_{s=L/\lambda^{3/2}}^{\infty} \frac{1}{s^{4/3}} \exp \left ( -\frac{3^{2/3}}{4} s^{2/3} \right ) {{\rm d}}s \right \}^{{-}1}. \end{equation}

The initial conditions are given by

(3.76a–d)\begin{equation} R_0(0,0,0) = 1, \quad \frac{\partial {R_0}}{\partial {t}}(0,0,0) = 0, \quad C_0 ( y, 0 ) = 0, \quad m_0(0,0) = 1, \end{equation}

for $0< y < L$.

4. Numerical method

4.1. Algorithm

The computational algorithm includes the integration of the Rayleigh–Plesset equation (3.9) on the shortest time scale $t$, the integration of the conservation of mass (3.25) on the intermediate time scale $\varLambda$ and solving the diffusion equation (3.37) on the longest time $\lambda$. The NAG routine D02EJF is employed to solve the Rayleigh–Plesset equation (3.9) on the shortest time scale (Numerical Algorithms Group 2019b). This is the only one of our equations on this shortest time scale which requires billions of time steps to simulate the many millions of bubble oscillations.

The mid-point rule is utilized to integrate the first derivative with respect to $\varLambda$ in (3.25) for the mass of gas in the bubble. The time integrals on the right-hand side of (3.25) and in the expression for the concentration gradient (3.55) are evaluated with the trapezoidal composite rule. The mass of gas in the bubble is updated on the completion of each cycle of bubble oscillation.

The method of lines reduces the diffusion equation (3.37) to a system of ordinary differential equations (Smith 1985) which may then be solved with the NAG routine D02EJF. The domain is chosen such that $L=50$. The concentration of gas in the outer liquid region is updated on the completion of each cycle of bubble oscillation.

The trapezoidal composite rule is applied to evaluate the spatial integral in the matching condition (3.44) in the interval $[0,L)$ (Burden & Faires 1985). The remainder of the integral in (3.44) on the semi-infinite interval $[L,\infty )$ uses the analytical solution (3.63) and the NAG routine D01AMF (Numerical Algorithms Group 2019a). The constant $A$ in the analytical solution (3.63) is evaluated from the value of $C_0$ at $y=L$. The matching condition $g(\lambda )$ is updated on the completion of each cycle of bubble oscillation.

The simulations were undertaken for the growth of bubbles over several million oscillations. The calculations were performed on the supercomputer BlueBEAR (University of Birmingham 2022). The basic mathematical model definitely cannot be traced with any CFD packages or known numerical methods, because it requires accurately solving the partial differential equation involved for billions of time steps in order to simulate millions of cycles of oscillation.

4.2. Comparison with experiment

To evaluate the present theoretical model, we will compare the computational results with two sets of experiments. Direct comparisons are undertaken for the maximum and equilibrium bubble radii during each cycle of oscillation. The growth and dissolution of bubbles subject to acoustic waves over millions of cycles are considered. This displays the excellent agreement of the theory with experiments.

The first set of experiments considered is for air-saturated water in a single bubble levitation cell (Lee et al. 2005). A hollow cylindrical lead zirconate titanate piezo ceramic transducer at the bottom of the cell was used to excite a standing acoustic wave at $22.35$ kHz and a strobe was used to capture the maximum bubble radius during its oscillations. A syringe was used to introduce a single bubble. The following values for gas bubbles in water are adopted following the experimental conditions $\Delta = 10^{5}$ kg m$^{-1}$s$^{-2}$, $\sigma = 0.072$ Nm$^{-1}$, $\kappa = 1.4$, $\mu = 0.001$ Pa s, $\rho = 10^{3}$ kg m$^{-3}$, $\bar {\nu }_{a} = 22.35 \times 10^{3}$ Hz, $\bar {p}_a = 2.2 \times 10^{4}$ kg m$^{-1}$ s$^{-2}$, $D = 2 \times 10^{-9}$ m$^{2}$ s$^{-1}$ and $k_{H} = 2.5 \times 10^{6}$ m$^{2}$ s$^{-2}$. The initial pressure of the bubble gas is given by $\bar {p}_{g0} = \Delta + 2 \sigma / \bar {R}_{eqi}$, being in the equilibrium state.

The initial mass of gas in the bubble is given by

\[ \bar{m}_i = \phi_{air} \rho_{air} \frac{4}{3} {\rm \pi}\bar{R}_{eqi}^{3}, \]

in which $\rho _{air}$ is the density and $\phi _{air}$ is the initial volume fraction of air in the bubble. As the solubility of oxygen in water is higher than the solubility of nitrogen, the air in water is assumed to be $36$% oxygen so that $\rho _{air}=1.3$ kg m$^{-3}$ (Engineering ToolBox 2004). The saturation concentration of the gas in the liquid, $\bar {C}_{sat} = 4.057 \times 10^{-2}$ kg m$^{-3}$, has been chosen to fit the experimental observations of bubble growth (Lee et al. 2005) for the initial maximum bubble radius of approximately $22$ $\mathrm {\mu }$m. We choose the initial uniform concentration of the gas in the liquid $\bar {C}_i = \bar {C}_{sat}$ (Lee et al. 2005). The only unknown parameter value is the initial volume fraction of air in the bubble which would almost certainly take different values in experiments.

In order to validate our simplified model, a numerical solution was obtained for the above parameter values, two values of the initial volume fraction of air in the bubble ($\phi _{air}=$ $0.8$ and $0.9$) and an initial equilibrium bubble radius of $\bar {R}_{eqi} = 32.5$ $\mathrm {\mu }$m. The initial equilibrium bubble radius of $\bar {R}_{eqi} = 32.5$ $\mathrm {\mu }$m corresponds to the initial maximum bubble radius of approximately $35$ $\mathrm {\mu }$m in the experimental results of Lee et al. (2005). We have $Pe \approx 1.6 \times 10^{5}$, which is very large. Lee et al. (2005) estimated that the bubble remains spherical for the first $150$ s of its growth, so we only simulate for this time period. Figure 3 compares two numerical simulations for the time history of the maximum bubble radius during each cycle of oscillation with experimental results. The computational results agree with the experimental results for the long time history, for both values of the initial volume fraction. The bubble grows nonlinearly with time for the period considered and the growth rate decreases with increases in the initial volume fraction of air in the bubble.

Figure 3. Comparison of the time histories of the maximum bubble radius, $\bar {R}_{max}$, using two numerical solutions of the simplified model in subsection 4.1 and experimental results (Lee et al. 2005). In this experiment, a bubble grows in air-saturated water subject to a standing acoustic wave at $22.35$ kHz. The numerical solution is shown for two values of the initial volume fraction of air in the bubble: $\phi _{air}=$ $0.8$ and $0.9$. Other parameter values are $Re \approx 330$, $We \approx 45$, $p_{g0} \approx 1.0$, $p_a \approx 0.22$, $\nu _a \approx 0.073$, $Pe \approx 1.6 \times 10^{5}$, $d \approx 0.99$, $M \approx 9.2 \times 10^{-3}$ and $\kappa =1.4$.

We note that, if we had simulated over a longer time period of $210$ s, then the bubble typically exhibits instability (variations between no oscillations and very large oscillations). This instability occurs as the equilibrium bubble radius approaches approximately $50$ $\mathrm {\mu }$m. Lamb (1932) determined the natural frequency $\bar {f}_n$ of shape mode $n$ to be given by

(4.1)\begin{equation} \bar{f}_n = \frac{1}{2 {\rm \pi}\bar{R}_{eq}} \sqrt{(n-1)(n+1)(n+2)\frac{\sigma}{\rho \bar{R}_{eq}}}, \end{equation}

for $n \geq 2$ (see also figure 1). Using (4.1), for a frequency of $22.35$ kHz, the resonant radius of shape mode $n=3$ occurs at $53$ $\mathrm {\mu }$m. The transition to non-spherical oscillations in figure 1 of Lee et al. (2005) may be discerned to be shape mode $n=3$, albeit for an aqueous surfactant solution. There may be a correspondence between our instability at approximately $50$ $\mathrm {\mu }$m and the transition to non-spherical oscillations in experiments.

The dimensional equilibrium bubble radius, $\bar {R}_{eq}$, is evaluated as the solution to the following algebraic equation

\[ \bar{p}_{g0} \left ( \frac{\bar{m}}{\bar{m}_i} \right )^{\kappa} \left( \frac{\bar{R}_{eqi}}{\bar{R}_{eq}} \right)^{3 \kappa} - \frac{2 \sigma}{\bar{R}_{eq}} -(\bar{p}_{\infty} - \bar{p}_v)=0. \]

Figure 3 also shows an oscillation at a maximum bubble radius of approximately $40$ $\mathrm {\mu }$m. The maximum bubble radius is not monotonically increasing. Figure 4 shows that this oscillation does not appear in the time history of the equilibrium bubble radius or the mass of gas in the bubble. The oscillation corresponds to the dynamics of the bubble and not its underlying growth. The experimental observations of Lee et al. (2005) exhibit a range of values at a maximum bubble radius of approximately $40$ $\mathrm {\mu }$m, but these occur at earlier times than the oscillation which we have discovered in the theoretical results. Figures 3 and 4 demonstrate that the maximum bubble radius represents changes in the underlying growth except when the maximum radius is approximately $40$ $\mathrm {\mu }$m. The growth rate is slowed significantly around 80–100 s for $\phi _{air}= 0.8$ and 95–115 s for $\phi _{air}= 0.9$. We will investigate these two features by studying the bubble oscillation over short periods.

Figure 4. The time histories of (a) the equilibrium bubble radius $\bar {R}_{eq}$ and (b) the dimensionless mass of gas in the bubble $m_0$. The numerical solutions are shown for two values of the initial volume fraction of air in the bubble corresponding to figure 3.

Approximately three million oscillations are simulated over the time period of $150$ s in figures 3 and 4. Figure 5 shows six short time periods of $10^{-4}$ s, three before and three after the oscillation in the maximum bubble radius in figure 3 for $\phi _{air}=0.8$. The periods of oscillation in figure 5 agree with the period of the applied external pressure of approximately $4.5 \times 10^{-5}$ s. A second common feature is that the bubble splits its time approximately equally between greater and less than the equilibrium radius except for figure 5(c). A more intriguing common feature in figure 5 is that the maxima are further above the equilibrium bubble radius than the minima are below. This third common feature supports the area and shell mechanisms discussed in the introduction. In order to quantify this feature, we introduce the ratio

\[ \delta = \frac{\bar{R}_{max}-\bar{R}_{eq}}{\bar{R}_{eq}-\bar{R}_{min}}, \]

where $\bar {R}_{min}$ in the dimensional minimum bubble radius. In figure 5( f), $\delta = 1.39$ at $\bar {t}=150$ s.

Figure 5. Short time periods of $10^{-4}$ s of the histories of the bubble radius $\bar {R}$ for (a) $\bar {t}=30$ s, (b) $\bar {t}=60$ s, (c) $\bar {t}=95$ s before the oscillation at approximately $40$ $\mathrm {\mu }$m, (d) $\bar {t}=105$ s after the oscillation at approximately $40$ $\mathrm {\mu }$m, (e) $\bar {t}=120$ s and ( f) $\bar {t}=150$ s. The equilibrium bubble radius is also shown for comparison. The volume fraction of the gas in the bubble $\phi _{air}=0.8$ and the remaining parameter values are the same as figure 3.

Figure 5(d) corresponds to a much greater bubble growth than figure 5(c). The ratio in figure 5(c) of $\delta =1.44$ is slightly greater than in figure 5(d) of $\delta =1.3$. However, figure 5(c) has a large high frequency mode (bubble resonance) which is not present in figure 5(d). The disappearance of this high frequency mode is responsible for the increased bubble growth rate at $\bar {t}=105$ s in comparison with $\bar {t}=95$ s.

The natural frequency for spherical oscillation of a bubble, $\bar {f}_0$, is given by Brennen (2013)

(4.2)\begin{equation} \bar{f}_0 = \frac{1}{2 {\rm \pi}\rho \bar{R}_{eq}} \sqrt{\rho \left[ 3 \kappa (\bar{p}_{\infty} - \bar{p}_v) + \frac{2 (3 \kappa - 1) \sigma}{\bar{R}_{eq}} \right ] - \frac{4 \mu^{2}}{\bar{R}_{eq}^{2}}}. \end{equation}

Using (4.2), the natural frequency corresponding to figure 5(c) is approximately $90$ kHz for $\bar {R}_{eq}=36.7$ $\mathrm {\mu }$m. The bubble resonance frequency is four times the acoustic frequency in figure 5(c). This drives the natural oscillation of the bubble, as observed in the figure.

The oscillation at a maximum bubble radius of approximately $40$ $\mathrm {\mu }$m occurs later in our theoretical results than in the experimental observations doubtless due to the proximity of the bubble resonance. As discussed in Appendix C of Fyrillas & Szeri (1994), the polytropic bubble model is inaccurate near to the bubble resonance which would be expected to influence the accuracy of the maximum bubble radius.

Figures 5(a) and 5(b) examine the evolution of the bubble oscillations from the initiation of the acoustic wave until the oscillation in the maximum bubble radius. The amplitude of the bubble oscillation ($\bar {R}_{max}-\bar {R}_{min}$) increases from less than 4 $\mathrm {\mu }$m in figure 5(a) to more than 5 $\mathrm {\mu }$m in figure 5(c). The high frequency mode is not visible in figure 5(a), but its growth is already evident in figure 5(b). It is noteworthy that the maxima are sharp in figure 5(a) whilst the minima are flat. In the times following the oscillation in the maximum bubble radius shown in figures 5(e) and 5( f), the bubble oscillation maintains the temporal shape in figure 5(d) but its amplitude increases.

Figure 6 shows the variation of the concentration of gas in the outer region for four different times (with the numerical solution being extended using the far-field analytical solution (3.63)). The spatial variation starts near to the bubble interface and extends to the constant far-field value for each time. The extent of the outer region increases as time progresses, which is consistent with the process of diffusion. As the dissolved gas enters the bubble, the gas concentration decreases at the bubble interface. The dissolved gas in the liquid diffuses towards the bubble interface. This results in a monotonic spatial decrease in the gas concentration towards the bubble interface. The spatial gradient of the gas concentration increases rapidly near the bubble interface. However, the mass flux from the outer region into the inner region does not vary monotonically in time. It has a local minimum at approximately $\bar {t}=95$ s corresponding to bubble resonance and a local maximum at approximately $\bar {t}=105$ s. The concentration of the gas in the liquid does not change significantly if the radius is larger than $0.0008$ m. The extent of the outer region is an important factor limiting the application of results on a single bubble to a cloud of bubbles. The distance between a bubble and its neighbouring bubbles should be less than the sum of the spatial lengths of their respective outer regions.

Figure 6. The concentration of gas in the outer region $\bar {c}$ at $\bar {t}=50$ s, $\bar {t}=95$ s, $\bar {t}=105$ s and $\bar {t}=150$ s. The bubble interface corresponds to the equilibrium bubble radius. The numerical solutions are shown for an initial volume fraction of air in the bubble $\phi _{air}=0.8$ and the remaining parameter values are the same as figure 3.

The second experimental case to be considered is air-saturated water in a hollow glass cylinder between two hollow cylindrical transducers (Crum 1980). Bubbles were produced by increasing the gain on the amplifier until the threshold for gaseous cavitation was reached. The desired initial radii were obtained in a few minutes by adjustment of the amplitude of the acoustic pressure. These bubbles were levitated near the anti-node of an acoustic stationary wave at a frequency of $22.1$ kHz. The equilibrium bubble radius was deduced by measuring its terminal velocity through the liquid. The following values for gas bubbles in water are adopted following the experimental conditions $\Delta = 10^{5}$ kg m$^{-1}$ s$^{-2}$, $\sigma = 0.068$ Nm$^{-1}$, $\kappa = 1.4$, $\mu = 0.001$ Pa s, $\rho = 10^{3}$ kg m$^{-3}$, $\bar {\nu }_{a} = 22.1 \times 10^{3}$ Hz, $\bar {p}_a = 2.7 \times 10^{4}$ kg m$^{-1}$ s$^{-2}$, $D = 2.4 \times 10^{-9}$ m$^{2}$ s$^{-1}$ and $k_{H} = 2.5 \times 10^{6}$ m$^{2}$ s$^{-2}$. The initial pressure of the bubble gas $\bar {p}_{g0}$ and the initial mass of gas in the bubble $\bar {m}_i$ are given by the same formulas as above. The initial uniform concentration of the gas in the liquid, $\bar {C}_i = \bar {C}_{sat} = 3.888 \times 10^{-2}$ kg m$^{-3}$, has been chosen to fit the observations for the initial maximum bubble radius of approximately $30$ $\mathrm {\mu }$m. In practice, it is very difficult to locate the exact value of the saturation concentration of the gas in the liquid. We choose the initial volume fraction of air in the bubble $\phi _{air}$ to achieve agreement with the experimental results. We have $Pe \approx 1.5 \times 10^{5}$, which is very large. Due to the invasive nature of the measurement of the equilibrium bubble radius, the concentration of gas in the liquid is reset to the initial uniform concentration each time that the experimental observations are provided in Crum (1980).

Further validation of our simplified model is sought via a comparison of the time histories of the equilibrium bubble radius with the experimental results of Crum (1980). Three representative initial equilibrium bubble radii $\bar {R}_{eqi}=35.6$, 24.7 or 29.4 $\mathrm {\mu }$m are considered, for which the bubble grows, dissolves or remains at a steady state, respectively, in the experiments. Our computations repeated these three processes remarkably, both qualitatively and quantitatively. The bubble growth is simulated for a time period of $160$ s (approximately three and a half million oscillations) in figure 7. After $160$ s, the equilibrium bubble radius exhibits instability, which may correspond to a transition to non-spherical oscillations. Figure 7 also shows bubble dissolution for a time period of $160$ s. The equilibrium bubble radius decreases faster as time increases in both the numerical simulation and the experiment. However, the numerical values decrease slightly more rapidly than the experimental observations.

Figure 7. Comparison of the time histories of the equilibrium bubble radius, $\bar {R}_{eq}$, using three numerical solutions of the simplified model in subsection 4.1 and experimental results (Crum 1980). In this experiment, a bubble grows, dissolves or remains at a steady state depending on the initial equilibrium bubble radius $\bar {R}_{eqi}=35.6$, 24.7 or 29.4 $\mathrm {\mu }$m, respectively. The dissolving bubble uses an initial volume fraction of air in the bubble $\phi _{air}=0.99$ and the other two simulations use $\phi _{air}=0.9$. Other parameter values are $Re \approx 360$, $We \approx 52$, $p_{g0} \approx 1.0$, $p_a \approx 0.27$, $\nu _a \approx 0.079$, $Pe \approx 1.5 \times 10^{5}$, $d \approx 1.0$, $M \approx 7.9 \times 10^{-3}$ and $\kappa =1.4$.

Figure 8 shows the variation of the concentration of gas in the outer region at times corresponding to the experimental observations of Crum (1980) for both growth and dissolution at approximately $40$ s in figure 7. The mass flux from the inner region into the outer region for bubble dissolution is greater than the mass flux from the outer region into the inner region for bubble growth. We also note that the spatial length of the outer region is similar for the growth and dissolution of bubbles.

Figure 8. The concentration of gas in the outer region $\bar {c}$ at $\bar {t}=43.1$ s for bubble growth and at $\bar {t}=41.0$ s for bubble dissolution. The bubble interface corresponds to the equilibrium bubble radius. The parameter values are the same as figure 7.

4.3. Convergence tests

In order to evaluate the convergence of the method of lines, tests were performed for the volume fraction of the gas in the bubble $\phi _{air}=0.8$ and the remaining parameter values as in figure 3. The growth of bubbles was simulated with $250$, $500$ and $1000$ mesh points. Figure 9(a) shows the results for the spatial distribution of the dimensionless concentration of gas in the outer region at $50$ s. Figure 9(b–d) shows the times histories of the maximum bubble radius, equilibrium bubble radius and dimensionless mass of gas in the bubble, respectively. These results converge well for a number of mesh points greater than or equal to $250$. All the other calculations in this article have been performed with $500$ mesh points.

Figure 9. Convergence tests for the mesh in the outer region: (a) the dimensionless concentration of gas in the outer region $C_0$ at $\bar {t}=50$ s with $0< y< L$, (b) the time history of the maximum bubble radius $\bar {R}_{max}$, (c) the time history of the equilibrium bubble radius $\bar {R}_{eq}$ and (d) the time history of the dimensionless mass of gas in the bubble $m_0$. The volume fraction of the gas in the bubble $\phi _{air}=0.8$ and the remaining parameter values are the same as figure 3 with number of mesh points $250$, $500$ and $1000$.

5. Parametric studies

In each of the following subsections, the effects of varying a single parameter will be examined. The parameters will be the applied pressure frequency, the applied pressure amplitude and the initial gas concentration in the liquid. Our reference will be the parameter values used in figure 3 with an initial volume fraction of the gas in the bubble $\phi _{air}=0.8$.

We must ensure that thermal effects remain negligible for our adopted parameter values. The temperature dependence of the gas solubility, the vapour pressure, the diffusivity and the surface tension may have significant influence on the growth of bubbles by rectified diffusion (Webb et al. 2010). The repercussion of not including thermal effects is that any theoretical model may only integrate over time periods such that the temperature rise is insignificant. The thermal source term for ultrasonic absorption varies quadratically with the applied pressure amplitude (Pierce 1981) and quadratically with the applied pressure frequency via the absorption coefficient (Liebermann 1948). Therefore, we examine a reduction in the applied pressure amplitude and a moderate increase in the applied pressure frequency to guarantee that the temperature does not rise by more than $0.4^{\circ }{\rm C}$, which was observed by Lee et al. (2005).

5.1. Effects of applied pressure frequency

In addition to the simulation at $\bar {\nu }_a=22.35$ kHz in figure 3, results are obtained at the applied pressure frequencies of $\bar {\nu }_a=62$ kHz, $\bar {\nu }_a=42$ kHz, $\bar {\nu }_a=17$ kHz and $\bar {\nu }_a=12$ kHz. The time histories of the equilibrium bubble radius are compared in figure 10. The bubble growth rate increases with the applied pressure frequency in this range of values. At the frequency of $\bar {\nu }_a=17$ kHz, the high frequency mode (or bubble resonance) appears later and at larger equilibrium bubble radius in comparison with the $\bar {\nu }_a=22.35$ kHz simulation (see also figure 11c), the time scales for the growth and elimination of the high frequency mode being comparable. Using (4.2), the natural frequency corresponding to figure 11(c) is approximately $85$ kHz for $\bar {R}_{eq}=38.8$ $\mathrm {\mu }$m. The bubble resonance frequency is five times the acoustic frequency in figure 11(c). The high frequency mode for the $\bar {\nu }_a=12$ kHz case has not appeared in the first $150$ s.

Figure 10. Comparison of the time histories of the equilibrium bubble radius, $\bar {R}_{eq}$, for five values of the applied pressure frequency $\bar {\nu }_a$. The initial volume fraction of the gas in the bubble $\phi _{air}=0.8$ and the remaining parameter values, with the exception of $\bar {\nu }_a$, are the same as figure 3.

Figure 11. Short time periods of $10^{-4}$ s of the histories of the bubble radius $\bar {R}$ for (a) $\bar {t}=30$ s, (b) $\bar {t}=60$ s, (c) $\bar {t}=120$ s and (d) $\bar {t}=150$ s. The equilibrium bubble radius is also shown for comparison. The applied pressure frequency $\bar {\nu }_a=17$ kHz, the initial volume fraction of the gas in the bubble $\phi _{air}=0.8$ and the remaining parameter values are the same as figure 3.

For the frequencies above $\bar {\nu }_a=22.35$ kHz, the amplitude of oscillation increases and, as a result, the effects of nonlinearity are more evident in the time histories of the equilibrium bubble radius. At the frequency $\bar {\nu }_a=42$ kHz, the high frequency mode (or bubble resonance) occurs much earlier than at the frequency $\bar {\nu }_a=22.35$ kHz, but at a similar equilibrium bubble radius. After the usual rapid increase following the bubble resonance, there is a small decrease in equilibrium bubble radius, before the new, faster growth rate. The nonlinear dynamics of the bubble radius must be responsible for this short period of dissolution. At the frequency $\bar {\nu }_a=62$ kHz, there are no bubble resonances. There is an initial period of rapid growth followed by a decrease in the equilibrium bubble radius and then a longer time period of slower growth.

Approximately two and a half million oscillations are simulated over the time period of $150$ s for $\bar {\nu }_a=17$ kHz. Figure 11 shows four short time periods of $10^{-4}$ s, three before and one after the resonant oscillation. The periods of oscillation in figure 11 agree with the period of the applied external pressure of approximately $5.9 \times 10^{-5}$ s. The results in figure 11 share the common features of figure 5; however, there are two notable differences. Firstly, the minima in figure 11(a,b) are not flat as in figure 5(a). Secondly, the minima become flat in figure 11(d), whereas the minima are sharp in figure 5(d–f). The former doubtless contributes to the ratio $\delta =1.37$ for figure 11(a) which is smaller than $\delta =1.45$ for figure 5(a), whilst the latter contributes to the ratio $\delta =1.41$ for figure 11(d) which is greater than $\delta =1.39$ for figure 5( f). Therefore, the area and shell mechanisms are responsible for the smaller growth rate in figure 10 for $\bar {\nu }_a=17$ kHz in comparison with $\bar {\nu }_a=22.35$ kHz at $\bar {t}=30$ s and the larger growth rate in figure 10 for $\bar {\nu }_a=17$ kHz in comparison with $\bar {\nu }_a=22.35$ kHz at $\bar {t}=150$ s.

Louisnard & Gomez (2003) state that the growth rate scales as the reciprocal of the frequency which is the opposite trend to the one which we have found in our simulations. They consider acoustic amplitudes in the range between $1$ and $5$ bars, whereas our results are for $0.22$ bar. Our results, before the shape mode and transient cavitation thresholds, are qualitatively different, but these results are also outside the range studied by Louisnard & Gomez (2003). There is no contradiction between the two conclusions. Further experiments over the range of applied pressure frequencies studied in these theoretical models would be interesting future work.

5.2. Effects of applied pressure amplitude

In addition to the simulation at $\bar {p}_a=0.22$ bar in figure 3, results are obtained at the applied pressure amplitudes of $\bar {p}_a=0.17$ bar and $\bar {p}_a=0.12$ bar. The time histories of the equilibrium bubble radius are compared in figure 12. The bubble growth rate increases with the applied pressure amplitude, the growth rate being considerably more sensitive to changes in the applied pressure amplitude than the applied pressure frequency. For $\bar {p}_a=0.22$ bar and $\bar {p}_a=0.17$ bar, the bubble grows; whereas, for $\bar {p}_a=0.12$ bar, the bubble dissolves. There is no oscillation in the maximum bubble radius (or bubble resonance) in the first $150$ s for either $\bar {p}_a=0.17$ bar or $\bar {p}_a=0.12$ bar.

Figure 12. Comparison of the time histories of the equilibrium bubble radius, $\bar {R}_{eq}$, for three values of the applied pressure amplitude $\bar {p}_a$. The initial volume fraction of the gas in the bubble $\phi _{air}=0.8$ and the remaining parameter values, with the exception of $\bar {p}_a$, are the same as figure 3.

Approximately three million oscillations are simulated over the time period of $150$ s for $\bar {p}_a=0.12$ bar. Figure 13 shows four short time periods of $10^{-4}$ s. The oscillations in figure 13 are quite unlike the oscillations in figure 5: the maxima and minima have similar temporal shapes at every time in figure 13. The values of the ratio $\delta$ in figure 13 corresponding to bubble dissolution are clearly less than those in figure 5 corresponding to bubble growth. In figure 13(d), the value of the ratio is $\delta =1.18$, which is much less than $\delta =1.39$ in figure 5( f). The bubble shrinks due to the area and shell mechanisms being no longer sufficient to counteract the natural process of bubble dissolution (Blake 1949).

Figure 13. Short time periods of $10^{-4}$ s of the histories of the bubble radius $\bar {R}$ for (a) $\bar {t}=30$ s, (b) $\bar {t}=60$ s, (c) $\bar {t}=120$ s and (d) $\bar {t}=150$ s. The equilibrium bubble radius is also shown for comparison. The applied pressure amplitude $\bar {p}_a=0.12$ bar, the initial volume fraction of the gas in the bubble $\phi _{air}=0.8$ and the remaining parameter values are the same as figure 3.

5.3. Effects of initial gas concentration

In addition to the simulation at $\bar {C}_i/\bar {C}_{sat}=1$ in figure 3, results are obtained at the initial uniform concentration of the gas in the liquid $\bar {C}_i/\bar {C}_{sat}=0.975$ and $\bar {C}_i/\bar {C}_{sat}=0.95$. The time histories of the equilibrium bubble radius are compared in figure 14. The bubble growth rate increases with the initial uniform concentration of the gas in the liquid, the growth rate being extremely sensitive to changes as small as one per cent. For $\bar {C}_i/\bar {C}_{sat}=1$, the bubble grows; whereas, for $\bar {C}_i/\bar {C}_{sat}=0.95$, the bubble rapidly dissolves. In the case of $\bar {C}_i/\bar {C}_{sat}=0.975$, the bubble grows very slowly indeed. The crucial factor appears to be the ratio of two pressures $\bar {p}_{g0}$ and $k_H \bar {C}_i$. We note that Church (1988) found that the thresholds of rectified diffusion are quite sensitive to changes in the initial gas concentration.

Figure 14. Comparison of the time histories of the equilibrium bubble radius, $\bar {R}_{eq}$, for three values of the initial uniform concentration of the gas in the liquid divided by the saturation concentration of the gas in the liquid $\bar {C}_i/\bar {C}_{sat}$. The initial volume fraction of the gas in the bubble $\phi _{air}=0.8$ and the remaining parameter values, with the exception of $\bar {C}_i$, are the same as figure 3.

Figure 14 also shows an oscillation in the maximum bubble radius (or bubble resonance) for $\bar {C}_i/\bar {C}_{sat}=0.95$ at approximately $\bar {t}=60$ s (see also the high frequency mode in figure 15b). This oscillation in the maximum bubble radius occurs long before the one in the case $\bar {C}_i/\bar {C}_{sat}=1$. Using (4.2), the natural frequency corresponding to figure 15(b) is approximately $112$ kHz for $\bar {R}_{eq}=29.6$ $\mathrm {\mu }$m. The bubble resonance frequency is five times the acoustic frequency in figure 15(b).

Figure 15. Short time periods of $10^{-4}$ s of the histories of the bubble radius $\bar {R}$ for (a) $\bar {t}=30$ s, (b) $\bar {t}=60$ s, (c) $\bar {t}=120$ s and (d) $\bar {t}=150$ s. The equilibrium bubble radius is also shown for comparison. The initial uniform concentration of the gas in the liquid $\bar {C}_i/\bar {C}_{sat}=0.95$, the initial volume fraction of the gas in the bubble $\phi _{air}=0.8$ and the remaining parameter values are the same as figure 3.

Figure 15 shows four short time periods of $10^{-4}$ s for $\bar {C}_i/\bar {C}_{sat}=0.95$. The oscillations in figure 15 have a number of similarities with the oscillations in figure 5: the sharp maxima and flat minima prior to the resonant oscillation in the maximum bubble radius, the growth of a high frequency mode and the bubble oscillation maintaining its temporal shape in the times following the oscillation in the maximum bubble radius. However, in figure 15(c,d), the oscillation amplitude ($\bar {R}_{max}-\bar {R}_{min}$) decreases, which is consistent with the steadily increasing bubble dissolution rate in figure 14. We recall that large nonlinear oscillations are required to counter the natural process of bubble dissolution (Blake 1949).

6. Summary and conclusions

Rectified diffusion is associated with wide, important and countless applications; however, a quantitative understanding of its development has remained largely elusive for over a century. A theoretical study has been undertaken to investigate the growth and dissolution of spherical bubbles subject to acoustic forcing. This study reveals the temporal and spatial structure of various physical processes involved. This is a multi-scaled problem with the shortest time scale, $t=\tau$, associated with bubble oscillation, the intermediate scale, $\varLambda = Pe^{-1/2} \tau$, associated with the variation of the mass of gas in the bubble and the longest time scale, $\lambda = Pe^{-1} \tau$, associated with the diffusion of gas in the liquid. Here, $Pe$ is the Péclet number for mass transfer and $\tau$ is dimensionless time. There are also two spatial length scales: the short scale, $Pe^{-1/2} \bar {R}_{eqi}$, adjacent to the bubble, and the long scale, $\bar {R}_{eqi}$, for the outer liquid region, where $\bar {R}_{eqi}$ is the initial equilibrium bubble radius.

A multi-scaled perturbation method is thus employed in time and matched asymptotic expansions in space. The liquid surrounding the spherical bubble is divided into an inner region adjacent to the bubble and an outer region. In the inner region, a Lagrangian change of variables, the Fourier sine transform and Fubini's theorem are applied to determine an exact solution for the leading-order concentration of gas in the liquid. Using this exact solution, the basic mathematical model is simplified by the removal of the convection–diffusion equation on the shortest time scale. In the outer region, the concentration of gas in the liquid is shown to vary on the longest time scale. A far-field similarity solution has been found for the leading-order concentration of gas in the liquid, so that the diffusion equation may be numerically integrated on a truncated domain. Asymptotic matching is achieved by exploiting a global conservation of mass. Our techniques result in a tractable and accurate mathematical model for rectified diffusion. It consists of the Rayleigh–Plesset equation for the spherical oscillations of the bubble on the shortest time scale; an ordinary differential equation for the mass of gas in the bubble on the intermediate time scale; and a diffusion equation for the concentration of gas in the liquid on the longest time scale and the long length scale. As the Péclet number for mass transfer is usually very large, asymptotic errors are small and this simplified model is essentially equivalent to the basic mathematical model.

Predictions of our simplified model agree with experimental results for a single spherical bubble in the bulk over millions of cycles of oscillation. Theoretical and numerical studies are undertaken for the growth rate of the bubble, the bubble oscillations and the concentration of gas in the liquid in terms of the applied pressure frequency (12–62 kHz), the applied pressure amplitude (0.12–0.22 bar) and the initial gas concentration in the liquid (0.95–1.0 of the saturation concentration). The following new phenomena/features are observed.

• • Bubble growth rate increases with the applied pressure frequency, the applied pressure amplitude and the initial gas concentration in the liquid. The bubble growth rate is considerably more sensitive to changes in the applied pressure amplitude than the applied pressure frequency, but it is extremely sensitive to changes in the initial gas concentration in the liquid. Our recommendation to optimize bubble growth in applications is to increase the concentration of gas in the liquid.

• • Bubble growth rate decreases with increases in the initial volume fraction of gas in the bubble. Depending on how a bubble is generated, the initial volume fraction of gas may differ. Therefore, it may be difficult to obtain consistent results for bubble growth from different experiments. However, the threshold conditions appear to be insensitive to changes in the initial volume fraction of gas in the bubble.

• • In the absence of surfactants, the area and shell effects are the primary mechanisms in rectified diffusion, but some modifications are required to our understanding in order to account for nonlinear bubble oscillations. In the nonlinear oscillations, the maxima are further above the equilibrium bubble radius than the minima are below the equilibrium bubble radius. Therefore, the area and shell effects are promoted over and above what has been predicted in the literature (Young 1989).

• • In several numerical results for the bubble growth, there is a fall in the growth rate after some time, followed by a rapid increase, then a new, faster growth rate. The fall in the growth rate is associated with the local minimum amplitude of the gradient of the gas concentration at the bubble surface. A decrease in the shell effect is responsible for this fall in the growth rate. The shell effect may be reduced by the growth of a high frequency mode in the nonlinear oscillations of the bubble radius (or bubble resonance). At the resonance, the bubble oscillates with the driving acoustic frequency as well as the natural frequency, the latter being typically several times the driving acoustic frequency. As a result, the concentration of gas in the inner region of the liquid does not have enough time to adjust in order to support the shell effect. After the resonance, the rapid increase is due to the local maximum amplitude of the gradient of the gas concentration at the bubble surface. The subsequent new, faster growth rate corresponds to the increased area effects as the bubble has a larger radius.

Our numerical results have justified the experimental use of the maximum bubble radius to determine the underlying bubble growth for the air–water system. The only exception was found during the bubble resonance when the dynamics of the bubble radius produces oscillatory behaviour in the maximum bubble radius which is not reflected in the equilibrium bubble radius or the mass of gas in the bubble.

The analysis in this article may be extended to consider thermal effects. Thermal effects play an important role in non-invasive medical procedures and sonochemistry. Future research will incorporate a spherically symmetric equation for conservation of energy in the liquid and a spatially uniform equation for conservation of energy in the bubble (Young 1989). This thermal problem has three space scales and five time scales. The framework of the asymptotic analysis described in this article is essential for the understanding of the thermal problem.

The theoretical research in this article allows the investigation of the growth of an air bubble in water under quite general conditions. Smaller bubbles, which are of interest in sonoluminescence, may be studied without the advanced instrumentation required for experimental measurements (Lauterborn & Kurz 2010). A physical model for sonoluminescence must incorporate into the model in this article: thermal effects, compressibility of the liquid and plasma physics in the bubble. Such a sonoluminescence model may be developed building on the model in this article.