1 Introduction
The diffusion-driven growth and dissolution dynamics of bubbles of a soluble gas are topics that, despite having been studied for a long time, still awaken the interest of scientists and engineers. In addition to numerous numerical studies on different aspects of the matter, exhaustive analytical treatment in the canonical scenario of the mass (and heat) diffusion-driven growth or dissolution of an isolated bubble has been given over the years. The methods employed are based on the quasi-static approximation (Epstein & Plesset Reference Epstein and Plesset1950), thin boundary layer approximation (Plesset & Zwick Reference Plesset and Zwick1954), perturbation techniques (Duda & Vrentas Reference Duda and Vrentas1969), infinite series (Tao Reference Tao1978), integral methods (Rosner & Epstein Reference Rosner and Epstein1972) and self-similar solutions for bubble growth starting from zero initial size (Birkhoff, Margulies & Horning Reference Birkhoff, Margulies and Horning1958; Scriven Reference Scriven1959), to cite a few.
All these solutions have in common that they assume
-
(i) constant ambient pressure during the entire process and
-
(ii) a uniform concentration (or temperature) field in the liquid as an initial condition.
Perhaps the most widely used ones are the Epstein–Plesset solutions, valid for growth and dissolution, based on the quasi-stationary approximation (Epstein & Plesset Reference Epstein and Plesset1950), and Scriven’s exact solution for growth that accounts for the advection term in the diffusion equation (Scriven Reference Scriven1959). The predicted growth rates have been experimentally validated in supersaturated
$\text{CO}_{2}$
–water solutions by the works of Barker, Jefferson & Judd (Reference Barker, Jefferson and Judd2002) and Enríquez et al. (Reference Enríquez, Sun, Lohse, Prosperetti and van der Meer2014), among others. Likewise, in the case of dissolution, experimental verification of the Epstein–Plesset equation has been shown for monocomponent (Kapodistrias & Dahl Reference Kapodistrias and Dahl2012) and multicomponent (Shim et al.
Reference Shim, Wan, Hilgenfeldt, Panchal and Stone2014) bubbles.
With respect to the assumption (i) stated above, it happens that the effect of a non-constant pressure time history has been somewhat overlooked in the derivation of these analytical solutions. In some situations of practical interest, soluble bubbles are subject to successive slow compression–expansion cycles. Tisato et al. (Reference Tisato, Quintal, Chapman, Podladchikov and Burg2015) have successfully proven that the growth–dissolution dynamics of bubbles of a soluble gas can significantly damp the amplitude of seismic waves. A second example is the observation of gas bubble disease in stranded cetaceans after being exposed to low-frequency high-intensity acoustic pulses emitted by sonars. Crum & Mao (Reference Crum and Mao1996) attribute this to bubble growth triggered by rectified diffusion. Houser, Howard & Ridgway (Reference Houser, Howard and Ridgway2001) suggested that the likelihood of successful triggering is strongly dependent on the previous history of dives undergone by the cetacean. More specifically, the dive history directly determines the initial supersaturation level of dissolved nitrogen gas in the body fluid surrounding the trapped microbubbles at the time of insonation.
Other scenarios may entail fast, isolated changes in pressure rather than cyclic acoustic forcing. For instance, the growth of gas bubbles trapped in extravascular tissue is aided by steady decompression during depressurisation events (Marzella & Yin Reference Marzella and Yin1994). Payvar (Reference Payvar1987) studied the mass-transfer-driven bubble growth during rapid liquid decompressions typically found in hydraulic power recovery turbines. A numerical treatment based on the integral method developed by Rosner & Epstein (Reference Rosner and Epstein1972) was employed. This method requires the imposition of a parabolic concentration profile within a thermal boundary layer of time-varying thickness
$L(t)$
. Theofanous et al. (Reference Theofanous, Biasi, Isbin and Fauske1969) also performed a similar treatment for the growth of vapour bubbles in a superheated fluid under the effect of a decreasing pressure field.
However, as Jones & Zuber (Reference Jones and Zuber1978) pointed out, the integral method is only suitable for a monotonically decreasing pressure time history and this approach cannot be of general utility. Jones & Zuber (Reference Jones and Zuber1978) also studied the growth of vapour bubbles in a superheated liquid under variable pressure. They used the thin boundary layer approximation, which assumes that the temperature spatial variation,
${\rm\Delta}T$
, from the bubble surface temperature to that of the bulk fluid occurs within a thin thermal boundary layer of thickness
$L$
, much smaller than the bubble radius
$R$
. As a consequence, the moving bubble boundary could be modelled as a fixed Cartesian plane. The spherical geometry of the bubble was later included in the solution via a correction factor. Little attention was given to the history integral, consequence of the pressure history in time, that appears in their derivations. Instead they focused on bubble growth under a linearly decreasing pressure field.
While suitable for heat-induced growth, however, in mass-diffusion-driven growth (and specially dissolution), the thin boundary layer approximation
$L(t)/R(t)\ll 1$
may not be always valid. The concentration boundary layer thickness often becomes comparable to the bubble in size, regardless of the value of the diffusion coefficient. To prove this, we can take the well-known asymptotic solutions of Plesset & Zwick (Reference Plesset and Zwick1954), Birkhoff et al. (Reference Birkhoff, Margulies and Horning1958) or Scriven (Reference Scriven1959) for thermal diffusion growth – under assumptions (i) and (ii) – driven by a temperature difference
${\rm\Delta}T$
between the bubble boundary and the bulk fluid. The bubble radius scales as
$$\begin{eqnarray}R(t)\sim \mathit{Ja}_{\mathit{th}}\sqrt{D_{\mathit{th }}t},\quad \text{with }\mathit{Ja}_{\mathit{th}}=\frac{{\it\rho}_{l}c_{l}{\rm\Delta}T}{{\it\rho}_{g}h_{\mathit{fg}}}.\end{eqnarray}$$
$D_{\mathit{th}}$
refers to the thermal diffusion coefficient of the liquid,
$\mathit{Ja}_{\mathit{th}}$
is the Jakob number for heat transfer, where
$c_{l}$
is the specific heat of the liquid,
${\it\rho}_{l}$
is the liquid density,
${\it\rho}_{g}$
is the gas (vapour) density and
$h_{\mathit{fg}}$
is the latent heat. The thermal boundary layer evolves as
$L\sim \sqrt{D_{\mathit{th }}t}$
. It follows that for moderate and high superheats,
$L/R\sim \mathit{Ja}_{\mathit{th}}^{-1}\ll 1$
.
Equivalently, for mass-diffusion-controlled growth driven by a (molar) concentration difference
${\rm\Delta}C$
between the bubble boundary and the bulk fluid, Epstein & Plesset (Reference Epstein and Plesset1950) and Scriven (Reference Scriven1959) among others obtained
$$\begin{eqnarray}R(t)\sim \mathit{Ja}_{m}\sqrt{D_{m}t},\quad \text{with }\mathit{Ja}_{m}=\frac{M_{g}{\rm\Delta}C}{{\it\rho}_{g}}.\end{eqnarray}$$
$D_{m}$
denotes the mass diffusion coefficient,
$M_{g}$
is the gas molar mass and
$\mathit{Ja}_{m}$
(Szekely & Martins Reference Szekely and Martins1971) may be regarded as the analogous Jakob number for mass transfer. For small to moderate supersaturations, at long times the boundary layer thickness is of the order of the bubble radius,
$L\sim R$
, and analytical treatment accounting for a non-thin boundary layer is therefore essential. It should be pointed out that the scale laws just described correspond to cases of pure phase-change or mass-transfer-driven dynamics. The case where both phenomena are present is more involved and lies beyond the scope of the present work. The interested reader is referred to a recent paper by Fuster & Montel (Reference Fuster and Montel2015).
Turning now to assumption (ii), the effect of the previous growth–dissolution history of the bubble has not, to the best of our knowledge, been explored in detail. Naturally, if the bubble has been exchanging mass with its surroundings as a result of previous variations of the ambient pressure, then the concentration field at the beginning of the growth or dissolution stage of interest will not, in general, be uniform. Instead, it is determined by the boundary layer that has grown during that history. The temporal bubble dynamics has been found to be very sensitive to changes in the initial gas concentration profile dissolved in the liquid (Webb et al. Reference Webb, Arora, Roy, Payne and Coussios2010).
Figure 1. Experimental plot depicting the growth and dissolution dynamics of a sessile
$\text{CO}_{2}$
spherical bubble growing from a
$50~{\rm\mu}\text{m}$
pit on a flat chip immersed in a pressurized
$\text{CO}_{2}$
–water solution. The details of the experiments, performed in the facility described by Enríquez et al. (Reference Enríquez, Hummelink, Bruggert, Lohse, Prosperetti, van der Meer and Sun2013) and Enríquez et al. (Reference Enríquez, Sun, Lohse, Prosperetti and van der Meer2014) will be described in a companion paper. It shows (a) the evolution in time of the measured bubble radius
$R$
corresponding to (b) an imposed variation in the ambient pressure
$P_{\infty }(t)$
. The reference pressure,
$P_{c}\simeq 5.9~\text{bar}$
, is chosen to be the saturation pressure. Initially, the bubble has a stable radius
$R_{c}=225~{\rm\mu}\text{m}$
. At time
$T_{1}$
, the ambient pressure is set to drop from
$P_{c}$
to
$\simeq 5.45$
bar for a period of
${\rm\Delta}t=180~\text{ s}$
. This results in imminent bubble growth (white markers). The same scenario is then repeated at time
$T_{2}$
, resulting in a second growth cycle (dark markers). A 5-point moving average filter was performed on the time derivative of the measured radii. This allows a cleaner comparison of (c) the rate of growth of the pressure-corrected radius for the two growth cycles. The time axis is initialized on
$T_{1}$
or
$T_{2}$
accordingly. The uncertainty in the growth rate is estimated in
$\pm 0.0625~{\rm\mu}\text{m}~\text{s}^{-1}$
, much smaller than the differences between both cycles observed initially, at times
$t-T_{i}\lesssim 40~\text{ s}$
.
The effect of the previous history of the bubble on its current dynamics, hereon referred to as the history effect, can be observed from the experimental results in figure 1. Figure 1(a) depicts the evolution of the radius of a bubble under two expansion–compression cycles. Before the first growth cycle at
$t=T_{1}$
, the bubble is in equilibrium with its surroundings, implying a uniform concentration field at
$t\leqslant T_{1}$
. This is not the case for the second growth cycle. A non-uniform concentration profile is expected at
$t=T_{2}$
. Consequently, the measured growth rate is different. In this case it is larger, as will be explained in this paper, hence the bubble grows to a bigger size in the same growth time period. The growth rates are shown in figure 1(c). We have chosen to plot the pressure-corrected radius, defined as
$$\begin{eqnarray}R_{\mathit{corr}}(t)=R(t)\left(\frac{P_{\infty }(t)}{P_{c}}\right)^{1/3}.\end{eqnarray}$$
This way, the volumetric expansion of the bubble that is solely caused by the pressure drop (a fully non-diffusive effect due to Boyle’s law) is removed. The growth rates are initially different, but converge in time. It will be shown in this work that the observed differences in growth rate can be elegantly explained via a history integral of the bubble interfacial concentration, or alternatively, of the ambient pressure.
The objective of the present study is twofold:
-
(i) to present a theory for the diffusive-driven growth and dissolution of a spherical, isolated bubble that accounts for variable pressure time history. The dynamics evolves around an expression of similar form to (34) in Epstein & Plesset (Reference Epstein and Plesset1950) with an additional memory integral term that accounts for the history effect;
-
(ii) to illustrate the importance and nature of the history effect in a couple of bubble mass transfer processes of practical interest through analytical and numerical solutions.
Attending to these ideas, the paper is structured as follows. In § 2 the mathematical formulation is presented and the history integral term is derived. The governing equation for the bubble dynamics is then developed in § 3. Section 4 focuses on a bubble exposed to a train of piecewise constant-pressure steps. An analytical solution is presented, in addition to an illustrative example and a brief comparison of the aforementioned solution with numerical simulation. Section 5 then describes the role of the history effect on the potential growth of an isothermally oscillating bubble under harmonic pressure forcing. Finally, § 6 summarizes the main conclusions.
2 Formulation and derivation of the history integral term
2.1 Formulation
The growth rate of an isolated spherical gas bubble of radius
$R(t)$
suspended in a quiescent, infinite liquid environment subject to a time-varying liquid ambient pressure,
$P_{\infty }(t)$
, is to be determined. The analysis shall be restricted to monocomponent gas bubbles. The problem has spherical symmetry and consequently
$r$
may be taken as the radial distance from the bubble centre, while
$t$
is the time variable.
The molar concentration field
$C(r,t)$
of dissolved gas in the liquid is to be solved for from the advection–diffusion equation with spherical symmetry,
$$\begin{eqnarray}\frac{\partial C}{\partial t}+\frac{{\dot{R}}R^{2}}{r^{2}}\frac{\partial C}{\partial r}=D_{m}\frac{1}{r^{2}}\frac{\partial }{\partial r}\left(r^{2}\frac{\partial C}{\partial r}\right),\end{eqnarray}$$
where the dot notation stands for
$\text{d}/\text{d}t$
and
$D_{m}$
is the coefficient of mass diffusion. Notice that the factor
$u_{r}(r,t)={\dot{R}}R^{2}/r^{2}$
corresponds to the radial velocity field in the liquid that, by virtue of the continuity equation, is induced by the rate of change of the bubble radius. The concentration field is subject to boundary and initial conditions
$$\begin{eqnarray}C(R,t)=C_{s}(t),\quad C(\infty ,t)=C_{\infty },\quad C(r>R_{0},0)=C_{\infty },\end{eqnarray}$$
where
$R_{0}=R(0)$
is the initial bubble radius. The initial concentration of dissolved gas is assumed uniform throughout the liquid and equal to the concentration at the far field,
$C_{\infty }$
. The concentration boundary condition at the bubble surface is given by Henry’s law,
$$\begin{eqnarray}C_{s}=k_{H}P_{g},\end{eqnarray}$$
where
$k_{H}$
is Henry’s coefficient (molar based) and
$P_{g}(t)$
is the total gas pressure inside the bubble. The saturation pressure may thus be defined as
$P_{\mathit{sat}}=C_{\infty }/k_{H}$
. The bubble pressure is determined considering the liquid–gas surface tension
${\it\gamma}_{lg}$
, but otherwise neglecting the liquid vapour pressure together with damping and inertial effects. Naturally, we will not consider the case of an inertially pulsating bubble, driven by strong periodic acoustic forcing (Louisnard & Gomez Reference Louisnard and Gomez2003) or by resonance-triggering frequencies often employed in rectified diffusion. Indeed, in most diffusion-driven processes, excluding the aforementioned scenario, the contributions of the inertial and viscous terms in the Rayleigh–Plesset equation are negligible since bubble radius accelerations and velocities are relatively small. Payvar (Reference Payvar1987) reported that for bubbles of 0.1–1 mm in diameter,
$P_{\infty }(t)$
still remains the dominant term even for fast liquid decompressions of up to 100 bar taking place in less than 1 s. The bubble pressure is then simply obtained from the Young–Laplace equation,
$$\begin{eqnarray}P_{g}=P_{\infty }+\frac{2{\it\gamma}_{lg}}{R}.\end{eqnarray}$$
An isothermal liquid at temperature
$T_{\infty }$
, well below the boiling point, is considered. The bubble volume and pressure are related through the equation of state for an ideal gas,
$$\begin{eqnarray}{\textstyle \frac{4}{3}}{\rm\pi}R^{3}P_{g}=nR_{u}T_{\infty },\end{eqnarray}$$
where
$n(t)$
is the number of moles inside the bubble and
$R_{u}$
denotes the universal gas constant. Finally, Fick’s first law sets the molar flow rate of gas across the bubble surface to be
$$\begin{eqnarray}{\dot{n}}=4{\rm\pi}R^{2}D_{m}\left.\frac{\partial C}{\partial r}\right|_{r=R}.\end{eqnarray}$$
The governing equations (2.1)–(2.6) are best treated in dimensionless form. Consequently, let us define the dimensionless radius, ambient pressure, concentration and interfacial concentration as follows:
$$\begin{eqnarray}a=\frac{R}{R_{c}},\quad p=\frac{P_{\infty }}{P_{c}},\quad c=\frac{C-C_{\infty }}{k_{H}P_{c}},\quad c_{s}=\frac{C_{s}-C_{\infty }}{k_{H}P_{c}}.\end{eqnarray}$$
Here,
$R_{c}$
denotes some characteristic bubble radius. Pressure
$P_{c}$
is a characteristic liquid pressure, usually taken as the initial ambient pressure
$P_{\infty }(0)$
. Additionally, let us introduce the following dimensionless parameters, which remain constant throughout the process:
$$\begin{eqnarray}{\it\Upsilon}=\frac{C_{\infty }}{k_{H}P_{c}},\quad {\it\Lambda}=k_{H}R_{u}T_{\infty },\quad {\it\sigma}=\frac{2{\it\gamma}_{\mathit{lg}}}{R_{c}P_{c}}.\end{eqnarray}$$
The parameter
${\it\Upsilon}$
refers to the level of saturation of gas in the liquid at the characteristic pressure
$P_{c}$
. In fact, the (dimensionless) saturation pressure is simply
$p_{\mathit{sat}}=P_{\mathit{sat}}/P_{c}={\it\Upsilon}$
. Saturation conditions at
$P_{c}$
are described by
${\it\Upsilon}=1$
, while
${\it\Upsilon}<1$
and
${\it\Upsilon}>1$
imply undersaturation and supersaturation respectively. Similarly, note that the saturation level at
$P_{\infty }$
is given by
${\it\Upsilon}/p$
. Moreover,
${\it\sigma}$
represents the characteristic ratio of the Laplace pressure and
$P_{c}$
, while
${\it\Lambda}$
serves as the solubility parameter. The latter represents the ratio between the bubble’s volume and the volume of liquid needed to dissolve, under saturation conditions, the gas it contains. Henry’s law yields a new expression for
$c_{s}$
, which, when written in our notation, reads
$$\begin{eqnarray}c_{s}=(p+{\it\sigma}/a)-{\it\Upsilon}.\end{eqnarray}$$
Finally, we shall introduce two different dimensionless time variables,
${\it\tau}$
and
$\tilde{{\it\tau}}$
, in addition to the radial coordinate
${\it\xi}$
,
$$\begin{eqnarray}{\it\tau}=\frac{D_{m}}{R_{c}^{2}}t,\quad \text{d}\tilde{{\it\tau}}=\frac{D_{m}}{R^{2}(t)}\,\text{d}t,\quad {\it\xi}=\frac{r}{R(t)}.\end{eqnarray}$$
The identity
$$\begin{eqnarray}\text{d}{\it\tau}/\text{d}\tilde{{\it\tau}}=a^{2}\end{eqnarray}$$
directly relates the ‘physical’ dimensionless time
${\it\tau}$
and nonlinear time
$\tilde{{\it\tau}}$
. In the upcoming § 2.2, it will be shown that the mathematical nature of the problem encourages us to dispose of the ‘physical’ time
${\it\tau}$
and work with the nonlinear time
$\tilde{{\it\tau}}$
instead. We advance that this is done in order to analytically solve for the concentration gradient at the bubble boundary from the diffusion equation, all the while accounting for large variations in the bubble radius. Finally, coordinate
${\it\xi}\in [1,\infty )$
scales with the instantaneous bubble radius
$R(t)$
. Hence, the moving bubble boundary is advantageously always mapped by
${\it\xi}=1$
. Coordinate
${\it\xi}$
is an alternative to the Lagrangian coordinate
${\it\eta}=(r^{3}-R^{3}(t))/3$
which also eliminates the moving boundary problem. The latter transformation is usually employed in the treatment of rectified diffusion (Eller & Flynn Reference Eller and Flynn1965; Fyrillas & Szeri Reference Fyrillas and Szeri1994). However, in the diffusion-dominant regime of interest, the resulting advection–diffusion equation in
${\it\eta}$
is much harder, if not impossible, to treat analytically compared to the analogous equation in
${\it\xi}$
.
2.2 The history term on the concentration gradient at the bubble surface
The advection–diffusion equation (2.1) in dimensionless form becomes
$$\begin{eqnarray}\frac{\partial c}{\partial \tilde{{\it\tau}}}+\frac{1}{a}\frac{\text{d}a}{\text{d}\tilde{{\it\tau}}}\left(\frac{1}{{\it\xi}^{2}}-{\it\xi}\right)\frac{\partial c}{\partial {\it\xi}}=\frac{1}{{\it\xi}^{2}}\frac{\partial }{\partial {\it\xi}}\left({\it\xi}^{2}\frac{\partial c}{\partial {\it\xi}}\right).\end{eqnarray}$$
Note that the nonlinear time
$\tilde{{\it\tau}}$
is the time variable of choice for reasons that will soon become apparent. The boundary and initial conditions (2.2) of the concentration field
$c({\it\xi},\tilde{{\it\tau}})$
become
$$\begin{eqnarray}c(1,\tilde{{\it\tau}})=c_{s}(\tilde{{\it\tau}})=p+{\it\sigma}/a-{\it\Upsilon},\quad c(\infty ,\tilde{{\it\tau}})=0,\quad c({\it\xi}>1,0)=0.\end{eqnarray}$$
Only the dimensionless advection term in (2.12) is explicitly dependent on the bubble dynamics through the prefactor
$(\text{d}a/\text{d}\tilde{{\it\tau}})/a$
. We may characterize this prefactor as a time-dependent Péclet number based on the velocity of the bubble boundary,
$$\begin{eqnarray}\mathit{Pe}(\tilde{{\it\tau}})=\frac{1}{a}\frac{\text{d}a}{\text{d}\tilde{{\it\tau}}}\equiv R{\dot{R}}/D_{m}\end{eqnarray}$$
whose magnitude is essentially associated with the instantaneous ratio of advective to diffusive transport. To illustrate the role of the history effect on the diffusive growth or dissolution rates of bubbles, we will restrict the analysis to mass transfer processes where the advection term in (2.12) is small compared to the diffusion term. We anticipate that the history effect is present in every transient diffusion–advection problem in which boundary conditions change in time. However, should the advective mass transport of species be dominant, the advective nature of the flow would then, at least partially, obscure the contribution of the history effect on the bubble dynamics as well as of any other effect of a diffusive nature.
When can we neglect the advection term? To find out, we can perform an order of magnitude analysis on (2.12), taking
$O(\partial c)={\rm\Delta}c$
,
$O({\it\xi})=1$
and denoting
$O(\partial {\it\xi})$
as
$l$
, a characteristic thickness of the boundary layer, defined as
$\partial c/\partial {\it\xi}|_{{\it\xi}=1}={\rm\Delta}c/l$
. The analysis reveals that the diffusion term magnitude is
${\rm\Delta}c/l^{2}$
, while the advection term has magnitude
$O(\mathit{Pe}){\rm\Delta}c$
. It follows from the Epstein–Plesset solution that the penetration length of the mass boundary layer will never be much greater than the size of the bubble:
$l\lesssim 1$
. Thus, the dimensionless advection term may be neglected provided
$|\mathit{Pe}(\tilde{{\it\tau}})|\ll 1$
at all times.
The negligible advection assumption is usually valid for small supersaturation or undersaturation ratios, defined as
$C_{\infty }/[k_{H}P_{\infty }(t)]-1\equiv {\it\Upsilon}/p-1$
. Larger ratios are allowed if the gas solubility is very poor (
${\it\Lambda}\ll 1$
). In fact, it will be shown in a companion paper that, in the case of a growing or dissolving bubble at pressure
$P_{\infty }$
, the Péclet number is related to the Jakob number, defined in (1.2b
), through:
$|\mathit{Pe}|\approx Ja_{m}^{2}={\it\Lambda}|{\it\Upsilon}/p-1|$
. Thus, a highly soluble gas such as
$\text{CO}_{2}$
gas in water (
${\it\Lambda}\approx 0.8$
at room temperature) will require that
$|{\it\Upsilon}/p-1|\ll 1$
for advection to be negligible.
It must be emphasized that the advection term has two components. The component containing
$1/{\it\xi}^{2}$
is the (dimensionless) physical advection term, associated with the radial component of the velocity field,
$u_{r}(r,t)=R^{2}{\dot{R}}/r^{2}$
. The component containing
$-{\it\xi}$
is the artificially induced advection that compensates for the scaling nature of
${\it\xi}$
with
$R(t)$
. It naturally arises from the non-dimensionalization process. Neglecting both terms essentially means that we are now rigorously solving the diffusion equation with an ever-present numerical advection term (associated to the non-physical velocity field
$u_{r}(r,t)={\dot{R}}r/R$
) that accounts for the moving boundary (Peñas López, Parrales & Rodríguez-Rodríguez Reference Peñas López, Parrales and Rodríguez-Rodríguez2015). Nonetheless, the effect of the artificial advection on the interfacial concentration gradient will still be small for
$|\mathit{Pe}(\tilde{{\it\tau}})|\ll 1$
.
The concentration field may then be sought by solving
$$\begin{eqnarray}\frac{\partial c}{\partial \tilde{{\it\tau}}}=\frac{1}{{\it\xi}^{2}}\frac{\partial }{\partial {\it\xi}}\left({\it\xi}^{2}\frac{\partial c}{\partial {\it\xi}}\right),\end{eqnarray}$$
which has no direct dependency on the bubble dynamics. The choice of
$\tilde{{\it\tau}}$
over
${\it\tau}$
is justified in that
$\tilde{{\it\tau}}$
allows us to go from (2.1) to (2.15) based on the single assumption that the advection term is small. Had we chosen to perform an equivalent non-dimensionalization with physical time
${\it\tau}$
(cf. later (5.6)), it is only possible to treat the arising
$a({\it\tau})$
-dependent diffusive term employing an additional approximation, namely treating the radius
$a$
as a constant. The essence of the so-called quasi-stationary approximation (Weinberg & Subramanian Reference Weinberg and Subramanian1980) behind the Epstein–Plesset equation resides in these two approximations, namely (i) dropping the advective transport term arising from the interface motion and (ii) treating the bubble radius as a constant in the concentration boundary condition at the interface. Approximation (ii), hereon referred to as the quasi-static radius approximation, is only suitable when considering small or slow changes in the radius size from the equilibrium or initial size. Otherwise, this approximation will inherently decrease the accuracy of the solution. It is then concluded that the purpose of
$\tilde{{\it\tau}}$
is to avert making use of the quasi-static radius approximation, thereby extending the parameter range in which the theory is valid. Additionally, writing the advection–diffusion equation in terms of
$\tilde{{\it\tau}}$
allows us to point out an interesting conclusion regarding the appropriateness of neglecting advection effects. Indeed, the advective term in (2.12) is exactly zero at the bubble interface, as
$(1/{\it\xi}^{2}-{\it\xi})=0$
there. Thus, in a region close to the bubble surface, it is reasonable to expect the advective term to play a small role even for moderate values of the Péclet number.
The solution to (2.15) for the concentration gradient across the bubble interface (see appendix A) reads
$$\begin{eqnarray}-\!\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}=c_{s}+\frac{c_{s0}}{\sqrt{{\rm\pi}\tilde{{\it\tau}}}}+\int _{0}^{\tilde{{\it\tau}}}\frac{1}{\sqrt{{\rm\pi}(\tilde{{\it\tau}}-\tilde{x})}}\frac{\text{d}c_{s}}{\text{d}\tilde{x}}\,\text{d}\tilde{x}.\end{eqnarray}$$
The initial interfacial concentration,
$c_{s0}$
, is given by
$$\begin{eqnarray}c_{s0}=c_{s}(0)=p_{0}+{\it\sigma}/a_{0}-{\it\Upsilon},\end{eqnarray}$$
where
$p_{0}=p(0)$
and
$a_{0}=a(0)$
denote the initial pressure and radius respectively. The first term in the right-hand side of (2.16) corresponds to the steady-state solution determined by the instantaneous interfacial concentration,
$c_{s}(\tilde{{\it\tau}})$
. The second term is the transient component exclusively associated with the initial conditions. The third and final term is the history or memory integral term. It depends on the time history of
$c_{s}(\tilde{{\it\tau}})$
caused by any prior variations in the ambient pressure
$p(\tilde{{\it\tau}})$
. It addresses the temporal delay required for the concentration boundary layer to become fully developed (physically due to the finite diffusivity
$D_{m}$
) as the boundary condition, in this case
$c_{s}$
, changes with time. In other words, the current concentration profile and therefore the mass flux across the bubble are conditioned by the preceding time history of the boundary condition.
The history integral may be evaluated analytically for the cases when
$c_{s}$
varies in sudden (step-like) jumps and when
$c_{s}$
varies harmonically in time. These will be addressed later in § 4 and § 5 respectively.
3 The Epstein–Plesset equation with history term
The ideal gas equation of state provided in (2.5) may be combined with Fick’s law, equation (2.6), to obtain the following mass conservation equation:
$$\begin{eqnarray}4{\rm\pi}R^{2}{\dot{R}}\left(P_{\infty }+\frac{2{\it\gamma}_{sl}}{R}\right)+\frac{4{\rm\pi}R^{3}}{3}\left({\dot{P}}_{\infty }-\frac{2{\it\gamma}_{sl}}{R^{2}}{\dot{R}}\right)=4{\rm\pi}R^{2}\bar{R}T_{\infty }D\left.\frac{\partial C}{\partial r}\right|_{r=R}.\end{eqnarray}$$
In dimensionless form, it reads
$$\begin{eqnarray}\frac{\text{d}a}{\text{d}\tilde{{\it\tau}}}\left(p+\frac{{\it\sigma}}{a}\right)+\frac{1}{3}a\left(\frac{\text{d}p}{\text{d}\tilde{{\it\tau}}}-\frac{{\it\sigma}}{a^{2}}\frac{\text{d}a}{\text{d}\tilde{{\it\tau}}}\right)={\it\Lambda}a\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}.\end{eqnarray}$$
Inserting the expression for the concentration gradient (2.16) into (3.2) yields the governing equation for the bubble radius dynamics. This may be regarded as the Epstein–Plesset with history term (EPH) equation, written below in terms of
$c_{s}$
:
$$\begin{eqnarray}\frac{\text{d}a}{\text{d}\tilde{{\it\tau}}}(c_{s}+{\it\Upsilon})+\frac{1}{3}a\frac{\text{d}c_{s}}{\text{d}\tilde{{\it\tau}}}=-{\it\Lambda}a\left[c_{s}+\frac{c_{s0}}{\sqrt{{\rm\pi}\tilde{{\it\tau}}}}+\int _{0}^{\tilde{{\it\tau}}}\frac{1}{\sqrt{{\rm\pi}(\tilde{{\it\tau}}-\tilde{x})}}\frac{\text{d}c_{s}}{\text{d}\tilde{x}}\,\text{d}\tilde{x}\right].\end{eqnarray}$$
The liquid pressure has been related to the surface concentration through
$p=c_{s}+{\it\Upsilon}-{\it\sigma}/a$
using the relation previously provided in (2.9). In the absence of surface tension,
${\it\sigma}=0$
, the EPH equation (3.3) may then be conveniently expressed directly in terms of the ambient pressure
$p$
,
$$\begin{eqnarray}\frac{1}{a}\frac{\text{d}a}{\text{d}\tilde{{\it\tau}}}+\frac{1}{3p}\frac{\text{d}p}{\text{d}\tilde{{\it\tau}}}=-\frac{{\it\Lambda}}{p}\left[p-{\it\Upsilon}+\frac{p_{0}-{\it\Upsilon}}{\sqrt{{\rm\pi}\tilde{{\it\tau}}}}+\int _{0}^{\tilde{{\it\tau}}}\frac{1}{\sqrt{{\rm\pi}(\tilde{{\it\tau}}-\tilde{x})}}\frac{\text{d}p}{\text{d}\tilde{x}}\,\text{d}\tilde{x}\right].\end{eqnarray}$$
To obtain the radius evolution in the physical time,
$a({\it\tau})$
, we must numerically integrate the differential EPH equation in (3.3) or (3.4) for
$a(\tilde{{\it\tau}})$
in addition to the differential equation (2.11) for
${\it\tau}(\tilde{{\it\tau}})$
.
3.1 Recovering the Epstein–Plesset equation from the EPH
We now show that the Epstein–Plesset equation may be recovered from the EPH equation. In order to do so we revert back to (3.2). Since the original Epstein–Plesset equation is formulated using the linear time
${\it\tau}=D_{m}t/R_{c}^{2}$
, it is convenient for both time derivatives in
$\tilde{{\it\tau}}$
to be replaced by derivatives in
${\it\tau}$
instead. Recalling that
$\text{d}{\it\tau}=a^{2}\,\text{d}\tilde{{\it\tau}}$
, we may recast (3.2) as
$$\begin{eqnarray}\frac{\text{d}a}{\text{d}{\it\tau}}\left(p+\frac{2{\it\sigma}}{3a}\right)+\frac{1}{3}a\frac{\text{d}p}{\text{d}{\it\tau}}=\frac{{\it\Lambda}}{a}\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}.\end{eqnarray}$$
The analytical expression for the interfacial concentration gradient is still a function of
$\tilde{{\it\tau}}$
as given by (2.16). The Epstein–Plesset equation assumes a constant pressure history, i.e.
$p=p_{0}$
and
$\text{d}p/\text{d}{\it\tau}=0$
. Equation (3.5) then reduces to
$$\begin{eqnarray}\frac{\text{d}a}{\text{d}{\it\tau}}\left(p_{0}+\frac{2{\it\sigma}}{3a}\right)=\frac{{\it\Lambda}}{a}\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}.\end{eqnarray}$$
The interfacial concentration gradient deserves special treatment. The history integral term vanishes, consequence of the constant-pressure condition. Moreover, unlike the EPH equations (3.3) or (3.4), the Epstein–Plesset equation is based on the quasi-static radius approximation (Weinberg & Subramanian Reference Weinberg and Subramanian1980). The concentration gradient is calculated for a fixed bubble boundary of radius
$a$
. The solution is then coupled with the equation for mass conservation (3.5), where
$a$
is now time dependent. This has two implications: (i) the bubble radius
$a$
is treated as constant when solving for the concentration gradient. Consequently, (2.11) simplifies to
${\it\tau}=a^{2}\tilde{{\it\tau}}$
. Additionally, (ii) the initial interfacial concentration calculated for a static
$a$
is in fact ‘reused’ for all its values. This amounts to setting
$a_{0}=a$
in (2.17), hence
$c_{s0}=c_{s}=p_{0}+{\it\sigma}/a-{\it\Upsilon}$
. Under the quasi-static radius approximation, the concentration gradient in (2.16) then becomes
$$\begin{eqnarray}-\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}=\left[c_{s}+\frac{c_{s0}}{\sqrt{{\rm\pi}\tilde{{\it\tau}}}}\right]\approx \left(p_{0}-{\it\Upsilon}+\frac{{\it\sigma}}{a}\right)\left[1+\frac{a}{\sqrt{{\rm\pi}{\it\tau}}}\right].\end{eqnarray}$$
Inserting this expression into (3.6), one finally recovers the Epstein–Plesset equation in dimensionless form (Epstein & Plesset Reference Epstein and Plesset1950, equation (34)):
$$\begin{eqnarray}\frac{\text{d}a}{\text{d}{\it\tau}}=-\frac{{\it\Lambda}(p_{0}-{\it\Upsilon}+{\it\sigma}/a)}{p_{0}+2{\it\sigma}/(3a)}\left[\frac{1}{a}+\frac{1}{\sqrt{{\rm\pi}{\it\tau}}}\right].\end{eqnarray}$$
An equivalent equation in
$\tilde{{\it\tau}}$
, which does not make use of the quasi-static radius approximation, may be shown to be
$$\begin{eqnarray}\frac{\text{d}a}{\text{d}\tilde{{\it\tau}}}=-\frac{{\it\Lambda}a}{p_{0}+2{\it\sigma}/(3a)}\left[(p_{0}-{\it\Upsilon})\left(1+\frac{1}{\sqrt{{\rm\pi}\tilde{{\it\tau}}}}\right)+\frac{{\it\sigma}}{a}\left(1+\frac{a/a_{0}}{\sqrt{{\rm\pi}\tilde{{\it\tau}}}}\right)\right].\end{eqnarray}$$
In the absence of surface tension,
${\it\sigma}=0$
, (3.9) becomes a separable differential equation in
$a$
and
$\tilde{{\it\tau}}$
. Its analytical solution, with
$a(0)=a_{0}$
, is easily found to be
$$\begin{eqnarray}a(\tilde{{\it\tau}})=a_{0}\exp \left\{-{\it\Lambda}(p_{0}-{\it\Upsilon})\left(\tilde{{\it\tau}}+2\sqrt{\tilde{{\it\tau}}/{\rm\pi}}\right)\right\}.\end{eqnarray}$$
The exact analytical solution of
$a({\it\tau})$
to (3.10) exists in parametric form, with
$a(\tilde{{\it\tau}})$
and
${\it\tau}(\tilde{{\it\tau}})$
. The physical time
${\it\tau}$
is related to
$\tilde{{\it\tau}}$
through
$$\begin{eqnarray}\displaystyle {\it\tau} & = & \displaystyle \frac{a_{0}}{2{\it\Lambda}({\it\Upsilon}-p_{0})}\left[1-\text{e}^{-2{\it\Lambda}({\it\Upsilon}-p_{0})(\tilde{{\it\tau}}+2\sqrt{\tilde{{\it\tau}}/{\rm\pi}})}+\sqrt{2{\it\Lambda}({\it\Upsilon}-p_{0})}\text{e}^{2{\it\Lambda}({\it\Upsilon}-p_{0})}\right.\nonumber\\ \displaystyle & & \displaystyle \times \!\left.\left\{\text{erf}\!\left(\!\sqrt{2{\it\Lambda}({\it\Upsilon}-p_{0})/{\rm\pi}}\right)-\text{erf}\!\left(\sqrt{2{\it\Lambda}({\it\Upsilon}-p_{0})\tilde{{\it\tau}}}+\sqrt{2{\it\Lambda}({\it\Upsilon}-p_{0})/{\rm\pi}}\right)\!\right\}\right]\!,\qquad ~\end{eqnarray}$$
which can be derived by solving (2.11).
4 Multiple step-like variations of the ambient pressure
This section intends to shed light on the history effect in bubble growth through illustrative examples based on analytical and numerical solutions. To this end, we shall consider the case where the time-dependent ambient pressure consists of
$N$
consecutive step-like jumps in pressure. At the
$n\text{th}$
jump taking place at nonlinear time
$\tilde{T}_{n}$
, the pressure changes by
${\rm\Delta}p_{n}$
. The pressure history and its time derivative may be modelled as
$$\begin{eqnarray}p=p_{0}+\mathop{\sum }_{n=1}^{N}{\rm\Delta}p_{n}H(\tilde{{\it\tau}}-\tilde{T}_{n}),\quad \frac{\text{d}p}{\text{d}\tilde{{\it\tau}}}=\mathop{\sum }_{n=1}^{N}{\rm\Delta}p_{n}{\it\delta}(\tilde{{\it\tau}}-\tilde{T}_{n}),\quad \text{for }n=1,\ldots ,N,\end{eqnarray}$$
where
$H$
denotes the Heaviside function and
${\it\delta}$
is the Dirac delta. As anticipated at the end of § 2, it is then possible to analytically evaluate the history integral term and consequently solve the EPH equation provided the Laplace pressure is neglected. The analytical solutions to the EPH equations (3.4) and (2.15) for the concentration field are derived in appendix B.
Table 1. Coefficients
$\tilde{T}_{n}$
and
${\rm\Delta}p_{n}$
of the pressure step function used in the example (see figure 2). The resulting coefficients
$T_{n}$
associated with the physical time
${\it\tau}$
are also included.
4.1 An example to illustrate the history effect
Let us consider the radius history of a bubble with negligible Laplace pressure (
${\it\sigma}=0$
) initially in equilibrium at the saturation pressure. Setting
$P_{c}=P_{\infty }(0)$
and
$R_{c}=R(0)$
renders
$p_{0}=1$
,
$a_{0}=1$
and
${\it\Upsilon}=1$
. The bubble is then exposed to
$N=6$
consecutive jumps in pressure. We prescribe coefficients
$\tilde{T}_{n}$
and
${\rm\Delta}p_{n}$
, tabulated in table 1. These are entered into (4.1), resulting in the pressure time history plotted in figure 2(b). An iterative procedure was employed to establish coefficients
$\tilde{T}_{n}$
such that the two expansion (growth) stages (
$T_{1}^{+}<{\it\tau}<T_{2}^{-}$
and
$T_{4}^{+}<{\it\tau}<T_{5}^{-}$
) have the same duration in
${\it\tau}$
(see figure 2
d). Moreover, both expansion stages are identical in the sense that
$a(T_{1})=a(T_{4})=a_{0}=1$
,
${\rm\Delta}p_{1}={\rm\Delta}p_{4}$
and
$p_{1}=p_{4}$
. The analytical solution for
$a(\tilde{{\it\tau}})$
provided by (B 3) and (B 4) is plotted in figure 2(b). The solution is more naturally interpreted when presented in
${\it\tau}$
(figure 2
c). This requires numerical integration of (2.11) once the analytical solution for
$a(\tilde{{\it\tau}})$
has been obtained.
At
${\it\tau}=T_{1}^{-}$
, right before the first growth stage, the history integral is identically zero, i.e. there is no previous history. This is obviously not the case at
${\it\tau}=T_{4}^{-}$
. Consequently, the memory effect must be behind the differences in bubble growth observed between these two pressure-wise identical expansion periods (figure 2
a,c).
Figure 2. Analytical solution for (a) the evolution of the dimensionless bubble radius,
$a$
, corresponding to (b) multiple jumps in the dimensionless ambient pressure
$p$
, both plotted against the dimensionless, nonlinear time
$\tilde{{\it\tau}}$
. Initial conditions correspond to perfect saturation. Additionally, equivalent plots of (c) the dimensionless bubble radius and (d) dimensionless ambient pressure are plotted against the dimensionless linear time
${\it\tau}$
. For reference, the physical parameters employed are those considering a
$\text{CO}_{2}$
bubble of size
$R_{c}=225~{\rm\mu}\text{m}$
in water (
$k_{H}=3.40\times 10^{-4}~\text{mol}~\text{N}^{-1}~\text{m}^{-1}$
,
$D_{m}=1.92\times 10^{-9}~\text{m}^{2}~\text{s}^{-1}$
) under conditions
$P_{c}=4.9$
bar and
$T_{\infty }=293~\text{ K}$
. Surface tension is neglected. The resulting solubility and saturation parameters are
${\it\Lambda}=0.828$
and
${\it\Upsilon}=1$
.
To highlight this effect, figure 3(a) compares the bubble radius history during both growth stages. At equal times, measured from the beginning of each expansion period, the bubble always exhibits a larger radius in the second growth cycle. This is consequence of the larger, history-augmented growth rate (figure 3 b). The second growth rate is most enhanced at the beginning. It then asymptotically converges in time with the first growth rate, by then confirming the complete dissipation of the history effect.
The physical explanation of the history effect lies in the concentration profiles near the bubble a short time after each jump (figure 3
b). Both concentration profiles are bounded by an identical value of
$c_{s}({\it\tau}^{\ast })<0$
, and the far-field concentration
$c(\infty ,{\it\tau})=0$
. The first concentration profile corresponds to a uniform initial condition
$c({\it\xi},T_{1}^{-})=0$
. On the other hand, the second concentration profile remembers that
$c_{s}>0$
during the previous dissolution period at
$T_{2}^{+}<{\it\tau}<T_{3}^{-}$
. Consequently, there exists a supersaturation region (
$c>0$
) near the bubble, containing the mass of gas that was transferred to the liquid during said dissolution period. The concentration boundary layer takes time to fully develop in response to any change in the interfacial concentration. This effective adaptation time, owing to the finite diffusivity, results in this case in a momentarily steeper interfacial concentration gradient that enhances the growth rate substantially.
Figure 3. The first (dashed line) and second (solid line) growth cycles of the bubble (see figure 2) are compared through (a) the dimensionless radius
$a$
and (b) its rate of change, computed numerically out of the analytical
$a({\it\tau})$
. These are plotted against the dimensionless linear time
${\it\tau}$
after each of the two negative pressure jumps (
${\rm\Delta}p_{1}$
and
${\rm\Delta}p_{4}$
) that lead to undersaturation. Inset: dimensionless concentration radial profiles evaluated at a short time
${\it\tau}^{\ast }$
immediately after each of the two jumps.
The effect described here is also present in the experimental data shown in figure 1. As a matter of fact, both the size and growth rate time histories are qualitatively identical to those shown in figure 2, predicted by the analytical solution. The quantitative analysis of these experiments will be carried out in a future companion paper, as it involves taking into account the presence of the plate and natural convection, among other effects.
Another manifestation of the history effect is observed for
${\it\tau}>T_{6}^{+}$
, after the expansion–compression cycles have ended (see figure 2). Even though the pressure returns to the saturation value,
$c_{s}({\it\tau}\geqslant T_{6}^{+})=0$
, static equilibrium is not instantaneously achieved. As before, a supersaturation region near the bubble (
$c>0$
) remains from the preceding dissolution stage. This amounts to a positive interfacial concentration gradient. The radius of the bubble grows from
$a({\it\tau}=T_{6}^{+})\approx 0.9$
to
$a({\it\tau}=30)\approx 1.1$
, certainly a non-negligible increment. This growth is entirely provided by the history effect. In other words, the interfacial concentration gradient in (2.16) is non-zero due to the sole contribution of the history integral term.
4.2 Comparison with numerical simulation
An important question that arises now is how important are history effects compared to others that have been neglected in our discussion thus far, such as advection or surface tension. To shed light on this question, we numerically solve the full mass transfer problem (2.12), (2.13) and (3.2) which takes these effects into account. Details of the numerical treatment of the problem may be found in § B.2.
Figure 4. Evolution of (a) the dimensionless bubble radius
$a$
in time
${\it\tau}$
according to the simulation scenarios (i)–(iv) and the analytical solution in (B 3) and (B 4). The bubble is initially in a 105 % supersaturated liquid. It is exposed to (b) a prescribed pressure time history consisting of two particular jumps. Inset: the dissolution dynamics for
${\it\tau}>5$
according to the Epstein–Plesset solution (3.8) is compared with solutions (iii) and (iv). For reference, the physical parameters employed are those considering a
$\text{CO}_{2}$
bubble of size
$R_{c}=175~{\rm\mu}\text{m}$
in water (
$k_{H}=3.40\times 10^{-4}~\text{mol}~\text{N}^{-1}~\text{ m}^{-1}$
,
$D_{m}=1.92\times 10^{-9}~\text{m}^{2}~\text{s}^{-1}$
,
${\it\gamma}_{\mathit{lg}}=0.07~\text{N}~\text{m}^{-1}$
) under conditions
$P_{c}=2$
bar and
$T_{\infty }=293~\text{K}$
. The resulting surface tension, solubility and saturation parameters are
${\it\sigma}=0.004$
,
${\it\Lambda}=0.828$
and
${\it\Upsilon}=1.05$
.
As an example, let us consider a bubble in an initially 105 % supersaturated liquid. The initial concentration field of dissolved gas is assumed uniform. Setting
$P_{c}=P_{\infty }(0)$
and
$R_{c}=R(0)$
entails
$p_{0}=1$
,
$a_{0}=1$
and
${\it\Upsilon}=1.05$
. The bubble is then exposed to the pressure time history shown in figure 4(b). There are three constant-pressure stages, the first corresponding to liquid supersaturation (expansion stage), the second to undersaturation (compression stage) and the third to perfect saturation (
$p_{\mathit{sat}}={\it\Upsilon}=1.05$
). The set of physical conditions are specified in the caption of figure 4. Four different versions of the problem are numerically solved, namely
-
(i) neglecting the dimensionless advection term, i.e. setting
$\mathit{Pe}=0$
in (B 13); -
(ii) neglecting surface tension, i.e. setting
${\it\sigma}=0$
in (B 12) and (B 14); -
(iii) neglecting both the advection term and surface tension (the same scenario as the one analytically solved, i.e.
$\mathit{Pe}=0$
,
${\it\sigma}=0$
) and -
(iv) taking into account both the advection term and Laplace pressure.
These solutions are compared in figure 4(a), along with the analytical solution (for which
$\mathit{Pe}={\it\sigma}=0$
). History effects are critical in order to properly describe this particular dissolution (compression) stage (
$5<{\it\tau}<10$
). While growth is well determined by existing solutions (Epstein & Plesset Reference Epstein and Plesset1950; Scriven Reference Scriven1959), the dissolution experienced by this bubble is greatly affected by the low-concentration boundary layer left by the preceding growth stage. It stands to reason that for the initial instants of dissolution, during which history effects are most important, the dissolution rate observed is much faster than that obtained under the assumption of a historyless, uniform concentration field at
${\it\tau}=5^{+}$
. This is corroborated in the inset plot of figure 4, where we plot the Epstein–Plesset solution (3.8), which assumes a uniform initial concentration field (with
$p_{0}=1.1$
). To enable proper comparison, the initial radius
$a_{0}$
at
${\it\tau}=5^{+}$
for the Epstein–Plesset solution with and without surface tension has been fitted to that predicted by curve (iv) and (iii) respectively. Note that curve (iii) is identical to the analytical solution.
When the pressure drops to the saturation value at
${\it\tau}=10$
, the bubble continues to grow for some time. As discussed in § 4.1, it is consequence of the high concentration boundary layer of dissolved gas in the nearby liquid that was left by the preceding dissolution stage. Since the bubble is nearly close to the equilibrium in this stage, both advection and surface tension are small compared to history, at least for a few time units. At later times however, surface tension breaks the diffusive equilibrium, and ultimately drives the bubble towards its total dissolution. Nonetheless, this occurs at times much longer than those of the pressure cycle considered here.
Under the conditions investigated here, the effect of advection and surface tension on the bubble radius dynamics nearly cancel each other out during the expansion–compression cycle. Consequently, it is observed that the full solution (iv) (which contains both surface tension and advection) is in better agreement with the analytical solution (iii) (which contains neither) than with curves (i) or (iii) during this period. In fact, a quantitative explanation may be obtained from a simple analysis. For small values of
$\mathit{Pe}$
and constant ambient pressure, the contribution of advection towards the growth or dissolution rate may be quantified using the asymptotic solution of Duda & Vrentas (Reference Duda and Vrentas1969, equation (43)). To do so, we subtract the Epstein–Plesset solution (3.8) with
${\it\sigma}=0$
from (43). For our particular set of initial conditions,
$p_{0}=a_{0}=1$
, we obtain
$$\begin{eqnarray}\left[\frac{\text{d}a}{\text{d}{\it\tau}}\right]_{\mathit{Pe}\neq 0}-\left[\frac{\text{d}a}{\text{d}{\it\tau}}\right]_{\mathit{Pe}=0}\approx {\it\Lambda}^{2}\left(1-{\it\Upsilon}\right)^{2}\left[1-\frac{2}{{\rm\pi}}+\frac{1}{\sqrt{{\rm\pi}{\it\tau}}}\left({\it\tau}+\frac{2}{{\rm\pi}}-\frac{1}{2}\right)\right].\end{eqnarray}$$
Conversely, a simple estimation of the effect of the Laplace pressure on the rate of change of the bubble radius can be done by subtracting (3.8) with
${\it\sigma}=0$
from itself. Provided
$p_{0}=a_{0}=1$
and
${\it\sigma}\ll 1$
, this results in
$$\begin{eqnarray}\left[\frac{\text{d}a}{\text{d}{\it\tau}}\right]_{{\it\sigma}\neq 0}-\left[\frac{\text{d}a}{\text{d}{\it\tau}}\right]_{{\it\sigma}=0}\approx -\frac{{\it\Lambda}{\it\sigma}}{3a}(1+2{\it\Upsilon})\left(\frac{1}{a}+\frac{1}{\sqrt{{\rm\pi}{\it\tau}}}\right).\end{eqnarray}$$
We may use these formulas to evaluate the characteristic percentage contribution to the growth rate of each effect with respect to that given by Epstein–Plesset solution (3.8) with
${\it\sigma}=0$
. Halfway through the initial growth stage (
${\it\tau}=2.5$
,
$a=1.18$
) we determine that advection enhances growth by approximately 5 % whereas surface tension slows it down by
$-7\,\%$
. Thus, this small net difference of 2 % (which actually decreases as the bubble grows) justifies the close agreement between curves (iii) and (iv).
5 Small-amplitude isothermal oscillations
This section aims to deliver insight on the role of the history effect in the problem concerning the mass transfer across a bubble that pulsates non-inertially under sinusoidal acoustic excitation. This constitutes a particular regime under the broad phenomenon known as rectified diffusion, associated with many practical applications. For a more general description of this complex phenomenon the reader is referred to the seminal works of Hsieh & Plesset (Reference Hsieh and Plesset1961), Eller & Flynn (Reference Eller and Flynn1965), Crum & Hansen (Reference Crum and Hansen1982) or the more recent work by Zhang & Li (Reference Zhang and Li2014a ). Here, we will only consider low-frequency forcing, which has special relevance in seismological events (Tisato et al. Reference Tisato, Quintal, Chapman, Podladchikov and Burg2015) and in the exposure of marine mammals to low-frequency sonars (Crum & Mao Reference Crum and Mao1996). More specifically, we shall restrict the present analysis to small-amplitude, isothermal oscillations consistent with (2.4) and (2.5).
Provided the bubble pressure varies harmonically in time, it is then possible to readily determine the role of the history effect on the mass transfer across such a bubble. In turn, it may be used to explain and predict the nature (phase and amplitude) of oscillation. To this end, let us consider a harmonic pressure forcing,
$$\begin{eqnarray}P(t)=P_{c}[1+{\it\varepsilon}\sin (2{\rm\pi}f_{c}t)],\end{eqnarray}$$
where
${\it\varepsilon}$
is the dimensionless forcing amplitude and
$f_{c}$
is the forcing frequency. This time it shall prove convenient to work with the linear time
${\it\tau}=D_{m}t/R_{c}^{2}$
. The dimensionless pressure
$p({\it\tau})$
is then
$$\begin{eqnarray}p({\it\tau})=1+{\it\varepsilon}\sin {\it\Omega}{\it\tau},\quad \text{with }{\it\Omega}=\frac{2{\rm\pi}f_{c}R_{c}^{2}}{D_{m}}.\end{eqnarray}$$
The angular forcing frequency
${\it\Omega}$
has been non-dimensionlized with the diffusive time scale. The frequency
${\it\Omega}/(2{\rm\pi})$
is in fact equivalent to the frequency-based Péclet number typically used in rectified diffusion problems,
$\mathit{Pe}_{f}=f_{c}R_{c}^{2}/D_{m}$
, as introduced by Fyrillas & Szeri (Reference Fyrillas and Szeri1994). We propose a solution for the bubble dynamics of the form
$$\begin{eqnarray}a({\it\tau})=\bar{a}({\it\tau})-{\it\delta}\sin ({\it\Omega}{\it\tau}+{\it\phi})+O({\it\delta}^{2}),\end{eqnarray}$$
where
$\bar{a}=\bar{R}/R_{c}$
is the dimensionless equilibrium radius, while
${\it\delta}$
is the oscillation amplitude. Lastly,
${\it\phi}={\it\phi}({\it\Omega})$
denotes the phase shift in oscillations compared to the phase expected for an equivalent non-soluble, isothermally contracting and expanding bubble under the sole effect of the oscillatory ambient pressure forcing.
Solution (5.3) and the analysis soon to follow make use of three underlying assumptions.
-
(i) The bubble remains isothermal throughout the oscillations. Inertial and viscous effects in the bubble dynamics are completely neglected. The effect of the Laplace pressure on the oscillatory problem is also ignored:
${\it\sigma}/\bar{a}\ll 1$
. -
(ii) The strain amplitude of the oscillations is small,
${\it\delta}/\bar{a}\ll 1$
. -
(iii) The equilibrium radius
$\bar{a}({\it\tau})$
is assumed to vary sufficiently slowly in time to be treated as constant within an individual oscillation period. This time scale for bubble growth/dissolution must be much larger than the period of oscillation,
${\it\Omega}^{-1}$
.
The derivation and discussion of the range of frequency
${\it\Omega}$
, pressure amplitude
${\it\varepsilon}$
and saturation level
${\it\Upsilon}$
for which assumptions (i)–(iii) hold is provided in appendix C. Here, we choose to highlight that assumption (i) requires that the value of
${\it\Omega}$
must satisfy
$$\begin{eqnarray}{\it\Omega}<D_{\mathit{th}}/D_{m},\quad 1\lesssim {\it\Omega}\ll {\it\Omega}_{\mathit{res}},\end{eqnarray}$$
where
$D_{\mathit{th}}$
is the thermal diffusivity of the gas at constant volume and
${\it\Omega}_{\mathit{res}}$
is the dimensionless Minnaert resonance frequency. Note that this range may span several orders of magnitude. Moreover, assumption (ii) requires that
${\it\varepsilon}\ll 1$
.
Assumption (i) allows for the gas pressure in the bubble to be the same as the ambient pressure
$p({\it\tau})$
and consequently
$c_{s}({\it\tau})=p({\it\tau})-{\it\Upsilon}$
. The radius
$a_{p}$
of the equivalent non-soluble bubble mentioned above (for which
${\it\phi}_{p}\equiv 0$
) then satisfies
$pa_{p}^{3}=\bar{a}^{3}$
. Linearization results in
$$\begin{eqnarray}a_{p}({\it\tau})=\bar{a}p^{-1/3}=\bar{a}(1+{\it\varepsilon}\sin {\it\Omega}{\it\tau})^{-1/3}=\bar{a}-\frac{{\it\varepsilon}\bar{a}}{3}\sin {\it\Omega}{\it\tau}+O({\it\varepsilon}^{2}\bar{a}^{2}).\end{eqnarray}$$
Naturally, (i) implies that
$a_{p}$
is in anti-phase with pressure
$p$
and interfacial concentration
$c_{s}$
.
To conclude, we would like to remark that assumptions (i)–(iii) essentially simplify this complex mass transfer problem to the point that the role of the history effect can be smoothly rooted out and analysed by means of simple analytical expressions. This is the purpose of the next subsection.
We anticipate that it is beyond the scope of this work to determine the threshold pressure amplitude for growth under rectified diffusion or to model the bubble growth rate under specific low-frequency scenarios of interest. These notions have been investigated in the case of volcanic systems (Ichihara & Brodsky Reference Ichihara and Brodsky2006), and in high supersaturation levels with potentially large-amplitude oscillations, as observed in the capillaries of marine mammals (Crum & Mao Reference Crum and Mao1996; Ilinskii, Wilson & Hamilton Reference Ilinskii, Wilson and Hamilton2008) in which the effect of the viscoelastic medium is undoubtedly important (Zhang & Li Reference Zhang and Li2014b ).
5.1 The oscillatory problem
The full advection–diffusion problem is usually tackled by splitting it into a smooth problem and an oscillatory problem following the work of Fyrillas & Szeri (Reference Fyrillas and Szeri1994). The smooth problem (diffusion time scale
${\it\tau}\sim 1$
) accounts for the steady part of the boundary condition and yields the net flux of dissolved gas across the oscillating bubble. The oscillatory problem (diffusion time scale
${\it\tau}\sim {\it\Omega}^{-1}$
) takes into account the unsteady part of the boundary condition and describes the zero-average-mass exchange occurring over one bubble oscillation. This approach is only valid on the assumption that the time scales of the two processes are well separated:
${\it\Omega}\gg 1$
. Since we are also considering frequencies
${\it\Omega}\sim 1$
, splitting the problem is no longer possible. Moreover, the smooth solution, whether following the approach of Eller & Flynn (Reference Eller and Flynn1965) or Fyrillas & Szeri (Reference Fyrillas and Szeri1994), is constructed from the mass transport equation in Lagrangian spherical coordinates (as opposed to (5.6)) and requires special time averaging of the radius dynamics to properly capture the effect of advection. The interested reader is referred to Ilinskii et al. (Reference Ilinskii, Wilson and Hamilton2008) for a brief review on these two forms of the smooth solution.
Fortunately, the history effect is only directly relevant in the oscillatory problem associated with the unsteady boundary condition. In other words, we are only interested in the solution during an individual period of oscillation about
$\bar{a}$
. Solving the governing mass transport equation directly in the form
$$\begin{eqnarray}\frac{\partial c}{\partial {\it\tau}}+\frac{\mathit{Pe}({\it\tau})}{a^{2}}\left(\frac{1}{{\it\xi}^{2}}-{\it\xi}\right)\frac{\partial c}{\partial {\it\xi}}=\frac{1}{a^{2}{\it\xi}^{2}}\frac{\partial }{\partial {\it\xi}}\left({\it\xi}^{2}\frac{\partial c}{\partial {\it\xi}}\right)\end{eqnarray}$$
is suitable to this end. Note that
$\mathit{Pe}({\it\tau})=(\text{d}a/\text{d}{\it\tau})a\equiv R{\dot{R}}/D_{m}$
is our
${\dot{R}}$
-based Péclet number (as previously defined in § 2.2). It follows that under the small-amplitude oscillation restriction, the magnitude of the dimensionless advection term in (5.6) is always smaller than the unsteady and diffusive terms by a factor of
$O({\it\varepsilon})$
, where
${\it\varepsilon}\ll 1$
(see appendix D). Note that this condition holds true regardless of the magnitude of
$\mathit{Pe}$
(cf.
$O(\mathit{Pe})\sim {\it\varepsilon}{\it\Omega}\bar{a}^{2}/3$
), which may not necessarily be smaller than unity. It is therefore not unreasonable to neglect the advection term when dealing just with the oscillatory problem consisting of an individual cycle. Note, nonetheless, that the small contribution of the advection term on the interfacial gradient (negligible over one individual cycle) is still crucial to providing net growth over many cycles.
Moreover, assumption (ii) allows us to make the approximation
$a({\it\tau})=\bar{a}({\it\tau})+O({\it\delta})\approx \bar{a}$
, which is constant over an oscillation period according to (iii). While this greatly simplifies the problem, it only induces very small errors of
$O({\it\varepsilon})\ll 1$
.
Additionally, we will consider saturation conditions
${\it\Upsilon}=1$
. This is justified in that the bulk concentration of the liquid comes into play in the smooth solution only (Fyrillas & Szeri Reference Fyrillas and Szeri1994). Similarly, surface tension and time-averaged bubble interface motion are mostly relevant in the smooth solution.
Implementing these notions, the solution to (5.6) for the concentration profile is found to be (see appendix D for the complete derivation)
$$\begin{eqnarray}c({\it\xi},{\it\tau})=\frac{{\it\varepsilon}}{{\it\xi}}\exp \left\{-\bar{a}\sqrt{{\it\Omega}/2}({\it\xi}-1)\right\}\sin \!\left\{{\it\Omega}{\it\tau}-\bar{a}\sqrt{{\it\Omega}/2}({\it\xi}-1)\right\},\end{eqnarray}$$
which renders the following concentration gradient at the wall:
$$\begin{eqnarray}\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}=-{\it\varepsilon}\left[\sin {\it\Omega}{\it\tau}+\bar{a}\sqrt{{\it\Omega}}\sin \!\left({\it\Omega}{\it\tau}+\frac{{\rm\pi}}{4}\right)\right].\end{eqnarray}$$
The first term on the right-hand side of (5.8) refers to the quasi-steady mass flux imposed by the current value of
$c_{s}({\it\tau})={\it\varepsilon}\sin {\it\Omega}{\it\tau}$
. The last term in the right-hand side constitutes the contribution from the history integral (derived in appendix D). In the scenario of small-amplitude oscillations considered here, the mass flux coming from the contribution of the history term is thus observed to be proportional to the pressure amplitude
${\it\varepsilon}$
. The history effect imposes a phase shift in the sinusoidal interfacial concentration gradient,
$0<{\it\phi}_{\mathit{grad}}({\it\Omega})\leqslant {\rm\pi}/4$
, with respect to the phase of
$-c_{s}({\it\Omega}{\it\tau})$
or the oscillatory component of
$a_{p}({\it\tau})$
.
As
${\it\Omega}\rightarrow 0$
, the history integral term vanishes and
${\it\phi}_{\mathit{grad}}\rightarrow 0$
. The mass flux is in phase with
$a_{p}$
(cf. (5.5)). Conversely, as
${\it\Omega}\rightarrow \infty$
, the instantaneous mass flux is entirely provided by the history integral. The maximum possible phase shift with respect to
$a_{p}$
is attained:
${\it\phi}_{\mathit{grad}}={\rm\pi}/4$
. In fact, the frequency dependency of this phase shift is independent of
${\it\varepsilon}$
and may be shown to be
$$\begin{eqnarray}{\it\phi}_{\mathit{grad}}({\it\Omega},\bar{a})=-2\arctan \!\left\{1+\frac{\sqrt{2}}{\bar{a}\sqrt{{\it\Omega}}}-\frac{\sqrt{2}}{\bar{a}\sqrt{{\it\Omega}}}\left(\bar{a}^{2}{\it\Omega}+\sqrt{2{\it\Omega}}\bar{a}+1\right)^{1/2}\right\}.\end{eqnarray}$$
Solution (5.7) may be described as a damped transverse wave with wavenumber
$k=\bar{a}\sqrt{{\it\Omega}/2}$
and phase velocity
$\sqrt{2{\it\Omega}}/\bar{a}$
attenuating over a penetration depth
$l\sim 1/(\bar{a}\sqrt{{\it\Omega}})$
. Figure 5 shows the close agreement of this (advectionless) analytical solution with the full numerical computation (taking into account advection) involving (5.6) and the mass conservation equation (3.5) with no surface tension for
${\it\Omega}=1$
and
${\it\Omega}=1000$
(with equilibrium radius
$\bar{a}=1$
). This verifies the negligible qualitative impact of the advection term in the oscillatory problem even for this moderate pressure amplitude of
${\it\varepsilon}=0.01$
. On another note, the boundary layer thickness
$h$
may be best estimated as the wavelength,
$h=2{\rm\pi}/k$
. Over a distance
$h$
, (horizontal span in figure 5) the concentration amplitude reduces by a factor of
$\text{e}^{2{\rm\pi}}\approx 540$
.
Figure 5. Radial concentration profiles at different times over an oscillation period for (a)
${\it\Omega}=1$
(b)
${\it\Omega}=1000$
taking
$\bar{a}=1$
and
${\it\varepsilon}=0.01$
. The radial coordinate has been normalized by the wavenumber and spans one wavelength. The circle represents the relative bubble size (drawn to scale with the horizontal axis) for reference.
For
${\it\Omega}=1000$
, figure 5(b) hints that the shift in the phase of the concentration gradient with respect to
$-c_{s}$
is practically that of the high-frequency limit,
${\it\phi}_{\mathit{grad}}={\rm\pi}/4$
. The interfacial gradient is flat at
${\it\Omega}{\it\tau}=3{\rm\pi}/4,7{\rm\pi}/4$
as opposed to
$c_{s}({\it\tau})=0$
taking place at
${\it\Omega}{\it\tau}={\rm\pi},2{\rm\pi}$
. For
${\it\Omega}=1$
, direct evaluation of (5.9) yields
${\it\phi}_{\mathit{grad}}={\rm\pi}/8$
(as later observed in figure 6
a).
Figure 6. Numerically obtained oscillation waveforms for
${\it\Omega}=1$
, with
$\bar{a}=1$
,
${\it\varepsilon}=0.01$
and
${\it\Lambda}=0.828$
. Phase shifts are annotated in (a) showing the normalized waveforms for strains
$(a-\bar{a})/{\it\delta}$
,
$(a_{\mathit{corr}}-\bar{a})/{\it\delta}_{\mathit{corr}}$
,
$(a_{p}-\bar{a}_{p})/{\it\delta}_{p}$
and for the interfacial concentration gradient
$(\partial c/\partial {\it\xi}|_{{\it\xi}=1})/{\it\delta}_{\mathit{grad}}$
. (b) Unscaled strain waveforms
$a/\bar{a}-1$
,
$a_{\mathit{corr}}/\bar{a}-1$
and
$a_{p}/\bar{a}-1$
.
It stands to reason that solutions (5.7) and (5.8) evaluated in the limit
${\it\Omega}\gg 1$
converge to those one may obtain following the approach of Fyrillas & Szeri (Reference Fyrillas and Szeri1994) (inherently valid for
${\it\Omega}\gg 1$
) subject to the quasi-static radius approximation, i.e. taking
$a({\it\tau})=\bar{a}$
as constant. The complete derivation may be found in § D.1.
5.2 Strain amplitudes and phase
The frequency-dependent
${\it\phi}_{\mathit{grad}}$
induced by the history effect implies that pressure-induced contraction of the bubble radius (
$p>1$
,
$c_{s}>0$
,
$\text{d}a_{p}/\text{d}{\it\tau}<0$
) is not in phase with dissolution (negative mass diffusion rate,
$\partial c/\partial {\it\xi}(1,{\it\tau})>0$
). Likewise, the pressure-induced expansion is not in phase with growth. We would like to conclude this work by determining the dependency on
${\it\Omega}$
of the overall strain amplitude and phase lag with respect to
$a_{p}({\it\phi}<0)$
of the isothermally pulsating bubble.
Let us then define the dimensionless pressure-corrected radius,
$a_{\mathit{corr}}({\it\tau})$
, in an analogous manner to
$R_{\mathit{corr}}$
in (1.3) as follows:
$$\begin{eqnarray}a=p^{-1/3}a_{\mathit{corr}}=a_{p}a_{\mathit{corr}}/\bar{a}.\end{eqnarray}$$
The pressure-corrected radius is associated purely with the mass transfer across the interface. It has the useful property that its rate of change,
$\text{d}a_{\mathit{corr}}/\text{d}{\it\tau}$
, always takes the same sign as the interfacial concentration gradient. Consequently,
${\it\phi}_{\mathit{corr}}={\it\phi}_{\mathit{grad}}-{\rm\pi}/2$
, as shown in appendix E.
The characteristic amplitude of oscillation for the interfacial concentration gradient may be estimated from (5.8),
$$\begin{eqnarray}O\!\left(\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}\right)\sim {\it\varepsilon}(1+\bar{a}\sqrt{{\it\Omega}})={\it\delta}_{\mathit{grad}}.\end{eqnarray}$$
Equation (5.10) may be rewritten in terms of the different oscillation amplitudes about the equilibrium radius
$\bar{a}$
,
$$\begin{eqnarray}\bar{a}+{\it\delta}\approx (\bar{a}+{\it\delta}_{p})(\bar{a}+{\it\delta}_{\mathit{corr}})/\bar{a}.\end{eqnarray}$$
Linearization provides
${\it\delta}={\it\delta}_{p}+{\it\delta}_{\mathit{corr}}+O({\it\delta}_{p}{\it\delta}_{\mathit{corr}})$
. Amplitudes
${\it\delta}_{p}$
and
${\it\delta}_{\mathit{corr}}$
may be reasonably approximated as the product of the characteristic time derivative and the characteristic time scale,
$1/{\it\Omega}$
. We obtain
$$\begin{eqnarray}\frac{1}{{\it\Omega}}\frac{\text{d}a_{p}}{\text{d}{\it\tau}}\sim \frac{{\it\varepsilon}\bar{a}}{3}={\it\delta}_{p},\end{eqnarray}$$
which is consistent with (5.5). Likewise, making use of (E 3),
$$\begin{eqnarray}\frac{1}{{\it\Omega}}\frac{\text{d}a_{\mathit{corr}}}{\text{d}{\it\tau}}\sim \frac{{\it\Lambda}{\it\varepsilon}}{\bar{a}{\it\Omega}}\left(1+\bar{a}\sqrt{{\it\Omega}}\right)={\it\delta}_{\mathit{corr}}.\end{eqnarray}$$
Two limiting cases arise. For small
${\it\Omega}\ll 1$
, we see that
${\it\delta}\approx {\it\delta}_{\mathit{corr}}\sim {\it\varepsilon}{\it\Lambda}/\bar{a}{\it\Omega}$
. The bubble oscillation amplitudes are provided entirely by mass transfer across the bubble surface, the pressure-induced expansion and contraction amplitude is negligible in comparison. The behaviour of the different phase shifts as
${\it\Omega}$
approaches this limit is:
$$\begin{eqnarray}\text{as}\quad {\it\Omega}\rightarrow 0,\quad {\it\phi}_{\mathit{grad}}\rightarrow 0,\quad {\it\phi}_{\mathit{corr}}\rightarrow -{\rm\pi}/2,\quad {\it\phi}\rightarrow {\it\phi}_{\mathit{corr}}\rightarrow -{\rm\pi}/2.\end{eqnarray}$$
Taking
${\it\Omega}=1$
,
$\bar{a}=1$
(figure 6), the dominant contribution in
${\it\delta}$
and hence
${\it\phi}$
still come from
${\it\delta}_{\mathit{corr}}$
and
${\it\phi}_{\mathit{corr}}$
. The relative contributions of
${\it\delta}_{\mathit{corr}}$
and
${\it\delta}_{p}$
may be estimated from
${\it\delta}_{\mathit{corr}}/{\it\delta}_{p}\sim 6{\it\Lambda}\approx 5$
.
For large
${\it\Omega}\gg 1$
, we see that
${\it\delta}\approx {\it\delta}_{p}\sim {\it\varepsilon}\bar{a}/3$
. There is negligible mass transfer during an individual oscillation due to the short oscillation period. The phase shifts behave according to:
$$\begin{eqnarray}\text{as}\quad {\it\Omega}\rightarrow \infty ,\quad {\it\phi}_{\mathit{grad}}\rightarrow {\rm\pi}/4,\quad {\it\phi}_{\mathit{corr}}\rightarrow -{\rm\pi}/4,\quad {\it\phi}\rightarrow {\it\phi}_{p}=0.\end{eqnarray}$$
Taking
${\it\Omega}=100$
,
$\bar{a}=1$
(figure 7),
$a_{p}$
now provides the main contribution since
${\it\delta}_{\mathit{corr}}/{\it\delta}_{p}\sim 3{\it\Lambda}/10\approx 1/4$
.
Figure 7. Oscillation waveforms for
${\it\Omega}=100$
, with
$\bar{a}=1$
,
${\it\varepsilon}=0.01$
(see caption of figure 6).
Ultimately, the difference in the phase lags
${\it\phi}$
and
${\it\phi}_{\mathit{corr}}$
partly induced by the history effect has a direct impact on the relative contributions of the area effect and the so-called shell effect (Eller & Flynn Reference Eller and Flynn1965), the two driving mechanisms behind bubble growth in rectified diffusion. The area effect refers to the increase of the bubble interfacial area during expansion and the subsequent decrease during contraction. The velocity field generated by the bubble oscillations also influences the gas transport through the advection term in the advection–diffusion equation. This influence is referred to as the shell effect. More specifically, we can regard the shell effect as the advection-induced squeezing of the concentration boundary layer when the bubble radius expands and similarly the stretching when the radius contracts. In turn, this effect yields either a steeper or shallower concentration gradient at the bubble surface compared to that of a motionless bubble of the same size. As a result, the area effect increases mass transfer with respect to the motionless equilibrium bubble when
$a>\bar{a}$
. Similarly, the shell effect has an amplifying effect on the interfacial mass flux when
$\text{d}a/\text{d}{\it\tau}>0$
.
Figure 8. Strain waveforms for
${\it\Omega}=5000$
, with
$\bar{a}=1$
,
${\it\varepsilon}=0.01$
and
${\it\Lambda}=0.828$
. (a) Normalized strain waveforms
$(a-\bar{a})/{\it\delta}$
and
$(a_{\mathit{corr}}-\bar{a})/{\it\delta}_{\mathit{corr}}$
. SE and AE refer to the ‘shell effect’ and ‘area effect’, which may increase (
$+$
) or hinder (
$-$
) mass transfer across the bubble, be it growth (
$\text{d}a_{\mathit{corr}}/\text{d}{\it\tau}>0$
, grey background) or dissolution. (b) Unscaled strain waveforms
$a/\bar{a}-1$
and
$a_{\mathit{corr}}/\bar{a}-1$
.
Let us consider the high-frequency limit where
${\it\phi}=0$
,
${\it\phi}_{\mathit{corr}}={\rm\pi}/4$
and
$a\approx a_{p}$
. The area and shell effects act non-symmetrically in the bubble dissolution (
$\text{d}a_{\mathit{corr}}/\text{d}{\it\tau}<0$
) and growth (
$\text{d}a_{\mathit{corr}}/\text{d}{\it\tau}>0$
) periods, which are of equal duration in time. Figure 8 shows that the area and shell effect increase mass transfer during three-quarters of the growth period and one-quarter of the dissolution period. Consequently, these effects reduce mass transfer during one-quarter of the growth period and three-quarters of the dissolution period. As a result, bubble growth is always promoted.
6 Conclusions
The contribution of any past mass transfer events between a gas bubble and its liquid surroundings towards the current diffusion-driven bubble growth or dissolution dynamics has been referred to as the ‘history effect’. The history effect arises from the non-instantaneous development of the concentration boundary layer in response to changes in the concentration at the bubble interface caused, for instance, by variations of the ambient pressure in time. To put in another way, the current state of the concentration profile must be naturally conditioned by the preceding time history of the concentration field itself. As a consequence, the mass flux across the bubble is conditioned by the preceding time history of the boundary condition. This very last notion is the essence of the history effect. It has been shown that the contribution of the history effect in the current interfacial concentration gradient is fully contained within a memory integral of the interface concentration.
Under the assumption that advection effects are small, the integral term has been used to derive the governing equation for bubble dynamics concerning the canonical case of a spherical bubble suspended in a quiescent liquid. It has been termed as the Epstein–Plesset with history term (EPH) equation. This equation is not restricted to a constant-pressure time history and does not make use of the quasi-static radius approximation. It is however expressed in nonlinear time
$\tilde{{\it\tau}}$
, which can be related to the standard time once the EPH equation has been solved for
$a(\tilde{{\it\tau}})$
, the time evolution of the bubble radius.
The EPH equation has been analytically solved for the case of multiple step-like jumps in pressure. The nature and relevance of the history effect in the bubble dynamics has been assessed through illustrative examples. A future companion paper shall deal with the experimental and numerical analysis of the history effect regarding a sessile bubble exposed to such step-like pressure time histories.
Finally, we have investigated the role of the history effect concerning the problem of mass transfer across a non-inertial bubble that pulsates under harmonic pressure forcing. The history effect has been shown to induce a phase shift in the interfacial concentration gradient with respect to the phase of the interfacial concentration. At sufficiently high forcing frequencies, the oscillatory mass flux across the bubble is entirely provided by the history integral term and the aforementioned phase shift asymptotically tends to
${\rm\pi}/4$
. This phase shift causes the shell effect and area effect (the two driving mechanisms behind rectified diffusion) to act non-symmetrically on growth and dissolution during each individual bubble oscillation period. Growth is always favoured as a result.
Acknowledgements
The authors acknowledge the support of the Spanish Ministry of Economy and Competitiveness through grant DPI2014-59292-C3-1-P, partly financed by European funds. The authors are deeply grateful to A. Moreno-Soto especially for his experimental guidance during the acquisition of figure 1. The authors are also grateful to O. Enríquez for promoting this study and D. Lohse for his helpful suggestions.
Appendix A. Derivation of the history integral term
Upon the substitution
$f({\it\xi},\tilde{{\it\tau}})={\it\xi}c({\it\xi},\tilde{{\it\tau}})$
, and furthermore neglecting any arising
$(\text{d}a/\text{d}\tilde{{\it\tau}})/a$
terms, equation (2.15) transforms to
$$\begin{eqnarray}\frac{\partial f}{\partial \tilde{{\it\tau}}}=\frac{\partial ^{2}f}{\partial {\it\xi}^{2}}.\end{eqnarray}$$
Taking the Laplace transform of (A 1) subject to the initial condition
$f({\it\xi}>1,0)=0$
yields
$$\begin{eqnarray}s\hat{f}=\frac{\text{d}^{2}\hat{f}}{\text{d}{\it\xi}^{2}}.\end{eqnarray}$$
A solution compatible with boundary conditions
$f(1,{\it\tau})=c_{s}(\tilde{{\it\tau}})$
and
$f(\infty ,\tilde{{\it\tau}})=0$
is found to be
$$\begin{eqnarray}\hat{f}({\it\xi};s)={\hat{c}}_{s}\text{e}^{-\sqrt{s}({\it\xi}-1)},\end{eqnarray}$$
where
${\hat{c}}_{s}$
is a function of
$s$
, to be determined from the boundary conditions. The transformed concentration gradient across the bubble interface is
$$\begin{eqnarray}\left.\frac{\partial {\hat{c}}}{\partial {\it\xi}}\right|_{{\it\xi}=1}=\left.\frac{\partial \hat{f}}{\partial {\it\xi}}\right|_{{\it\xi}=1}-\hat{f}(1;s)=-(1+\sqrt{s}){\hat{c}}_{s}.\end{eqnarray}$$
The initial value of
$c_{s}$
is required. It is given by
$$\begin{eqnarray}c_{s}(0)=p_{0}+{\it\sigma}/a_{0}-{\it\Upsilon}=c_{s0},\end{eqnarray}$$
where
$p_{0}=p(0)$
and
$a_{0}=a(0)$
denote the initial pressure and radius respectively. Note that in perfect initial saturation conditions, the surface concentration is the same as the gas concentration in the liquid,
$C_{s}(0)=C_{\infty }$
. This amounts to
${\it\Upsilon}=p_{0}$
and
$c_{s0}={\it\sigma}/a_{0}>0$
, corresponding to slow dissolution purely driven by the Laplace pressure. Taking then the inverse Laplace transform of (A 4) results in
$$\begin{eqnarray}-\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}=c_{s}+\int _{0}^{\tilde{{\it\tau}}}\frac{1}{\sqrt{{\rm\pi}(\tilde{{\it\tau}}-\tilde{x})}}\left[\frac{\text{d}c_{s}}{\text{d}\tilde{x}}+c_{s0}{\it\delta}(\tilde{x})\right]\text{d}\tilde{x}.\end{eqnarray}$$
The term containing the Dirac delta
${\it\delta}(\tilde{x})$
may be directly integrated. This finally sets the concentration gradient at the interface to be
$$\begin{eqnarray}-\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}=c_{s}+\frac{c_{s0}}{\sqrt{{\rm\pi}\tilde{{\it\tau}}}}+\int _{0}^{\tilde{{\it\tau}}}\frac{1}{\sqrt{{\rm\pi}(\tilde{{\it\tau}}-\tilde{x})}}\frac{\text{d}c_{s}}{\text{d}\tilde{x}}\,\text{d}\tilde{x}.\end{eqnarray}$$
It is not surprising that Jones & Zuber (Reference Jones and Zuber1978) obtained an analogous expression for the heat flux across a vapour bubble. In fact, the history integral provides the general solution (for the spatial derivative at the boundary) for the canonical diffusion problem consisting of a scalar field
$U(z,t)$
obeying
$$\begin{eqnarray}\frac{\partial U}{\partial t}={\it\alpha}\frac{\partial ^{2}U}{\partial z^{2}},\end{eqnarray}$$
with boundary and initial conditions
$U(0,t)=U_{s}(t)$
and
$U(z,0)=U(\infty ,t)=0$
. The solution for the gradient at the boundary
$z=0$
(see, for example, Landau & Lifshitz Reference Landau and Lifshitz1987) reads
$$\begin{eqnarray}-\left.\frac{\partial U}{\partial z}\right|_{z=0}=\frac{1}{\sqrt{{\it\alpha}}}\int _{0}^{t}\frac{1}{\sqrt{{\rm\pi}(t-x)}}\frac{\text{d}U_{s}}{\text{d}x}\,\text{d}x.\end{eqnarray}$$
The history integral is therefore essential in determining the heat flux across a plate or the shear stress exerted by a viscous fluid on a plate moving in its plane (often referred to as the Basset force). While the history integral has been identified to play an important role in the rectilinear motion of bubbles in viscous flows (Magnaudet & Legendre Reference Magnaudet and Legendre1998), its effect has been overlooked in diffusion-driven bubble dissolution and growth.
Appendix B. Solutions for step changes in pressure
B.1 Analytical solution
Let
$p_{n}$
be the pressure at the
$n\text{th}$
plateau (constant-pressure segment) after jump
$n$
. Its value is thus
$$\begin{eqnarray}p_{n}=p_{0}+\mathop{\sum }_{j=1}^{n}{\rm\Delta}p_{j}\quad \text{for }\tilde{T}_{n}^{+}\leqslant \tilde{{\it\tau}}\leqslant \tilde{T}_{n+1}^{-}.\end{eqnarray}$$
The time coefficient
$\tilde{T}_{n}^{-}$
describes the instant in time right before the
$n\text{th}$
pressure jump at
$\tilde{{\it\tau}}=\tilde{T}_{n}$
, while
$\tilde{T}_{n}^{+}$
refers to the moment right after.
B.1.1 Solution for the Epstein–Plesset with history term equation
The EPH equation (3.4) may be integrated in time within the
$n\text{th}$
segment as
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \int _{a(\tilde{T}_{n}^{-})}^{a(\tilde{{\it\tau}})}\,\text{d}\ln a+\frac{1}{3}\int _{p_{n-1}}^{p_{n}}\,\text{d}\ln p=-\frac{{\it\Lambda}}{p_{n}}\left[\int _{\tilde{T}_{n}^{+}}^{\tilde{{\it\tau}}}\left(p_{n}-{\it\Upsilon}+\frac{p_{0}-{\it\Upsilon}}{\sqrt{{\rm\pi}{\tilde{y}}}}\right)\,\text{d}{\tilde{y}}\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\left.\displaystyle \int _{\tilde{T}_{n}^{+}}^{\tilde{{\it\tau}}}\left\{\int _{0}^{{\tilde{y}}}\frac{1}{\sqrt{{\rm\pi}({\tilde{y}}-\tilde{x})}}\mathop{\sum }_{j=1}^{n}{\rm\Delta}p_{j}{\it\delta}(\tilde{x}-\tilde{T}_{j})\,\text{d}\tilde{x}\right\}\text{d}{\tilde{y}}\right].\end{eqnarray}$$
Evaluating this integral finally yields, for segments
$n=1$
to
$n=N$
,
$$\begin{eqnarray}\displaystyle a(\tilde{{\it\tau}}) & = & \displaystyle \displaystyle a(\tilde{T}_{n}^{-})\left(\frac{p_{n}}{p_{n-1}}\right)^{-1/3}\exp \left\{-\frac{{\it\Lambda}}{p_{n}}\left[(\tilde{{\it\tau}}-\tilde{T}_{n})(p_{n}-{\it\Upsilon})+\frac{2}{\sqrt{{\rm\pi}}}(p_{0}-{\it\Upsilon})\right.\right.\nonumber\\ \displaystyle & & \displaystyle \displaystyle \times \left(\sqrt{\tilde{{\it\tau}}}-\sqrt{\tilde{T}_{n}}\right)+\left.\left.\frac{2}{\sqrt{{\rm\pi}}}\mathop{\sum }_{j=1}^{n}{\rm\Delta}p_{j}\left(\sqrt{\tilde{{\it\tau}}-\tilde{T}_{j}}-\sqrt{\tilde{T}_{n}-\tilde{T}_{j}}\right)\right]\right\}\nonumber\\ \displaystyle & & \displaystyle \qquad \text{for }\tilde{T}_{n}^{+}\leqslant \tilde{{\it\tau}}\leqslant \tilde{T}_{n+1}^{-}.\end{eqnarray}$$
The initial condition is
$a(0)=a_{0}$
. The end radius from any segment
$a(\tilde{T}_{n}^{-})$
must be first computed before moving on to the next. We begin with segment
$n=0$
, i.e. before the first pressure jump. The radius dynamics in this initial segment is evidently identical to the solution in (3.10) and are straightforwardly given by
$$\begin{eqnarray}a(\tilde{{\it\tau}})=a_{0}\exp \{-{\it\Lambda}(p_{0}-{\it\Upsilon})(\tilde{{\it\tau}}+2\sqrt{\tilde{{\it\tau}}/{\rm\pi}})\}\quad \text{for }0\leqslant \tilde{{\it\tau}}\leqslant T_{1}^{-}.\end{eqnarray}$$
B.1.2 Solution for the concentration field
From (A 3), the Laplace transformed concentration field
${\hat{c}}({\it\xi};s)$
may be conveniently split as the product of two separate functions
$\hat{f}(s)$
and
${\hat{g}}({\it\xi};s)$
as follows:
$$\begin{eqnarray}{\hat{c}}({\it\xi};s)={\hat{g}}(s)\times {\hat{h}}({\it\xi};s)=s{\hat{c}}_{s}\times \frac{\text{e}^{-\sqrt{s}({\it\xi}-1)}}{{\it\xi}s}.\end{eqnarray}$$
The inverse Laplace transforms of these two functions may be shown to be
$$\begin{eqnarray}g(\tilde{{\it\tau}})=\frac{\text{d}c_{s}}{\text{d}{\it\tau}}+c_{s0}{\it\delta}(\tilde{{\it\tau}}),\quad h(\tilde{{\it\tau}})=\frac{1}{{\it\xi}}\text{erfc}\left(\frac{{\it\xi}-1}{2\sqrt{\tilde{{\it\tau}}}}\right).\end{eqnarray}$$
Using these results, the concentration field in the time domain
$c({\it\xi},\tilde{{\it\tau}})$
may be then computed by means of the convolution theorem. This gives
$$\begin{eqnarray}c({\it\xi},\tilde{{\it\tau}})=\frac{c_{s0}}{{\it\xi}}\text{erfc}\left(\frac{{\it\xi}-1}{2\sqrt{\tilde{{\it\tau}}}}\right)+\frac{1}{{\it\xi}}\int _{0}^{\tilde{{\it\tau}}}\frac{\text{d}c_{s}}{\text{d}\tilde{x}}\text{erfc}\left(\frac{{\it\xi}-1}{2\sqrt{\tilde{{\it\tau}}-\tilde{x}}}\right)\,\text{d}\tilde{x}.\end{eqnarray}$$
Surface tension is to be neglected once again. This renders
$c_{s0}=p_{0}-{\it\Upsilon}$
and
$\text{d}c_{s}/\text{d}\tilde{{\it\tau}}=\text{d}p/\text{d}\tilde{{\it\tau}}$
, the latter modelled in (4.1). After inserting these expressions into (B 7), the time evolution of the concentration field at every
$n\text{th}$
segment is finally obtained:
$$\begin{eqnarray}c({\it\xi},\tilde{{\it\tau}})=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{p_{0}-{\it\Upsilon}}{{\it\xi}}\text{erfc}\left(\frac{{\it\xi}-1}{2\sqrt{\tilde{{\it\tau}}}}\right) & \text{for }0\leqslant \tilde{{\it\tau}}\leqslant T_{1}^{-},\\ \displaystyle \frac{p_{0}-{\it\Upsilon}}{{\it\xi}}\text{erfc}\left(\frac{{\it\xi}-1}{2\sqrt{\tilde{{\it\tau}}}}\right)+\frac{1}{{\it\xi}}\mathop{\sum }_{j=1}^{n}{\rm\Delta}p_{j}\,\text{erfc}\left(\frac{{\it\xi}-1}{2\sqrt{\tilde{{\it\tau}}-\tilde{T}_{j}}}\right) & \text{for }\tilde{T}_{n}^{+}\leqslant \tilde{{\it\tau}}\leqslant \tilde{T}_{n+1}^{-}.\end{array}\right.\end{eqnarray}$$
B.2 Numerical model
The pressure step functions are analytically modelled by the logistic function
$$\begin{eqnarray}p(\tilde{{\it\tau}})=p_{0}+\frac{1}{2}\mathop{\sum }_{n=1}^{N}{\rm\Delta}p_{n}\{1+\tanh [\tilde{k}(\tilde{{\it\tau}}-\tilde{T}_{n})]\},\end{eqnarray}$$
where
$\tilde{k}$
is a large constant. The Heaviside functions are recovered as
$\tilde{k}\rightarrow \infty$
, whilst here the value
$\tilde{k}=1000$
was deemed as sufficiently large. The time derivative of the pressure reads
$$\begin{eqnarray}\frac{\text{d}p}{\text{d}\tilde{{\it\tau}}}=\frac{\tilde{k}}{2}\mathop{\sum }_{n=1}^{N}{\rm\Delta}p_{n}\text{sech}^{2}[\tilde{k}(\tilde{{\it\tau}}-\tilde{T}_{n})].\end{eqnarray}$$
The governing equations will be numerically integrated in nonlinear time
$\tilde{{\it\tau}}$
, purely for consistency with the analytical derivation. However, we shall conveniently establish the pressure time history input and present the bubble size history in the linear time
${\it\tau}$
. To this end, the physical time
${\it\tau}$
was computed at each time step by integrating
$\text{d}{\it\tau}=a^{2}\text{d}\tilde{{\it\tau}}$
. Provided
$\tilde{k}$
is large enough, (B 9) and (B 10) may be computed from
$$\begin{eqnarray}p(\tilde{{\it\tau}}-\tilde{T}_{n})=p({\it\tau}-T_{n})=p_{0}+\frac{1}{2}\mathop{\sum }_{n=1}^{N}{\rm\Delta}p_{n}\{1+\tanh [k({\it\tau}-T_{n})]\},\quad \frac{\text{d}p}{\text{d}\tilde{{\it\tau}}}=a^{2}\frac{\text{d}p}{\text{d}{\it\tau}},\end{eqnarray}$$
where use of
$k=\tilde{k}$
has been made. Note that the precise value of
$k$
is not relevant provided it is large. The radius dynamics can then be obtained by integrating (3.2), namely
$$\begin{eqnarray}\frac{\text{d}a}{\text{d}\tilde{{\it\tau}}}=a\left({\it\Lambda}\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}-\frac{1}{3}\frac{\text{d}p}{\text{d}\tilde{{\it\tau}}}\right)\left(p+\frac{2{\it\sigma}}{3a}\right)^{-1},\end{eqnarray}$$
subject to the prescribed initial size
$a_{0}$
. The concentration gradient at the interface is solved for numerically from (2.12) using a finite-differences scheme. Equation (2.12) is written here again for the reader’s convenience:
$$\begin{eqnarray}\frac{\partial c}{\partial \tilde{{\it\tau}}}+\mathit{Pe}(\tilde{{\it\tau}})\left(\frac{1}{{\it\xi}^{2}}-{\it\xi}\right)\frac{\partial c}{\partial {\it\xi}}=\frac{1}{{\it\xi}^{2}}\frac{\partial }{\partial {\it\xi}}\left({\it\xi}^{2}\frac{\partial c}{\partial {\it\xi}}\right),\end{eqnarray}$$
where
$\mathit{Pe}(\tilde{{\it\tau}})=(\text{d}a/\text{d}\tilde{{\it\tau}})/a$
. The required boundary and initial conditions are
$$\begin{eqnarray}c(1,\tilde{{\it\tau}})=p(\tilde{{\it\tau}})+{\it\sigma}/a-{\it\Upsilon},\quad \left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=\infty }=0,\quad c({\it\xi}>1,0)=0.\end{eqnarray}$$
Appendix C. Conditions for small-amplitude isothermal oscillations
Following the work of Prosperetti (Reference Prosperetti1977) (see figure 1 of that paper), the assumption of isothermal oscillations (polytropic exponent equal to unity) may be safely assumed, provided the thermal Péclet number based on the oscillation frequency is smaller than one:
${\it\Omega}_{\mathit{th}}\equiv 2{\rm\pi}f_{c}R_{c}^{2}/D_{\mathit{th}}<1$
, with
$D_{\mathit{th}}=k_{g}/{\it\rho}_{g}c_{v,g}$
. Symbols
$k_{g}$
,
${\it\rho}_{g}$
and
$c_{v,g}$
denote the gas thermal conductivity, density and specific heat at constant volume respectively. Consequently,
${\it\Omega}$
does not need to be small for this assumption to hold, as long as
${\it\Omega}<D_{\mathit{th}}/D_{m}$
. As an example, for air at 300 K,
$D_{\mathit{th}}/D_{m}\sim 10^{4}$
.
Assumption (i) in § 5 additionally implies that the oscillation frequency
$f_{c}$
must be much smaller than the resonance frequency of the bubble. The low-frequency limit, for history effects to be visible, is given by
$f_{c}\gtrsim 1/t_{m}$
, where
$t_{m}=R_{c}^{2}/D_{m}$
is the time scale of mass transfer by diffusion. When
$f_{c}\ll 1/t_{m}$
, the rate of change of the bubble interfacial concentration is very slow compared to mass diffusion. The concentration field has enough time to reach the quasi-steady solution characterized by the absence (full dissipation) of the history effect. Summarizing these last ideas,
${\it\Omega}$
must additionally satisfy
$$\begin{eqnarray}1\lesssim {\it\Omega}\ll {\it\Omega}_{\mathit{res}},\end{eqnarray}$$
where
${\it\Omega}_{\mathit{res}}=\sqrt{3P_{c}/{\it\rho}_{l}}R_{c}/D_{m}$
, is the Péclet number based on Minnaert’s resonance frequency. As an example, consider bubbles trapped in magma chambers exposed to earthquake-induced shaking with typical periods of 1–20 s (Ichihara & Brodsky Reference Ichihara and Brodsky2006). This results in
${\it\Omega}\sim 1$
–
$100$
for a
$100~{\rm\mu}\text{m}$
bubble. Likewise, a 10
${\rm\mu}\text{m}$
bubble in the tissue of a marine mammal exposed to 5–5000 Hz sonar (Crum & Mao Reference Crum and Mao1996) corresponds to
${\it\Omega}\sim 30$
–
$3000$
.
Except for small values of
${\it\Omega}\ll 1$
outside our range of interest (where the oscillation period is slow enough so that gas diffusion across the bubble surface results in large volumetric oscillations)
${\it\delta}$
will always be of the same order of magnitude as the pressure-induced isothermal expansion and contraction amplitude. Comparison of (5.3) and (5.5) gives
${\it\delta}/\bar{a}\sim {\it\varepsilon}/3$
. Therefore, (ii) requires that
${\it\varepsilon}\ll 1$
.
Turning now to (iii), consider a bubble in a gas-saturated liquid
${\it\Upsilon}=1$
, with negligible Laplace pressure
${\it\sigma}=0$
. Under these conditions, the growth rate purely due to rectified diffusion is expected to be largest. For such a case, the asymptotic growth rate of
$\bar{a}$
due to rectified diffusion under assumptions (i)–(iii) was found to be reasonably well approximated by Hsieh & Plesset (Reference Hsieh and Plesset1961),
$$\begin{eqnarray}\text{d}\bar{a}/\text{d}{\it\tau}=(2/3){\it\Lambda}{\it\varepsilon}^{2}.\end{eqnarray}$$
The inverse of the characteristic time scale for growth is then
${\it\tau}_{b}^{-1}\sim {\it\Lambda}{\it\varepsilon}^{2}/\bar{a}$
, which must be much smaller than
${\it\Omega}$
. Let us now consider low forcing frequencies in addition to a highly supersaturated (say
${\it\Upsilon}\sim 2$
–3) or undersaturated (
${\it\Upsilon}\approx 0$
) liquid together with a non-negligible Laplace pressure. The likely existence of high diffusion rates will result in fast diffusion-driven bubble growth or dissolution. Consequently, the equilibrium radius
$\bar{a}$
may no longer be constant over an individual oscillation period (Ilinskii et al.
Reference Ilinskii, Wilson and Hamilton2008), i.e. assumption (iii) no longer holds. The inverse of the characteristic time for bubble growth or dissolution may be estimated from (3.8) as
${\it\tau}_{b}^{-1}\sim |{\it\Upsilon}-1-{\it\sigma}/\bar{a}|{\it\Lambda}/\bar{a}^{2}$
. Thus, (iii) imposes that the solution saturation level, pressure amplitude and frequency must fulfil the following inequalities:
$$\begin{eqnarray}|{\it\Upsilon}-1-{\it\sigma}/\bar{a}|\ll {\it\Omega}\bar{a}^{2}/{\it\Lambda},\quad {\it\varepsilon}^{2}\ll \bar{a}{\it\Omega}/{\it\Lambda}.\end{eqnarray}$$
Finally note that if surface tension and undersaturation are to be overcome, equating Hsieh–Plesset solution (C 2) with the Epstein–Pesset solution (3.8) leads to an approximate, simplistic threshold condition for growth (Safar Reference Safar1968):
${\it\varepsilon}^{2}=(3/2)(1-{\it\Upsilon}+{\it\sigma}/\bar{a})$
. As noted by Safar (Reference Safar1968), this threshold is only valid for large isothermal bubbles in a liquid close to saturation (
${\it\Upsilon}\approx 1$
) under sufficiently high ambient pressures such that the Laplace pressure is comparatively small:
${\it\sigma}/\bar{a}\ll 1$
. As an example,
${\it\varepsilon}\sim 0.01$
would be required to overcome the surface-tension-driven dissolution of a
$\text{CO}_{2}$
bubble of size
$R_{c}\sim 100{\rm\mu}\text{m}$
in saturated water at
$P_{c}\sim 20~\text{MPa}$
.
Appendix D. The oscillatory problem
Valuable insight on the nature of the concentration field may be gained by performing an order of magnitude analysis on the governing mass transport equation:
$$\begin{eqnarray}\frac{\partial c}{\partial {\it\tau}}+\frac{\mathit{Pe}({\it\tau})}{a^{2}}\left(\frac{1}{{\it\xi}^{2}}-{\it\xi}\right)\frac{\partial c}{\partial {\it\xi}}=\frac{1}{a^{2}{\it\xi}^{2}}\frac{\partial }{\partial {\it\xi}}\left({\it\xi}^{2}\frac{\partial c}{\partial {\it\xi}}\right),\end{eqnarray}$$
where
$\mathit{Pe}({\it\tau})=(\text{d}a/\text{d}{\it\tau})a=R{\dot{R}}/D_{m}$
. Its magnitude is given by
$O(\mathit{Pe})\sim {\it\varepsilon}{\it\Omega}\bar{a}^{2}/3$
since
$O(a)\sim \bar{a}$
and
$O(\text{d}a/\text{d}{\it\tau})=O(\text{d}a_{p}/\text{d}{\it\tau})\sim {\it\varepsilon}{\it\Omega}\bar{a}/3$
. Bear in mind that
$\bar{a}$
will often be of order unity (provided
$R(t)$
is comparable to the chosen characteristic radius
$R_{c}$
), but we shall carry the analysis allowing for any magnitude of
$\bar{a}$
.
Evaluation of (D 1) on
${\it\xi}=1$
makes the advection term become identically zero. The diffusive term must therefore be of leading order in a layer bounded by
$1<{\it\xi}<1+l$
where
$O(\partial c)=O({\rm\Delta}c_{s})\sim {\it\varepsilon}$
. Parameter
$l$
is the dimensionless penetration depth of diffusion using
$R(t)$
as the length scale, i.e.
$l=L/R(t)$
. The relevant length scales are
$O({\it\xi})\sim 1$
and
$O(\partial {\it\xi})\sim l$
, while the characteristic time scale per oscillation is
$O(\partial {\it\tau})\sim {\it\Omega}^{-1}$
. The magnitudes of the unsteady and diffusive terms in (D 1) are given by
$$\begin{eqnarray}O\left(\frac{\partial c}{\partial {\it\tau}}\right)\sim {\it\varepsilon}{\it\Omega},\quad O\left(\frac{1}{a^{2}{\it\xi}^{2}}\frac{\partial }{\partial {\it\xi}}\left({\it\xi}^{2}\frac{\partial c}{\partial {\it\xi}}\right)\right)\sim \frac{{\it\varepsilon}}{\bar{a}^{2}l^{2}}.\end{eqnarray}$$
Balancing these two terms yields
$l\sim 1/(\bar{a}\sqrt{{\it\Omega}})$
. Note that the equivalent dimensional penetration depth is
$L\sim \sqrt{D_{m}/(2{\rm\pi}f_{c})}$
, and it is independent of the bubble size. In the low-frequency limit, it follows from the Epstein–Plesset solution that
$l\sim 1$
. Hence,
$(\bar{a}^{2}{\it\Omega})^{-1}\sim 1$
is of the order of the diffusion time scale associated with a length scale equal to the current bubble size. Length scales are now
$O({\it\xi})=O(\partial {\it\xi})\sim 1$
and it is easy to show that the advection term scales as
${\sim}{\it\varepsilon}^{2}{\it\Omega}$
. Approaching the high-frequency limit, the penetration depth is much smaller than current bubble radius (
$l\ll 1$
), i.e. when
$\bar{a}^{2}{\it\Omega}\gg 1$
. A series expansion of the advection term taking
${\it\xi}=1+l$
reveals that the magnitude of the dimensionless advection term is independent of
$l$
since:
$$\begin{eqnarray}O\left(\frac{\mathit{Pe}({\it\tau})}{a^{2}}\left(\frac{1}{{\it\xi}^{2}}-{\it\xi}\right)\frac{\partial c}{\partial {\it\xi}}\right)\sim \frac{{\it\varepsilon}{\it\Omega}}{3}(3l)\frac{{\it\varepsilon}}{l}\sim {\it\varepsilon}^{2}{\it\Omega}.\end{eqnarray}$$
We conclude that under the small-amplitude oscillation restriction, the magnitude of the dimensionless advection term in (D 1) is always smaller than the unsteady and diffusive terms by a factor of
${\it\varepsilon}\ll 1$
. Therefore, the dimensionless advection term may be neglected in the oscillatory problem. Moreover, following the discussion presented in § 5.1, we may assume
$a({\it\tau})=\bar{a}({\it\tau})+O({\it\delta})\approx \bar{a}$
to remain constant over an oscillation period and additionally take the liquid to be saturated:
${\it\Upsilon}=1$
.
Hence, letting
$f={\it\xi}c$
, equation (D 1) is reduced to a parabolic equation,
$$\begin{eqnarray}\frac{\partial f}{\partial {\it\tau}}=\frac{1}{\bar{a}^{2}}\frac{\partial ^{2}f}{\partial {\it\xi}^{2}}\end{eqnarray}$$
together with boundary conditions
$$\begin{eqnarray}c(1,{\it\tau})=f(1,{\it\tau})=c_{s}({\it\tau})={\it\varepsilon}\sin {\it\Omega}{\it\tau},\quad c(\infty ,{\it\tau})=f(\infty ,{\it\tau})=0.\end{eqnarray}$$
Resorting to the same treatment given to (A 1) involving Laplace transforms, we arrive from (D 4) at an analogous expression for (2.16), the concentration gradient evaluated at the interface:
$$\begin{eqnarray}-\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}=c_{s}+\frac{c_{s0}\bar{a}}{\sqrt{{\rm\pi}{\it\tau}}}+\bar{a}\int _{0}^{{\it\tau}}\frac{1}{\sqrt{{\rm\pi}({\it\tau}-x)}}\frac{\text{d}c_{s}}{\text{d}x}\,\text{d}x.\end{eqnarray}$$
Under our particular conditions,
$c_{s0}=0$
. The integral term, with
$c_{s}({\it\tau})={\it\varepsilon}\sin {\it\Omega}{\it\tau}$
, may be evaluated from an identity provided by Stepanyants & Yeoh (Reference Stepanyants and Yeoh2009), yielding
$$\begin{eqnarray}\left.-\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}={\it\varepsilon}\sin {\it\Omega}{\it\tau}+\bar{a}{\it\varepsilon}{\it\Omega}\!\int _{0}^{{\it\tau}}\!\frac{\cos {\it\Omega}x}{\sqrt{{\rm\pi}({\it\tau}-x)}}\,\text{d}x={\it\varepsilon}\left[\sin {\it\Omega}{\it\tau}+\bar{a}{\it\Omega}\text{Re}\left\{\frac{\text{e}^{\text{i}{\it\Omega}{\it\tau}}\,\text{erf}\sqrt{\text{i}{\it\Omega}{\it\tau}}}{\sqrt{\text{i}{\it\Omega}}}\right\}\right].\end{eqnarray}$$
The asymptotic expansion as
${\it\tau}\rightarrow \infty$
(since we are interested in the steady periodic state solution) of the last term is
$$\begin{eqnarray}\text{Re}\left\{\frac{\text{e}^{\text{i}{\it\Omega}{\it\tau}}\,\text{erf}\sqrt{\text{i}{\it\Omega}{\it\tau}}}{\sqrt{\text{i}{\it\Omega}}}\right\}\sim \frac{1}{\sqrt{2{\it\Omega}}}\left(\cos {\it\Omega}{\it\tau}+\sin {\it\Omega}{\it\tau}\right).\end{eqnarray}$$
Finally, inserting this result into (D 7) we obtain
$$\begin{eqnarray}\displaystyle \left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}=-{\it\varepsilon}\left[\sin {\it\Omega}{\it\tau}+\bar{a}\sqrt{{\it\Omega}}\sin \!\left({\it\Omega}{\it\tau}+\frac{{\rm\pi}}{4}\right)\right]. & & \displaystyle\end{eqnarray}$$
The last term is hence the contribution from the history integral. As reiterated by Fyrillas & Szeri (Reference Fyrillas and Szeri1994), the problem of the oscillating bubble is completely analogous to Stokes’ second problem of viscous flow near an oscillating flat plate (see for example Landau & Lifshitz Reference Landau and Lifshitz1987). The concentration profile must be a harmonic function in
${\it\tau}$
with the same frequency as the boundary condition
$c_{s}({\it\tau})={\it\varepsilon}\sin {\it\Omega}{\it\tau}$
. Equation (D 4) may be easily solved for harmonic motion considering a solution of the form
$$\begin{eqnarray}f({\it\xi},{\it\tau})={\it\xi}c({\it\xi},{\it\tau})=\text{Im}\{\hat{f}({\it\xi})\text{e}^{\text{i}{\it\Omega}{\it\tau}}\},\end{eqnarray}$$
compatible with the boundary conditions in (D 5). It then follows that
$\hat{f}(1)={\it\varepsilon}$
and
$\hat{f}(\infty )=0$
. The solution for the concentration profile is
$$\begin{eqnarray}c({\it\xi},{\it\tau})=\frac{{\it\varepsilon}}{{\it\xi}}\exp \{-\bar{a}\sqrt{{\it\Omega}/2}({\it\xi}-1)\}\sin \{{\it\Omega}{\it\tau}-\bar{a}\sqrt{{\it\Omega}/2}({\it\xi}-1)\},\end{eqnarray}$$
which renders the following concentration gradient profile:
$$\begin{eqnarray}\frac{\partial c}{\partial {\it\xi}}=-{\it\varepsilon}\exp \{-\bar{a}\sqrt{{\it\Omega}/2}({\it\xi}-1)\}\!\left[\frac{1}{{\it\xi}^{2}}\sin {\it\Omega}{\it\tau}+\frac{\bar{a}\sqrt{{\it\Omega}}}{{\it\xi}}\sin \!\left({\it\Omega}{\it\tau}-\bar{a}\sqrt{{\it\Omega}/2}({\it\xi}-1)+\frac{{\rm\pi}}{4}\right)\right]\!.\end{eqnarray}$$
Evaluating (D 12) on
${\it\xi}=1$
identically results in (D 9).
D.1 Treatment following Fyrillas & Szeri (Reference Fyrillas and Szeri1994)
As demonstrated by Fyrillas & Szeri (Reference Fyrillas and Szeri1994), when considering a large frequency-based Péclet number
${\it\Omega}/(2{\rm\pi})\gg 1$
, to obtain the asymptotic solution for the oscillatory concentration field one must first solve
$$\begin{eqnarray}\frac{\partial c}{\partial \hat{{\it\tau}}}=\frac{\partial ^{2}c}{\partial {\it\eta}^{2}}.\end{eqnarray}$$
The Lagrangian coordinate
${\it\eta}$
, linear time variable
${\it\tau}^{\prime }$
and nonlinear time
$\hat{{\it\tau}}$
are defined as follows:
$$\begin{eqnarray}{\it\eta}=\left(\frac{{\it\Omega}}{2{\rm\pi}}\right)^{1/2}\frac{a^{3}}{3}\left({\it\xi}^{3}-1\right),\quad {\it\tau}^{\prime }=f_{c}t=\frac{{\it\Omega}}{2{\rm\pi}}{\it\tau},\quad \hat{{\it\tau}}({\it\tau}^{\prime })=\int _{0}^{{\it\tau}^{\prime }}a^{4}(x^{\prime })\,\text{d}x^{\prime }.\end{eqnarray}$$
The boundary condition at the interface associated with the oscillatory problem reads
$$\begin{eqnarray}c({\it\eta}=0,\hat{{\it\tau}})=p_{g}({\it\tau}^{\prime })-\langle p_{g}({\it\tau}^{\prime })\rangle _{\hat{{\it\tau}}}.\end{eqnarray}$$
Neglecting inertial and viscous effects,
$p_{g}({\it\tau}^{\prime })$
is defined as
$$\begin{eqnarray}p_{g}({\it\tau}^{\prime })=\frac{p({\it\tau}^{\prime })+{\it\sigma}/a({\it\tau}^{\prime })}{1+{\it\sigma}/\bar{a}},\end{eqnarray}$$
while time averaging is computed according to
$$\begin{eqnarray}\langle f({\it\eta},{\it\tau}^{\prime })\rangle _{\hat{{\it\tau}}}=\frac{1}{\hat{{\it\tau}}(T)}\int _{0}^{T}f({\it\eta},{\it\tau}^{\prime })a^{4}({\it\tau}^{\prime })\,\text{d}{\it\tau}^{\prime },\end{eqnarray}$$
where
$T$
denotes the dimensionless period of oscillation. Expressing this boundary condition as a Fourier series,
$$\begin{eqnarray}c({\it\eta}=0,\hat{{\it\tau}})=\mathop{\sum }_{m=1}^{\infty }\left[a_{m}\cos {\it\omega}_{m}\hat{{\it\tau}}+b_{m}\sin {\it\omega}_{m}\hat{{\it\tau}}\right],\quad \text{with }{\it\omega}_{m}=\frac{2{\rm\pi}m}{\hat{{\it\tau}}(T)},\end{eqnarray}$$
a compatible solution is found to be
$$\begin{eqnarray}c({\it\eta},\hat{{\it\tau}})=\mathop{\sum }_{m=1}^{\infty }\exp \left(-\sqrt{{\it\omega}_{m}/2}{\it\eta}\right)\left[a_{m}\cos \left({\it\omega}_{m}\hat{{\it\tau}}-\sqrt{{\it\omega}_{m}/2}{\it\eta}\right)+b_{m}\sin \!\left({\it\omega}_{m}\hat{{\it\tau}}-\sqrt{{\it\omega}_{m}/2}{\it\eta}\right)\right].\end{eqnarray}$$
The interfacial concentration gradient is thus
$$\begin{eqnarray}\frac{\partial c}{\partial {\it\eta}}({\it\eta}=0,\hat{{\it\tau}})=-\mathop{\sum }_{m=1}^{\infty }{\it\omega}_{m}\left[a_{m}\cos \left({\it\omega}_{m}\hat{{\it\tau}}+{\rm\pi}/4\right)+b_{m}\sin \!\left({\it\omega}_{m}\hat{{\it\tau}}+{\rm\pi}/4\right)\right].\end{eqnarray}$$
Comparing (D 19) and (D 20), it is inferred that the history effect shifts every frequency component of the (negative) interfacial concentration gradient profile with respect to the interfacial concentration by
${\rm\pi}/4$
when
${\it\Omega}/(2{\rm\pi})\gg 1$
.
Next, we shall prove that the solutions (D 19) and (D 20) above will converge with those given in (D 11) and (D 12) given the right set of assumptions. We are considering harmonic pressure forcing:
$p({\it\tau}^{\prime })=1+{\it\varepsilon}\sin 2{\rm\pi}{\it\tau}^{\prime }$
and
$T=1$
. Making the quasi-static radius approximation, i.e. taking
$a({\it\tau}^{\prime })=\bar{a}+O({\it\delta})\approx \bar{a}$
as constant implies that
$\hat{{\it\tau}}=\bar{a}^{4}{\it\tau}^{\prime }$
. Provided the Laplace pressure is small (
${\it\sigma}/\bar{a}\ll 1$
), then
$p_{g}({\it\tau}^{\prime })\approx 1+{\it\varepsilon}\sin 2{\rm\pi}{\it\tau}^{\prime }-{\it\sigma}/\bar{a}$
and the boundary condition (D 15) becomes
$$\begin{eqnarray}c({\it\eta}=0,\hat{{\it\tau}})={\it\varepsilon}\sin 2{\rm\pi}{\it\tau}^{\prime }.\end{eqnarray}$$
From (D 18), we infer that
$a_{1}=0$
,
$b_{1}={\it\varepsilon}$
,
${\it\omega}_{1}=2{\rm\pi}/\bar{a}^{4}$
and
$a_{m}=b_{m}=0$
for
$m>1$
. The solution in (D 19) then simplifies to
$$\begin{eqnarray}c({\it\eta},{\it\tau}^{\prime })={\it\varepsilon}\exp \left(-\frac{\sqrt{{\rm\pi}}}{\bar{a}^{2}}{\it\eta}\right)\sin \!\left(2{\rm\pi}{\it\tau}^{\prime }-\frac{\sqrt{{\rm\pi}}}{\bar{a}^{2}}{\it\eta}\right).\end{eqnarray}$$
Since
${\it\Omega}/(2{\rm\pi})\gg 1$
, we have seen that the boundary layer thickness is very small in comparison with the bubble radius. We are in fact in the limit
${\it\xi}\rightarrow 1^{+}$
. Applying this limit to the Lagrangian coordinate
${\it\eta}$
, we obtain the following identity:
$$\begin{eqnarray}\lim _{{\it\xi}\rightarrow 1^{+}}\{{\it\eta}\}=\lim _{{\it\xi}\rightarrow 1^{+}}\left\{\left(\frac{{\it\Omega}}{2{\rm\pi}}\right)^{1/2}\frac{\bar{a}^{3}}{3}({\it\xi}-1)({\it\xi}^{2}+{\it\xi}+1)\right\}=\left(\frac{{\it\Omega}}{2{\rm\pi}}\right)^{1/2}\bar{a}^{3}({\it\xi}-1).\end{eqnarray}$$
Using this result, we may rewrite the solution for
$c({\it\eta},{\it\tau}^{\prime })$
in (D 22) as a function of our original variables,
$c({\it\xi},{\it\tau})$
. This gives
$$\begin{eqnarray}c({\it\xi},{\it\tau})={\it\varepsilon}\exp \left\{-\bar{a}\sqrt{{\it\Omega}/2}({\it\xi}-1)\right\}\sin \{{\it\Omega}{\it\tau}-\bar{a}\sqrt{{\it\Omega}/2}({\it\xi}-1)\},\end{eqnarray}$$
which in turn yields the following gradient profile:
$$\begin{eqnarray}\frac{\partial c}{\partial {\it\xi}}({\it\xi},{\it\tau})=-{\it\varepsilon}\bar{a}\sqrt{{\it\Omega}}\sin \!\left({\it\Omega}{\it\tau}-\bar{a}\sqrt{{\it\Omega}/2}({\it\xi}-1)+\frac{{\rm\pi}}{4}\right).\end{eqnarray}$$
Applying the limit
${\it\xi}\rightarrow 1^{+}$
alongside
${\it\Omega}\gg 1$
to (D 11) and (D 12), those solutions reduce to expressions (D 24) and (D 25) above. Notice that effect of the bubble curvature is lost (
${\it\xi}$
in the denominator of (5.7) vanishes) consequence of the thin boundary layer approximation. Consequently, the steady-state (first) term contributing to the gradient in (D 9), which dominates for small values of
${\it\Omega}$
, is also lost.
Appendix E. Pressure-corrected radius
We begin by rewriting the mass conservation equation (3.5) without surface tension,
$$\begin{eqnarray}\frac{1}{a}\frac{\text{d}a}{\text{d}{\it\tau}}+\frac{1}{3p}\frac{\text{d}p}{\text{d}{\it\tau}}=\frac{{\it\Lambda}}{a^{2}p}\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}.\end{eqnarray}$$
Integrating this equation in time results in
$$\begin{eqnarray}ap^{1/3}=a_{0}\exp \left\{\int _{0}^{{\it\tau}}\frac{{\it\Lambda}}{a^{2}p}\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}\,\text{d}{\it\tau}^{\prime }\right\}=a_{\mathit{corr}},\end{eqnarray}$$
where
$a_{\mathit{corr}}({\it\tau})$
has been identified using (5.10). Its time derivative is then
$$\begin{eqnarray}\frac{\text{d}a_{\mathit{corr}}}{\text{d}{\it\tau}}=\frac{{\it\Lambda}a_{\mathit{corr}}}{a^{2}p}\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}.\end{eqnarray}$$
Since the prefactor multiplying the interfacial concentration gradient in (E 3) is positive at all times, we can write
$$\begin{eqnarray}\text{sign}\!\left(\frac{\text{d}a_{\mathit{corr}}}{\text{d}{\it\tau}}\right)=\text{sign}\!\left(\left.\frac{\partial c}{\partial {\it\xi}}\right|_{{\it\xi}=1}\right).\end{eqnarray}$$
It follows that the oscillating part of
$a_{\mathit{corr}}$
, namely
$a_{\mathit{corr}}-\bar{a}$
, is in phase with
$$\begin{eqnarray}\int _{0}^{{\it\tau}}\frac{\partial c}{\partial {\it\xi}}({\it\xi}=1,x)\,\text{d}x=-\frac{{\it\varepsilon}}{{\it\Omega}}\left[\sin \!\left({\it\Omega}{\it\tau}-\frac{{\rm\pi}}{2}\right)+\bar{a}\sqrt{{\it\Omega}}\sin \!\left({\it\Omega}{\it\tau}-\frac{{\rm\pi}}{4}\right)\right].\end{eqnarray}$$
Thus, equivalently,
${\it\phi}_{\mathit{corr}}={\it\phi}_{\mathit{grad}}-{\rm\pi}/2$
.
