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 Reference Wang and Blake2010, Reference Wang and Blake2011). 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 Reference Brennen2002). 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 Reference Young1989). 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 (Reference Louisnard and Gomez2003) 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 Reference Wang1998, Reference Wang2014; Tiwari, Pentano & Freund Reference Tiwari, Pentano and Freund2015; Lechner et al. Reference Lechner, Koch, Lauterborn and Mettin2017; Zeng et al. Reference Zeng, Gonzalez-Avila, Dijkink, Koukouvinis, Gavaises and Ohl2018), limited by computational resources, accumulated numerical errors and/or robustness of numerical modelling. Wang et al. (Reference Wang, Liu, Corbett and Smith2022) 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 Reference Fyrillas and Szeri1994). 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. Reference Peñas-López, Parrales, Rodríguez-Rodríguez and van der Meer2016), 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 Reference Young1989). 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 Reference Young1989).

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 Reference Wang2016; Smith & Wang Reference Smith and Wang2018). 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 Reference Fyrillas and Szeri1995).

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 Reference Eller and Flynn1965); (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 (Reference Fyrillas and Szeri1994) 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 Reference Prosperetti1984)

(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 Reference Eller and Flynn1965). 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 Reference Eller and Flynn1965; Fyrillas & Szeri Reference Fyrillas and Szeri1994; Brennen Reference Brennen2013). 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 Reference Epstein and Plesset1950; Eller & Flynn Reference Eller and Flynn1965; Fyrillas & Szeri Reference Fyrillas and Szeri1994). (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 (Reference Fyrillas and Szeri1994). 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.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 (Reference Plesset and Zwick1952):

(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 (Reference Fyrillas and Szeri1994) obtained the same equation for their smooth problem on the same length and time scales. Fyrillas & Szeri (Reference Fyrillas and Szeri1994) 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 (Reference Eller1969). The theoretical predictions of threshold (Fyrillas & Szeri Reference Fyrillas and Szeri1994) have also been successful in applications related to stable single bubble sonoluminescence (Brenner, Hilgenfeldt & Lohse Reference Brenner, Hilgenfeldt and Lohse2002).

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 (Reference Smith and Wang2022). 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. Reference Erdélyi, Magnus, Oberhettinger and Tricomi1954), 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. Reference Peñas-López, Parrales, Rodríguez-Rodríguez and van der Meer2016), 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.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 Reference Smith1985) 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 Reference Burden and Faires1985). 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. Reference Lee, Kentish and Ashokkumar2005). 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. Reference Lee, Kentish and Ashokkumar2005) 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. Reference Lee, Kentish and Ashokkumar2005). 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. (Reference Lee, Kentish and Ashokkumar2005). We have $Pe \approx 1.6 \times 10^{5}$, which is very large. Lee et al. (Reference Lee, Kentish and Ashokkumar2005) 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. Reference Lee, Kentish and Ashokkumar2005). 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 (Reference Lamb1932) 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. (Reference Lee, Kentish and Ashokkumar2005) 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. (Reference Lee, Kentish and Ashokkumar2005) 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 (Reference Brennen2013)

(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 (Reference Fyrillas and Szeri1994), 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 Reference Crum1980). 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 (Reference Crum1980).

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 (Reference Crum1980). 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 Reference Crum1980). 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 (Reference Crum1980) 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.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 (Reference Louisnard and Gomez2003) 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 (Reference Louisnard and Gomez2003). 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 Reference Blake1949).

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 (Reference Church1988) 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 Reference Blake1949).