A theoretical model for the growth of spherical bubbles by rectified diffusion

Abstract

Rectified diffusion is a bubble growth phenomenon that occurs in acoustic fields. Despite the existence of a well-established spherically symmetric mathematical model, theoretical results have been unsuccessful in reproducing the bubble growth even in the case of a single spherical bubble in the bulk. In the latter case, the influence of surfactants and acoustic microstreaming have been speculated as the explanation for this disagreement. In this article, an exact solution for the leading-order concentration of gas in the liquid is determined. Using this exact solution, the well-established mathematical model is reduced to a system of two ordinary differential equations for the spherical bubble radius and the mass of gas in the bubble. This simplified model predicts the rapid bubble growth observed in experiments of a single spherical bubble in the bulk. The new results show that the bubble growth is asymptotically larger than at later times when the mass flux is limited by the slow diffusion of gas in a much larger region of the liquid surrounding the bubble. The new results also show that this bubble growth is relatively insensitive to a reduction in surface tension.

1. Introduction

The pressure in a non-oscillatory bubble must exceed the pressure in the liquid due to the effect of surface tension, thus there is a flow of gas from the bubble into the liquid and the bubble slowly dissolves. However, when a bubble undergoes large-amplitude oscillations in a sound field and the acoustic pressure amplitude exceeds a threshold value, the mass flow may be reversed (Blake 1949), this phenomenon being known as rectified diffusion. Rectified diffusion is associated with important applications in sonochemistry (Mason 2003), ultrasonic cleaning (Reuter et al. 2017) and biomedical ultrasound (Coussios & Roy 2008). The well-established spherically symmetric mathematical model for rectified diffusion in the bulk has been known for five decades (Eller & Flynn 1965). Unfortunately, this mathematical model has proved to be intractable over many cycles of bubble oscillation, the key challenge being the evaluation of the gas flux at the bubble interface. An additional theoretical difficulty is to develop a rigorous means by which to deduce the influence of convection associated with the bubble volume oscillations.

Evidence for the challenges in solving the well-established mathematical model may be found in the inconsistency with experimental results. Eller (1969) noted that his experimental rates of bubble growth exceed his theoretical values by a factor of more than 20. Acoustic microstreaming was suggested as being partly responsible for this disparity due to the increased availability of gas in the liquid for diffusion into the bubble. Crum (1980) further speculated that, at reduced surface tension, his experimental and theoretical growth rates diverge. Church (1988) found agreement with the experimental results of Eller (1969) at initial bubble radii below $35\ \mathrm {\mu }\rm {m}$, but his theoretical values were considerably less than experimental values for larger bubbles. Church (1988) attributed this discrepancy to acoustic microstreaming and surfactants on the bubble surface. The theoretical work of Fyrillas & Szeri (1994), which is still considered to be the state-of-the-art in this field (Lauterborn & Kurz 2010), quite seriously underpredicts the growth rate for larger bubbles in the experimental results of Eller (1969). Fyrillas & Szeri (1994) explained this slow growth rate in terms of the importance of surfactants even in relatively clean laboratory situations.

In the range close to a bubble's resonant radius, mass transfer due to rectified diffusion is assumed to increase. Furthermore, if acoustic microstreaming or non-spherical surface oscillations occurs during resonance, then the mass transfer increases by an order of magnitude (Gould 1974). Leong et al. (2010, 2011) showed an increase in microstreaming whenever surfactants are employed, the greatest increase in microstreaming being observed for the bulkiest surfactant headgroup. Furthermore, surface mode oscillations were shown to cause a dramatic increase in microstreaming. Soto et al. (2020) have experimentally proved that acoustic microstreaming is the leading-order contributor to the growth rate during resonance. In this case, the proximity of the wall triggers non-spherical bubble oscillations, thereby producing microstreaming during resonance that results in a convection-dominated growth. However, for air–water systems in figure 2 of Lee, Kentish & Ashokkumar (2005), it is possible to discern that the mass transfer slightly decreases during resonance at approximately $40\ \mathrm {\mu }\rm {m}$.

Peñas-López et al. (2016, 2017) and Soto et al. (2020) analysed the history effect in bubble growth and dissolution. By neglecting the convection of gas in the liquid, an exact solution for the gas flux at the bubble interface was obtained. Using this analytical solution, they studied the contribution of past mass transfer events between a gas bubble and its liquid surroundings towards the current diffusion-dominated growth or dissolution. In these articles, Peñas-López et al. (2016, 2017) and Soto et al. (2020) have shifted the focus of the literature towards the use of analytical solution methods rather than numerical approaches (Hilgenfeldt, Lohse & Brenner 1996). In this article, our approach is to obtain an analytical solution for the gas flux at the bubble interface without neglecting the influence of convection associated with the bubble volume oscillations.

The remainder of the paper is organized as follows. The well-established mathematical model for the rectified diffusion is described in § 2. In § 3, the framework of matched asymptotic expansions is employed and a Lagrangian change of variables is employed to transform the convection–diffusion equation into a diffusion equation. The leading-order problem for the concentration of gas in the inner liquid region is solved analytically using the Fourier sine transform and Fubini's theorem. Using this exact solution, the well-established mathematical model is reduced to a system of two ordinary differential equations for the bubble radius and the mass of gas in the bubble. In § 4, the numerical solutions are firstly compared with experimental observations. A parametric analysis is then undertaken with the above theory for the surface tension and the initial uniform concentration of the gas in the liquid. Finally, in § 5, this study is summarized and the key outcomes are identified.

2. Mathematical model

We now summarize the well-established mathematical model. This model assumes that bubbles are spherical, since their sizes 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. Furthermore, the flow of liquid surrounding the bubble is assumed to be spherically symmetric. The gas is assumed to be adiabatic and we adopt the polytropic relation with the polytropic index taken to be the ratio of specific heats. The viscous heating of the liquid is assumed to take place on a longer time scale than those considered herein. Phase change at the bubble interface is also neglected. The Rayleigh–Plesset equation for a spherical gas bubble in an incompressible liquid under adiabatic conditions is given by

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

where the pressure of the liquid at the bubble surface, $\bar {p}_l$, is given as follows:

(2.2)\begin{equation} \bar{p}_l = \bar{p}_{g0} \left ( \frac{\bar{m}}{\bar{m}_i} \right )^{\kappa} \left( \frac{\bar{R}_{eqi}}{\bar{R}} \right)^{3 \kappa} - \frac{2 \sigma}{\bar{R}} -(\bar{p}_{\infty} - \bar{p}_v) - \frac{4 \mu}{\bar{R}} \frac{\mbox{d} {\bar{R}}}{\mbox{d} {\bar{t}}}-\bar{p}_a \sin (2 {\rm \pi}\bar{\nu}_a \bar{t}), \end{equation}

in which $\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 pressure of the bubble gas at the mass of gas in the bubble $\bar {m}_i$ and bubble radius $\bar {R}_{eqi}$; $\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 applied pressure.

The concentration of gas in the liquid, $\bar {c}(\bar {r}, \bar {t})$, satisfies a spherically symmetric convection–diffusion equation of the form

(2.3)\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 \mbox{for} \ \bar{r} > \bar{R}(\bar{t}), \end{equation}

in which $\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.4)\begin{equation} \bar{c} \rightarrow \bar{C}_i \quad \mbox{as} \ \bar{r} \rightarrow \infty \end{equation}

and the boundary condition at the bubble surface is

(2.5)\begin{equation} \bar{c} (\bar{R}(\bar{t}), \bar{t}) = \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}

in which $\bar {C}_i$ is the initial uniform concentration of the gas in the liquid and $k_H$ is Henry's law constant (Eller & Flynn 1965; Fyrillas & Szeri 1994). The mass of gas in the bubble is governed by the conservation law

(2.6)\begin{equation} \frac{\mbox{d} {\bar{m}}}{\mbox{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.7a–d)\begin{equation} \bar{R}(0) = \bar{R}_{eqi}, \quad \frac{\mbox{d} {\bar{R}}}{\mbox{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}$. The only significant difference between this mathematical model and the model solved by Fyrillas & Szeri (1994) is that the partial pressure of gas in the bubble above depends on the mass of gas in the bubble.

Rectified diffusion requires the oscillation amplitudes to be large (Blake 1949). Experimental results in figure 6 of Leong et al. (2010) show the radial oscillation width to be between 10 % and 20 % of the bubble radius. Following Fyrillas & Szeri (1994), we scale the bubble radius by the initial equilibrium bubble radius because the amplitude of the radial oscillation and bubble radius itself are of the same order in the nonlinear model. We scale (2.1)–(2.7a–d) 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.8a,b)\begin{equation} \bar{t} = \frac{\bar{R}_{eqi}}{U} \tau, \quad \bar{r} = \bar{R}_{eqi} r, \end{equation}

in which $\varDelta = \bar {p}_{\infty } - \bar {p}_v$ and $U = \sqrt {{\varDelta }/{\rho }}$, where $\varDelta$ 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 becomes

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

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

(2.10)\begin{equation} \frac{\partial {\hat{c}}}{\partial {\tau}} + \frac{\mbox{d} {\hat{R}}}{\mbox{d} {\tau}} \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 \mbox{for} \ r > \hat{R}(\tau), \end{equation}

in which $Pe = U \bar {R}_{eqi} / D$ is the large Péclet number for mass transfer. The far-field boundary condition is

(2.11)\begin{equation} \hat{c} \rightarrow 0 \quad \mbox{as} \ r \rightarrow \infty \end{equation}

and the boundary condition at the bubble surface is

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

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

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

in which

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

The initial conditions are given by

(2.15a–d)\begin{equation} \hat{R}(0) = 1, \quad \frac{\mbox{d} {\hat{R}}}{\mbox{d} {\tau}}(0) = 0, \quad \hat{c} ( r, 0 ) = 0, \quad \hat{m}(0) = 1, \end{equation}

for $r > 1$. Equations (2.9)–(2.15a–d) are the dimensionless mathematical model to be studied.

3. Analysis

The typical value of the Péclet number for mass transfer, $Pe$, is greater than 100 000; therefore, this problem is susceptible to perturbation theory. In the following, we use the framework of matched asymptotic expansions. Our approach is to find an exact solution to the convection–diffusion equation for the leading-order concentration of gas in the inner liquid region. Using this exact solution, the flux of gas at the bubble surface may be evaluated. The mathematical model is then simplified to a system of two nonlinear ordinary differential equations. Two ordinary differential equations may be accurately numerically integrated even over many bubble oscillation cycles.

3.1. Lagrangian change of variables

We introduce the following change of independent variables which was first introduced by Plesset & Zwick (1952):

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

in which the new dependent variables are $R(t) = \hat {R}(\tau )$, $c(x,t) = \hat {c}(r,\tau )$ and $m(t) = \hat {m}(\tau )$. The dependent variables without bars or hats correspond to dimensionless variables in the Lagrangian frame of reference. Equation (2.10) may be rewritten as

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

The Rayleigh–Plesset equation (2.9) becomes

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

and the conservation of the mass of gas in the bubble (2.13) becomes

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

Equations (3.2)–(3.4) are singularly perturbed in $1/Pe^{1/2}$ having both an inner expansion when $x \ll Pe^{1/2}$ ($r=\hat {R}(\tau )+O(Pe^{-1/2})$) and an outer expansion when $x \sim Pe^{1/2}$ ($r=\hat {R}(\tau )+O(1)$). There are also multiple time scales. However, we will consider matched asymptotic expansions when $t \ll Pe^{1/2}$.

3.2. Inner region

We introduce asymptotic expansions taking the first two orders for the mass of gas in the bubble

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

as $Pe \rightarrow \infty$. Therefore, the leading-order problem for the concentration of gas in the liquid on this length scale and time scale becomes

(3.6)\begin{equation} \frac{\partial {c_0}}{\partial {t}} = R_0^{4} \frac{\partial^{2} {c_0}}{\partial {x}^{2}} \quad \mbox{for} \ x > 0 \end{equation}

subject to the boundary condition at the bubble surface

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

the initial condition

(3.8)\begin{equation} c_0(x,0)=0 \quad \mbox{for} \ x>0 \end{equation}

and a far-field condition which will follow from matching with the outer region.

At leading order, (3.4) becomes

(3.9)\begin{equation} \frac{\mbox{d} {m_0}}{\mbox{d} {t}} = 0, \end{equation}

with $m_0(0)=1$. We deduce that the mass of gas in the bubble does not change on this time scale at leading order, so $m_0 = 1$. At next order, we have

(3.10)\begin{equation} \frac{\mbox{d} {m_1}}{\mbox{d} {t}} = 4 {\rm \pi}M R_0^{4} \frac{\partial {c_0}}{\partial {x}} (0,t), \end{equation}

with the initial condition $m_1(0)=0$.

3.3. Outer region

We introduce the following change of variables ($y=O(1)$):

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

The concentration of gas in the outer liquid region is denoted $C$. The equation for mass transport (2.10) for $C(y, t) = \hat {c}(r, \tau )$ becomes

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

with far-field boundary condition

(3.13)\begin{equation} C \rightarrow 0 \quad \mbox{as} \ y \rightarrow \infty \end{equation}

and the initial condition

(3.14)\begin{equation} C = 0 \quad \mbox{on} \ t=0. \end{equation}

We introduce an expansion of the following form:

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

as $Pe \rightarrow \infty$. At leading order in (3.12), we obtain

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

with the initial condition $C_0(0)=0$ which is readily integrated to yield $C_0 = 0$.

3.4. Matching

We use limit matching technique in the form

(3.17)\begin{equation} \lim_{x \rightarrow \infty} c_0 = \lim_{y \rightarrow 0} C_0 \end{equation}

to obtain the matching condition on the inner region

(3.18)\begin{equation} c_0 \rightarrow 0 \quad\mbox{as} \ x \rightarrow \infty. \end{equation}

3.5. Leading-order inner solution

After the application of a Fourier sine transform to (3.6)–(3.8) and (3.18), the following initial value problem is obtained:

(3.19)\begin{equation} \frac{\partial {\tilde{c}}}{\partial {t}} + s^{2} R_0^{4} \tilde{c} = s R_0^{4} f \quad \mbox{with} \ \tilde{c}(s,0)=0, \end{equation}

where $f(t)$ was defined in (3.7) and

(3.20)\begin{equation} \tilde{c} (s,t) = \int_{x=0}^{\infty} c_0(x,t) \sin (s x) \,\mbox{d} x. \end{equation}

We solve this initial value problem to obtain

(3.21)\begin{equation} \tilde{c}(s,t) = \int_{\xi=0}^{t} f(\xi) R_0(\xi)^{4} s \exp \left({-}s^{2} \int_{\xi}^{t} R_0(\eta)^{4} \,\mbox{d} \eta \right) \mbox{d} \xi. \end{equation}

We then apply the inverse Fourier sine transform to find a solution for the concentration of gas in the liquid in the form

(3.22)\begin{equation} c_0(x,t) = \frac{2}{\rm \pi} \int_{s=0}^{\infty} \int_{\xi=0}^{t} f(\xi) R_0(\xi)^{4} \exp \left({-}s^{2} \int_{\xi}^{t} R_0(\eta)^{4} \,\mbox{d} \eta \right) \mbox{d} \xi s \sin( s x ) \,\mbox{d} s. \end{equation}

It is desirable to change the order of the integration in (3.22), but it is unclear whether the integral of the absolute value exists which is a requirement to apply Fubini's theorem. Therefore, integration by parts is applied to the inner integral which yields

(3.23)\begin{align} c_0(x,t) &= \frac{2}{\rm \pi} \int_{s=0}^{\infty} \left \{ f(t) - f(0) \exp \left({-}s^{2} \int_{0}^{t} R_0(\eta)^{4} \,\mbox{d} \eta \right) \right. \nonumber\\ &\quad \left. -\int_{\xi=0}^{t} f^{\prime}(\xi) \exp \left({-}s^{2} \int_{\xi}^{t} R_0(\eta)^{4} \,\mbox{d} \eta \right) \mbox{d} \xi \right \} \frac{\sin( s x )}{s}\, \mbox{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 determine the integrals in $s$ to yield

(3.24)\begin{align} c_0(x,t) &= f(t) - f(0) \,\mbox{erf} \left ( \frac{x}{2} \left \{ \int_{0}^{t} R_0(\eta)^{4} \,\mbox{d} \eta \right \}^{{-}1/2} \right ) \nonumber\\ &\quad -\int_{\xi=0}^{t} f^{\prime}(\xi)\, \mbox{erf} \left( \frac{x}{2} \left \{ \int_{\xi}^{t} R_0(\eta)^{4} \,\mbox{d} \eta \right \}^{{-}1/2} \right) \mbox{d} \xi, \end{align}

where

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

The solution (3.24) may be seen to be an exact solution of the leading-order inner approximation by substituting this formula back into (3.6)–(3.8) and (3.18).

Fyrillas & Szeri (1994) adopted the idea of linear superposition when they split their model into an oscillatory problem and a smooth problem. Their solution was a linear sum of the solutions to the oscillatory and smooth problems. Linear sums of solutions are inappropriate for nonlinear problems, they are an engineering approximation at best. Fyrillas & Szeri (1994) then solved (3.6) with Fourier series, but this problem has a continuous spectrum due to the semi-infinite spatial domain, and Fourier series are more appropriate for problems with a discrete spectrum on a finite spatial domain. They concluded that the mass transport associated with the transient of the oscillatory problem is negligible, but this is clearly incorrect as we describe in the mass flux below. Furthermore, we note that their solution to the smooth problem makes sense only near the threshold when the bubble gas in the liquid averaged over a cycle of oscillation must neither increase nor decrease.

The concentration gradient at the bubble interface located at $x=0$ is determined by taking the partial derivative of this analytical solution with respect to $x$, as follows:

(3.26)\begin{align} \frac{\partial {c_0}}{\partial {x}}(0,t) ={-}\frac{f(0)}{\sqrt{\rm \pi}} \left \{ \int_{0}^{t} R_0(\eta)^{4} \,\mbox{d} \eta \right \}^{{-}1/2} - \frac{1}{\sqrt{\rm \pi}} \int_{\xi=0}^{t} f^{\prime}(\xi) \left \{ \int_{\xi}^{t} R_0(\eta)^{4} \,\mbox{d} \eta \right \}^{{-}1/2} \,\mbox{d} \xi. \end{align}

The first term on the right-hand side of (3.26) corresponds to a slowly fading memory of the initial condition. The second term is the history term as described by Peñas-López et al. (2016). It represents a memory of the entire preceding time history of the first time derivative of the concentration at the bubble surface. The distinction between this expression for the mass flux and the corresponding one in Peñas-López et al. (2016) is that the influence of convection associated with the bubble volume oscillations is accounted for here. Each contribution of the first derivative of the concentration of gas at the bubble surface is weighted by the integral of the fourth power of the bubble radius. Therefore, the flux of dissolved gas has been determined as an exact solution of the leading-order inner problem. The key theoretical difficulty in this topic has been overcome (Fyrillas & Szeri 1994). The expression for the concentration gradient at the bubble interface (3.26) allows us to reduce the system of (2.9)–(2.15a–d) to a system of two ordinary differential equations for the growth and dissolution.

3.6. Simplified mathematical model

The simplified mathematical model is the dimensionless Rayleigh–Plesset equation

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

and conservation of the mass of gas in the bubble

(3.28)\begin{align} \frac{\mbox{d} {m}}{\mbox{d} {t}} ={-}\frac{4 {\rm \pi}^{1/2} M}{Pe^{1/2}} R^{4} \left [ f(0) \left \{ \int_{0}^{t} R(\eta)^{4} \,\mbox{d} \eta \right \}^{{-}1/2} + \int_{\xi=0}^{t} f^{\prime}(\xi) \left \{ \int_{\xi}^{t} R(\eta)^{4} \,\mbox{d} \eta \right \}^{{-}1/2} \,\mbox{d} \xi \right], \end{align}

with the initial conditions

(3.29a–c)\begin{equation} R(0) = 1, \quad \frac{\mbox{d} {R}}{\mbox{d} {t}}(0) = 0, \quad m(0) = 1, \end{equation}

in which

(3.30a,b)\begin{equation} f(0) ={-} 1 + {\rm d} p_{g0},\quad f^{\prime}(\xi) ={-} 3 \kappa \,{\rm d} p_{g0} \frac{m^{\kappa}}{R^{3 \kappa + 1}} \frac{\mbox{d} {R}}{\mbox{d} {\xi}}. \end{equation}

We remain consistent by neglecting terms of the order of $1/ Pe^{1/2}$ in the expression for $f^{\prime }(\xi )$. Equation (3.28) has reverted to the form (3.4) as it is equivalent to (3.9) and (3.10). The fourth power of $R$ on the right-hand side of (3.28) encapsulates both the area and shell effects responsible for the growth of the bubble. In the remainder of this article, we study this simplified initial value problem (3.27)–(3.29a–c) using numerical methods.

4. Numerical results

4.1. Algorithm

The system of ordinary differential equations (3.27)–(3.29a–c) is integrated numerically with the Numerical Algorithms Group (NAG) routine D02EJF. The time integrals on the right-hand side of (3.28) are evaluated with the trapezoidal composite rule and Bode's composite rule. Bode's rule is exploited to integrate over longer time intervals, which reduces the computer memory requirements and speeds up the evaluation of history effects whilst retaining the accuracy of the trapezoidal rule on shorter time intervals.

4.2. Validation

Experimental results do not exist in the published literature on the time scale of a few hundred bubble oscillations on which the simplified model (3.27)–(3.29a–c) is valid. In order to compare with the experiments of Eller (1969), we note the following. The main change in the physical balances occurs when time is of the order of the Péclet number for mass transfer, the outer expansion in the above matched asymptotic expansions does not become active until this time. The introduction of multiple scales when time is of the order of the square root of the Péclet number for mass transfer is to maintain the same physics. The mass of gas in the bubble is being forced to retain a short memory of the concentration gradient on the order one time scale and not being allowed to develop a long memory. There is no new physical balance on this time scale of the order of the square root of the Péclet number for mass transfer. Therefore, the underlying physics on the time scale less than the Péclet number for mass transfer will be unchanged. If we are able to determine the underlying physics during the times when our simplified model is valid, then this same physics will be valid for all time less than the order of the Péclet number for mass transfer.

We consider the case of air-saturated water in a vertical pipe (Eller 1969). A piston at the bottom of the pipe was used to excite a standing acoustic wave at $26.6$ kHz and a hypodermic needle was used to inject air bubbles into the pipe with a range of initial equilibrium radii. The following values for gas bubbles in water are adopted following the experimental conditions: $\varDelta = 10^{5}\ {\rm kg~m}^{-1}\ {\rm s}^{-2}$; $\sigma = 0.073\ {\rm Nm}^{-1}$; $\kappa = 1.4$, $\mu = 0.001\ {\rm Pas}$; $\rho = 10^{3}\ {\rm kg~m}^{-3}$; $\bar {\nu }_{a} = 26.6 \times 10^{3}$ Hz; $\bar {p}_a = 2 \times 10^{4}\ {\rm kg~m}^{-1}\ {\rm s}^{-2}$; $D = 2 \times 10^{-9}\ {\rm m}^{2}\ {\rm s}^{-1}$; and $k_{H} = 2.5 \times 10^{6}\ {\rm m}^{2}\ {\rm s}^{-2}$. The initial pressure of the bubble gas is given by $\bar {p}_{g0} = \varDelta + 2 \sigma / \bar {R}_{eqi}$ and the initial mass of gas in the bubble is given by

(4.1)\begin{equation} \bar{m}_i = \phi_{air} \rho_{air} \frac{4}{3} {\rm \pi}\bar{R}_{eqi}^{3}, \end{equation}

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 giving $\rho _{air}=1.3\ {\rm kg~m}^{-3}$. We adopt a reasonable estimate of the initial volume fraction of air in the bubble $\phi _{air}=0.9$. The saturation concentration of the gas in the liquid, $\bar {C}_{sat} = 4.15 \times 10^{-2}\ {\rm kg~m}^{-3}$, has been chosen to fit the threshold bubble radius of approximately $27\ \mathrm {\mu }\rm {m}$. We choose the initial uniform concentration of the gas in the liquid $\bar {C}_i = \bar {C}_{sat}$ (Eller 1969). The Péclet number for mass transfer takes values in the range $8 \times 10^{4}$ to $3 \times 10^{5}$. In order to validate the simplified model (3.27)–(3.29a–c), numerical solutions were obtained for these parameter values.

Based on the bubble buoyancy, Eller (1969) measured the equilibrium bubble radius rather than the bubble radius itself. The dimensional equilibrium bubble radius, $\bar {R}_{eq}$, is evaluated as the solution to the following algebraic equation:

(4.2)\begin{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. \end{equation}

We firstly consider two typical situations where the equilibrium bubble radius is larger and smaller than the threshold bubble radius, respectively. An example of the bubble growth is shown in figure 1(a) for an initial equilibrium bubble radius of $\bar {R}_{eqi} = 30\ \mathrm {\mu }\rm {m}$. The equilibrium bubble radius oscillates around a rapidly increasing mean value. The oscillations initially exhibit beating but eventually the frequency of the oscillations settles down to the applied pressure frequency. An example of the bubble dissolution is shown in figure 1(b) for an initial equilibrium bubble radius of $\bar {R}_{eqi} = 20\ \mathrm {\mu }\rm {m}$. In this case, the equilibrium bubble radius oscillates around a rapidly decreasing mean value and the frequency of the oscillations almost immediately follows the applied pressure frequency. The mean values of the equilibrium bubble radius, $\bar {a}$, are well approximated by a square root dependence in time. The square root dependence indicates the relative significance of the diffusion of air in the water for the bubble growth and dissolution as it corresponds to the scaling in the Boltzmann variable ($x/t^{1/2}$). The approximation by a square root holds except for values of the initial equilibrium bubble radius of $\bar {R}_{eqi}$ between $55\ \mathrm {\mu }{\rm m}$ and $70\ \mathrm {\mu }\rm {m}$ which will be discussed below.

Figure 1. The equilibrium bubble radius, $\bar {R}_{eq}$, using a numerical solution of the simplified mathematical model (3.27)–(3.29a–c) and an initial equilibrium bubble radius of (a) $\bar {R}_{eqi}=30\ \mathrm {\mu }\rm {m}$ and (b) $\bar {R}_{eqi}=20\ \mathrm {\mu }\rm {m}$. Power law fits for the mean value of the equilibrium bubble radius of (a) $\bar {a}=\bar {R}_{eqi}+7 \times 10^{-7} (\bar {t} U \bar {R}_{eqi})^{1/2}$ and (b) $\bar {a}=\bar {R}_{eqi}-2.4 \times 10^{-6} (\bar {t} U \bar {R}_{eqi})^{1/2}$ are compared with the numerical solution.

Figure 2 shows the oscillations of the bubble radius corresponding to figure 1(a). The periods of oscillation in figure 2 agree with the period of the applied external pressure of approximately $3.8 \times 10^{-5}$ s. We note that the maxima are farther above the equilibrium bubble radius than the minima are below. This feature supports the area and shell mechanisms of rectified diffusion (Young 1989).

Figure 2. The bubble radius, $\bar {R}$, and the equilibrium bubble radius, $\bar {R}_{eq}$, using a numerical solution of the simplified mathematical model (3.27)–(3.29a–c) and an initial equilibrium bubble radius of $\bar {R}_{eqi}=30\ \mathrm {\mu }\rm {m}$ corresponding to figure 1(a).

In order to compare with the time-averaged growth rates in the experiments of Eller (1969), we require to average the growth rate of the mean value of the equilibrium bubble radius. As the mean values of the equilibrium bubble radius follow a square root dependence in time, the time-averaged growth rate varies with the reciprocal of the square root of the time period. The time period over which the growth rate is averaged has been chosen to be $0.05$ s to fit with the experiments. The change in the mean value of the equilibrium bubble radius over this time period is of the order of a tenth of a micron, which is small in comparison with the equilibrium bubble radius itself. Numerical simulations are undertaken for a range of initial equilibrium radii $\bar {R}_{eqi}$. Figure 3 compares the time-averaged rate of growth of the mean value of the equilibrium bubble radius with the experimental results, the agreement being excellent up to an initial equilibrium bubble radius of $\bar {R}_{eqi} = 55\ \mathrm {\mu }\rm {m}$.

Figure 3. Comparison of the time-averaged rate of growth of the mean value of the equilibrium bubble radius, $\mbox {d} \bar {a}/ \mbox {d} \bar {t}$, using a numerical solution of the simplified mathematical model (3.27)–(3.29a–c), the experimental results of Eller (1969) and the theoretical results of Fyrillas & Szeri (1994). In the experiment, a bubble in air-saturated water is subject to a sound field of 26.6 kHz.

The bubble radius has a resonance at approximately $\bar {R}_{eqi}=60\ \mathrm {\mu }{\rm m}$ in our results, which has also been observed in the theoretical results of Fyrillas & Szeri (1994). Bubble resonance corresponds to the appearance of a high frequency mode in the bubble oscillations. Due to the inaccuracy of the polytropic bubble models near resonance as discussed in Appendix C of Fyrillas & Szeri (1994) and the absence of the square root dependence in time, results are omitted for the initial equilibrium bubble radius $\bar {R}_{eqi}$ between $55\ \mathrm {\mu }\rm {m}$ and $70\ \mathrm {\mu }\rm {m}$. In comparison with the predictions of Fyrillas & Szeri (1994), which are also shown in figure 3, our resonance affects a larger range of the initial equilibrium bubble radius. Between $55\ \mathrm {\mu }\rm {m}$ and $70\ \mathrm {\mu }{\rm m}$, it is possible that microstreaming governs the bubble growth rate in the experimental results of Eller (1969). For values of the initial equilibrium bubble radius greater than $70\ \mathrm {\mu }\rm {m}$, the square root dependence in time returns and the bubble growth rate is again dominated by the diffusion of air in the water. The bubble growth rate for values above the resonance are also shown in figure 3. The present theory displays a much better agreement with the experimental data.

Despite the comparison in figure 3, we are not saying that Fyrillas & Szeri (1994) is without its applications. For example, stable single bubble sonoluminescence requires that the bubble is diffusively stable (Brenner, Hilgenfeldt & Lohse 2002). The bubble gas averaged over a cycle of oscillation must neither increase nor decrease which corresponds to operating exactly on the threshold. Fyrillas & Szeri (1994) accurately predicts the threshold as shown in figure 3.

4.3. Parametric analysis

As discussed in the introduction, many authors have explained the underprediction of their theoretical growth rates in terms of the importance of surfactants. We now investigate how the time-averaged growth rates at $\bar {R}_{eqi}=40\ \mathrm {\mu }\rm {m}$ in figure 3 would vary in response to a reduction in surface tension. Figure 4(a) shows that the growth rate increases linearly in response to a reduction in surface tension. It is unclear how small quantities of surfactant could significantly increase the growth rate in terms of a small reduction in surface tension. Figure 5 of Lee et al. (2005) has also reported an approximately linear decrease of the bubble growth rate with surface tension for the addition of the surfactant sodium dodecyl sulphate.

Figure 4. The time-averaged rate of growth of the mean value of the equilibrium bubble radius, $\mbox {d} \bar {a}/ \mbox {d} \bar {t}$, using a numerical solution of the simplified mathematical model (3.27)–(3.29a–c) with $\bar {R}_{eqi}=40\ \mathrm {\mu }\rm {m}$ for (a) a range of values for the surface tension and (b) a range of values of the initial uniform concentration of the gas in the liquid.

Figure 4(b) shows that the bubble growth rate increases linearly with the initial uniform concentration of the gas in the liquid. The bubble growth rate is significantly more sensitive to changes in the initial concentration of the gas than it is to the surface tension.

5. Summary and conclusions

A theoretical study has been carried out to investigate the growth and dissolution of spherical bubbles subject to acoustic forcing. An exact leading-order solution is obtained for the gas flux at the bubble interface without neglecting the influence of convection associated with the bubble volume oscillations. The well-established mathematical model has been reduced to a system of two ordinary differential equations. Provided the Péclet number for mass transfer is sufficiently large, any errors in approximating the well-established mathematical model are negligibly small. Predictions of our simplified model have excellent agreement with experimental results. This new system of equations is important for the following reasons.

1. (1) These ordinary differential equations are straightforward to integrate numerically, whereas the well-established mathematical model has proved to be intractable.

2. (2) The spatial extent of the diffusion boundary layer in the liquid adjacent to the bubble interface is sufficiently small ($\bar {R}_{eqi}/Pe^{1/2}$), that the model applies to bubbles even in a dense cloud.

3. (3) This model, together with our recent research on strongly nonlinear analysis (Smith & Wang 2017, 2018, 2021), is the basis of our objective to simulate bubble growth over many millions of cycles of oscillation.

The rapid growth and dissolution has been shown to be dominated by the diffusion of gas in the liquid. The mass flux in this period is asymptotically larger than at later times when the flux is limited by the slow diffusion of gas in a much larger region of the liquid surrounding the bubble.

The time-averaged rate of growth of the mean value of the equilibrium bubble radius has been shown to have a slow linear decrease with the surface tension and a rapid linear increase with the initial uniform concentration of the gas in the liquid.