The history effect in bubble growth and dissolution. Part 1. Theory

TL;DR

This paper develops a theoretical framework to account for the history effect in bubble growth and dissolution, highlighting its influence on diffusion dynamics and rectified diffusion under varying pressure conditions.

Contribution

It introduces a differential equation incorporating a memory integral for bubble dynamics without relying on the quasi-static radius approximation.

Findings

Derived an analytical solution for pressure step changes.

Showed the history effect significantly influences bubble growth.

Analyzed the role of history effect in rectified diffusion.

Abstract

The term `history effect' refers to the contribution of any past mass transfer events between a gas bubble and its liquid surroundings towards the current diffusion-driven growth or dissolution dynamics of that same bubble. The history effect arises from the (non-instantaneous) development of the dissolved gas concentration boundary layer in the liquid in response to changes in the concentration at the bubble interface caused, for instance, by variations of the ambient pressure in time. Essentially, the history effect amounts to the acknowledgement that at any given time the mass flux across the bubble is conditioned by the preceding time-history of the concentration at the bubble boundary. Considering the canonical problem of an isolated spherical bubble at rest, we show that the contribution of the history effect in the current interfacial concentration gradient is fully contained within a memory integral of the interface concentration. Retaining this integral term, we formulate a governing differential equation for the bubble dynamics, analogous to the well-known Epstein-Plesset solution. Our equation does not make use of the quasi-static radius approximation. An analytical solution is presented for the case of multiple step-like jumps in pressure. The nature and relevance of the history effect is then assessed through illustrative examples. Finally, we investigate the role of the history effect in rectified diffusion for a bubble that pulsates under harmonic pressure forcing in the non-inertial, isothermal regime.

Tables

Table 1. Table 1: Coefficients T ~ n subscript ~ 𝑇 𝑛 \tilde{T}_{n} and Δ ​ p n Δ subscript 𝑝 𝑛 \Delta p_{n} of the pressure step function used in the example (see figure 2 ). The resulting coefficients T n subscript 𝑇 𝑛 T_{n} associated to the physical time τ 𝜏 \tau are also included.

$n$	1	2	3	4	5	6
${\tilde{T}}_n$	2	5	6.9	7.4	10.2	12.5
$T_n$	2	7.5	9.8	10.3	15.8	18.6
$\Delta p_n$	-0.1	0.3	-0.2	-0.1	0.3	-0.2

Equations

$$R(t)\sim \text{mbox}{Ja}_{th}\sqrt{D_{th}t},\text{mbox}with\text{mbox}{Ja}_{th}=\frac{{\rho }_l{c}_l\Delta T}{{\rho }_g{h}_{fg}}.$$

$$R(t)\sim \text{mbox}{Ja}_m\sqrt{D_{m}t},\text{mbox}with\text{mbox}{Ja}_m=\frac{M_g\Delta C}{{\rho }_g}.$$

$$R_{corr}(t)=R(t){\left(\frac{P_{\infty }(t)}{P_c}\right)}^{1/3}.$$

$$\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),$$

$$C(R,t)=C_s(t),C(\infty ,t)=C_{\infty },C(r>R_0,0)=C_{\infty },$$

$$C_s=k_H{P}_g,$$

$$P_g=P_{\infty }+\frac{2{\gamma }_{lg}}{R}.$$

$$\frac{4}{3}\text{uppi}{R}^3{P}_g=n{R}_u{T}_{\infty },$$

$$\dot{n}=4\text{uppi}{R}^2{D}_m{\frac{\partial C}{\partial r}|}_{r=R}.$$

$$a=\frac{R}{R_c},p=\frac{P_{\infty }}{P_c},c=\frac{C-C_{\infty }}{k_H{P}_c},c_s=\frac{C_s-C_{\infty }}{k_H{P}_c}.$$

$${\rm Y}=\frac{C_{\infty }}{k_H{P}_c},\Lambda =k_H{R}_u{T}_{\infty },\sigma =\frac{2{\gamma }_{lg}}{R_c{P}_c}.$$

$$c_s=\left(p+\sigma /a\right)-{\rm Y}.$$

$$\tau =\frac{D_m}{R_c^2}t,\mathrm{d}\tilde{\tau }=\frac{D_m}{R^2(t)}\mathrm{d}t,\xi =\frac{r}{R(t)}.$$

$$\mathrm{d}\tau /\mathrm{d}\tilde{\tau }=a^2$$

$$\frac{\partial c}{\partial \tilde{\tau }}+\frac{1}{a}\frac{\mathrm{d}a}{\mathrm{d}\tilde{\tau }}\left(\frac{1}{{\xi }^2}-\xi \right)\frac{\partial c}{\partial \xi }=\frac{1}{{\xi }^2}\frac{\partial }{\partial \xi }\left({\xi }^2\frac{\partial c}{\partial \xi }\right).$$

$$c(1,\tilde{\tau })=c_s(\tilde{\tau })=p+\sigma /a-{\rm Y},c(\infty ,\tilde{\tau })=0,c(\xi >1,0)=0.$$

$$\text{mbox}Pe(\tilde{\tau })=\frac{1}{a}\frac{\mathrm{d}a}{\mathrm{d}\tilde{\tau }}\equiv R\dot{R}/D_m$$

$$\frac{\partial c}{\partial \tilde{\tau }}=\frac{1}{{\xi }^2}\frac{\partial }{\partial \xi }\left({\xi }^2\frac{\partial c}{\partial \xi }\right),$$

$$-{\frac{\partial c}{\partial \xi }|}_{\xi =1}=c_s+\frac{c_{s0}}{\sqrt{\text{uppi}\tilde{\tau }}}+{\int }_0^{\tilde{\tau }}\frac{1}{\text{sqrtsign}\text{uppi}(\tilde{\tau }-\tilde{x})}\frac{\mathrm{d}{c}_s}{\mathrm{d}\tilde{x}}\mathrm{d}\tilde{x}.$$

$$c_{s0}=c_s(0)=p_0+\sigma /a_0-{\rm Y},$$

$$4\text{uppi}{R}^2\dot{R}\left(P_{\infty }+\frac{2{\gamma }_{sl}}{R}\right)+\frac{4\text{uppi}{R}^3}{3}\left({\dot{P}}_{\infty }-\frac{2{\gamma }_{sl}}{R^2}\dot{R}\right)=4\text{uppi}{R}^2\bar{R}{T}_{\infty }D{\frac{\partial C}{\partial r}|}_{r=R}.$$

$$\frac{\mathrm{d}a}{\mathrm{d}\tilde{\tau }}\left(p+\frac{\sigma }{a}\right)+\frac{1}{3}a\left(\frac{\mathrm{d}p}{\mathrm{d}\tilde{\tau }}-\frac{\sigma }{a^2}\frac{\mathrm{d}a}{\mathrm{d}\tilde{\tau }}\right)=\Lambda a{\frac{\partial c}{\partial \xi }|}_{\xi =1}.$$

$$\frac{\mathrm{d}a}{\mathrm{d}\tilde{\tau }}\left(c_s+{\rm Y}\right)+\frac{1}{3}a\frac{\mathrm{d}{c}_s}{\mathrm{d}\tilde{\tau }}=-\Lambda a\left[c_s+\frac{c_{s0}}{\sqrt{\text{uppi}\tilde{\tau }}}+{\int }_0^{\tilde{\tau }}\frac{1}{\text{sqrtsign}\text{uppi}(\tilde{\tau }-\tilde{x})}\frac{\mathrm{d}{c}_s}{\mathrm{d}\tilde{x}}\mathrm{d}\tilde{x}\right].$$

$$\frac{1}{a}\frac{\mathrm{d}a}{\mathrm{d}\tilde{\tau }}+\frac{1}{3p}\frac{\mathrm{d}p}{\mathrm{d}\tilde{\tau }}=-\frac{\Lambda }{p}\left[p-{\rm Y}+\frac{p_0-{\rm Y}}{\sqrt{\text{uppi}\tilde{\tau }}}+{\int }_0^{\tilde{\tau }}\frac{1}{\text{sqrtsign}\text{uppi}(\tilde{\tau }-\tilde{x})}\frac{\mathrm{d}p}{\mathrm{d}\tilde{x}}\mathrm{d}\tilde{x}\right].$$

$$\frac{\mathrm{d}a}{\mathrm{d}\tau }\left(p+\frac{2\sigma }{3a}\right)+\frac{1}{3}a\frac{\mathrm{d}p}{\mathrm{d}\tau }=\frac{\Lambda }{a}{\frac{\partial c}{\partial \xi }|}_{\xi =1}.$$

$$\frac{\mathrm{d}a}{\mathrm{d}\tau }\left(p_0+\frac{2\sigma }{3a}\right)=\frac{\Lambda }{a}{\frac{\partial c}{\partial \xi }|}_{\xi =1}.$$

$$-{\frac{\partial c}{\partial \xi }|}_{\xi =1}=\left[c_s+\frac{c_{s0}}{\sqrt{\text{uppi}\tilde{\tau }}}\right]\approx \left(p_0-{\rm Y}+\frac{\sigma }{a}\right)\left[1+\frac{a}{\sqrt{\text{uppi}\tau }}\right].$$

$$\frac{\mathrm{d}a}{\mathrm{d}\tau }=-\frac{\Lambda (p_0-{\rm Y}+\sigma /a)}{p_0+2\sigma /(3a)}\left[\frac{1}{a}+\frac{1}{\sqrt{\text{uppi}\tau }}\right].$$

$$\frac{\mathrm{d}a}{\mathrm{d}\tilde{\tau }}=-\frac{\Lambda a}{p_0+2\sigma /(3a)}\left[(p_0-{\rm Y})\left(1+\frac{1}{\sqrt{\text{uppi}\tilde{\tau }}}\right)+\frac{\sigma }{a}\left(1+\frac{a/a_0}{\sqrt{\text{uppi}\tilde{\tau }}}\right)\right].$$

$$a(\tilde{\tau })=a_0\exp \left\{-\Lambda (p_0-{\rm Y})\left(\tilde{\tau }+2\text{sqrtsign}\tilde{\tau }/\text{uppi}\right)\right\}.$$

Full text

††1 Fluid Mechanics Group, Universidad Carlos III de Madrid, Avda. de la Universidad 30, 28911 Leganés (Madrid), Spain††2 Physics of Fluids Group, Faculty of Science and Technology, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands

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, 1950), thin boundary layer approximation (Plesset & Zwick, 1954), perturbation techniques (Duda & Vrentas, 1969), infinite series (Tao, 1978), integral methods (Rosner & Epstein, 1972), and self-similar solutions for bubble growth starting from zero initial size (Birkhoff et al., 1958; Scriven, 1959), to cite a few.

All these solutions have in common that they assume

constant ambient pressure during the entire process, and 2. 2.

a uniform concentration (or temperature) field in the liquid as initial conditions.

Perhaps the most widely used ones are the Epstein-Plesset solution, valid for growth and dissolution, based on the quasi-stationary approximation (Epstein & Plesset, 1950), and Scriven’s exact solution for growth that accounts for the advection term in the diffusion equation (Scriven, 1959). The predicted growth rates have been experimentally validated in supersaturated CO2-water solutions by the works of Barker et al. (2002) and Enríquez et al. (2014), among others. Likewise, in the case of dissolution, experimental verification of the Epstein-Plesset equation has been shown for monocomponent (Kapodistrias & Dahl, 2012) and multicomponent (Shim et al., 2014) 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. (2015) 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 (1996) attribute this to bubble growth triggered by rectified diffusion. Houser et al. (2001) 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, 1994). Payvar (1987) 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 (1972) 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. (1969) 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 (1978) 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 (1978) 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, $\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 (1954), Birkhoff et al. (1958) or Scriven (1959) for thermal diffusion growth—under assumptions (i) and (ii)—driven by a temperature difference $\Delta T$ between the bubble boundary and the bulk fluid. The bubble radius scales as

[TABLE]

$D_{\mathit{th}}$ refers to the thermal diffusion coefficient of the liquid; $\mbox{{Ja}}_{\mathit{th}}$ is the Jakob number for heat transfer, where $c_{l}$ is the specific heat of the liquid, $\rho_{l}$ is the liquid density, $\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\mbox{{Ja}}_{\mathit{th}}^{-1}\ll 1$.

Equivalently, for mass-diffusion-controlled growth driven by a (molar) concentration difference $\Delta C$ between the bubble boundary and the bulk fluid, Epstein & Plesset (1950) and Scriven (1959) among others obtained

[TABLE]

$D_{m}$ denotes to the mass diffusion coefficient, $M_{g}$ is the gas molar mass and $\mbox{{Ja}}_{m}$ (Szekely & Martins, 1971) 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 (2015).

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 have been found to be very sensitive to changes in the initial gas concentration profile dissolved in the liquid (Webb et al., 2010).

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\leq 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

[TABLE]

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:

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 evolve around an expression of similar form to Equation (34) in Epstein & Plesset (1950) with an additional memory integral term that accounts for the history effect; 2. 2.

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 section 2, the mathematical formulation is presented and the history integral term is derived. The governing equation for the bubble dynamics is then developed in section 3. Section 4 focuses on a bubble exposed to a train of piece-wise 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, section 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,

[TABLE]

where the dot notation stands for $\mathrm{d}/\mathrm{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

[TABLE]

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,

[TABLE]

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 $\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, 2003) 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 (1987) 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,

[TABLE]

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,

[TABLE]

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

[TABLE]

The governing equations (4)–(9) are best treated in dimensionless form. Consequently, let us define the dimensionless radius, ambient pressure, concentration and interfacial concentration as follows:

[TABLE]

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:

[TABLE]

The parameter $\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}=\Upsilon$. Saturation conditions at $P_{c}$ are described by $\Upsilon=1$, while $\Upsilon<1$ and $\Upsilon>1$ imply undersaturation and supersaturation respectively. Similarly, note that the saturation level at $P_{\infty}$ is given by $\Upsilon/p$. Moreover, $\sigma$ represents the characteristic ratio of the Laplace pressure and $P_{c}$, while $\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

[TABLE]

Finally, we shall introduce two different dimensionless time variables, $\tau$ and $\tilde{\tau}$, in addition to the radial coordinate $\xi$,

[TABLE]

The identity

[TABLE]

directly relates the ‘physical’ dimensionless time $\tau$ and nonlinear time $\tilde{\tau}$. In the upcoming section 2.2, it will be shown that the mathematical nature of the problem encourages to dispose of the ‘physical’ time $\tau$ and work with the nonlinear time $\tilde{\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 $\xi\in[1,\infty)$ scales with the instantaneous bubble radius $R(t)$. Hence, the moving bubble boundary is advantageously always mapped by $\xi=1$. Coordinate $\xi$ is an alternative to the Lagrangian coordinate $\eta=\frac{1}{3}(r^{3}-R^{3}(t))$ which also eliminates the moving boundary problem. The latter transformation is usually employed in the treatment of rectified diffusion (Eller & Flynn, 1965; Fyrillas & Szeri, 1994). However, in the diffusion-dominant regime of interest, the resulting advection-diffusion equation in $\eta$ is much harder, if not impossible, to treat analytically than the analogous equation in $\xi$.

2.2 The history term on the concentration gradient at the bubble surface

The advection-diffusion equation (4) in dimensionless form becomes

[TABLE]

Note that the nonlinear time $\tilde{\tau}$ is the time variable of choice for reasons that will soon become apparent. The boundary and initial conditions ([math]) of the concentration field $c(\xi,\tilde{\tau})$ become

[TABLE]

Only the dimensionless advection term in (15) is explicitly dependent on the bubble dynamics through the prefactor $(\mathrm{d}a/\mathrm{d}\tilde{\tau})/a$. We may characterize this prefactor as a time-dependent Péclet number based on the velocity of the bubble boundary,

[TABLE]

whose magnitude is essentially associated to 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 (15) 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 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)=\Delta c$, $O(\xi)=1$ and denoting $O(\partial\xi)$ as $l$, a characteristic thickness of the boundary layer, defined as $\partial c/\partial\xi|_{\xi=1}=\Delta c/l$. The analysis reveals that the diffusion term magnitude is $\Delta c/l^{2}$, while the advection term has magnitude $O(\mbox{{Pe}})\>\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 $|\mbox{{Pe}}(\tilde{\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\Upsilon/p-1$. Larger ratios are allowed if the gas solubility is very poor ($\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 (1b), through: $|\mbox{{Pe}}|\approx Ja_{m}^{2}=\Lambda|\Upsilon/p-1|$. Thus, a highly soluble gas like CO2 gas in water, ($\Lambda\approx 0.8$ at room temperature) will require that $|\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/\xi^{2}$ is the (dimensionless) physical advection term, associated to the radial component of the velocity field, $u_{r}(r,t)=R^{2}\dot{R}/r^{2}$. The component containing $-\xi$ is the artificially induced advection that compensates for the scaling nature of $\xi$ with $R(t)$. It naturally arises from the nondimensionalization 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 et al., 2015). Nonetheless, the effect of the artificial advection on the interfacial concentration gradient will still be small for $|\mbox{{Pe}}(\tilde{\tau})|\ll 1$.

The concentration field may then be sought by solving

[TABLE]

which has no direct dependency on the bubble dynamics. The choice of $\tilde{\tau}$ over $\tau$ is justified in that $\tilde{\tau}$ allows us to arrive from (4) to (18) based on the single assumption that the advection term is small. Had we chosen to perform an equivalent nondimensionalization with physical time $\tau$ (cf. later Equation (40)), it is only possible to treat the arising $a(\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, 1980) behind the Epstein-Plesset equation resides in these two approximations, namely (1) dropping the advective transport term arising from the interface motion and (2) treating the bubble radius as a constant in the concentration boundary condition at the interface. Approximation (2), 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{\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{\tau}$ allows us to point out an interesting conclusion regarding the appropriateness of neglecting advection effects. Indeed, the advective term in Equation (15) is exactly zero at the bubble interface, as $\left(1/\xi^{2}-\xi\right)=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 (18) for the concentration gradient across the bubble interface (see Appendix A) reads

[TABLE]

The initial interfacial concentration, $c_{s0}$, is given by

[TABLE]

where $p_{0}=p(0)$ and $a_{0}=a(0)$ denote the initial pressure and radius respectively. The first term in the RHS of (19) corresponds to the the steady-state solution determined by the instantaneous interfacial concentration, $c_{s}(\tilde{\tau})$. The second term is the transient component exclusively associated to 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{\tau})$ caused by any prior variations in the ambient pressure $p(\tilde{\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 sections 4 and 5 respectively.

3 The Epstein-Plesset equation with history term

The ideal gas equation of state provided in (8) may be combined with Fick’s law, Equation (9), to obtain the following mass conservation equation:

[TABLE]

In dimensionless form, it reads

[TABLE]

Inserting the expression for the concentration gradient (19) into (22) 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}$:

[TABLE]

The liquid pressure has been related to the surface concentration through $p=c_{s}+\Upsilon-\sigma/a$ through the relation previously provided in (12). In the absence of surface tension, $\sigma=0$, the EPH equation (23) may then be conveniently expressed directly in terms of the ambient pressure $p$,

[TABLE]

To obtain the radius evolution in the physical time, $a(\tau)$, we must numerically integrate the differential EPH equation in (23) or (24) for $a(\tilde{\tau})$ in addition to the differential equation (14) for $\tau(\tilde{\tau})$.

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 (22). Since the original Epstein-Plesset equation is formulated using the linear time $\tau=D_{m}t/R_{c}^{2}$, it is convenient for both time derivatives in $\tilde{\tau}$ to be replaced by derivatives in $\tau$ instead. Recalling that $\mathrm{d}\tau=a^{2}\mathrm{d}\tilde{\tau}$, we may recast (22) as

[TABLE]

The analytical expression for the interfacial concentration gradient is still a function of $\tilde{\tau}$ as given by (19). The Epstein-Plesset equation assumes a constant pressure-history, i.e. $p=p_{0}$ and $\mathrm{d}p/\mathrm{d}\tau=0$. Equation (25) then reduces to

[TABLE]

The interfacial concentration gradient deserves special treatment. The history integral term vanishes, consequence of the constant pressure condition. Moreover, unlike the EPH equations (23) or (24), the Epstein-Plesset equation is based on the quasi-static radius approximation (Weinberg & Subramanian, 1980). The concentration gradient is calculated for a fixed bubble boundary of radius $a$. The solution is then coupled with the equation for mass conservation (25), 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, (14) simplifies to $\tau=a^{2}\tilde{\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 (20), hence $c_{s0}=c_{s}=p_{0}+\sigma/a-\Upsilon$. Under the quasi-static radius approximation, the concentration gradient in (19) then becomes

[TABLE]

Inserting this expression into (26), one finally recovers the Epstein-Plesset equation in dimensionless form (Epstein & Plesset, 1950, Equation (34)):

[TABLE]

An equivalent equation in $\tilde{\tau}$, which does not make use of the quasi-static radius approximation, may be shown to be

[TABLE]

In the absence of surface tension, $\sigma=0$, (29) becomes a separable differential equation in $a$ and $\tilde{\tau}$. Its analytical solution, with $a(0)=a_{0}$, is easily found to be

[TABLE]

The exact analytical solution of $a(\tau)$ to (30) exists in parametric form, with $a(\tilde{\tau})$ and $\tau(\tilde{\tau})$. The physical time $\tau$ is related to $\tilde{\tau}$ through

[TABLE]

which can be derived by solving (14).

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$-th jump taking place at nonlinear time $\tilde{T}_{n}$, the pressure changes by $\Delta p_{n}$. The pressure history and its time derivative may be modelled as

[TABLE]

where $H$ denotes the Heaviside function and $\delta$ is the Dirac delta. As anticipated at the end of section 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 Equation (24) and Equation (18) for the concentration field are derived in Appendix B.

4.1 An example to illustrate the history effect

Let us consider the radius history of a bubble with negligible Laplace pressure ($\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 $\Upsilon=1$. The bubble is then exposed to $N=6$ consecutive jumps in pressure. We prescribe coefficients $\tilde{T}_{n}$ and $\Delta p_{n}$, tabulated in table 1. These are entered into ([math]), 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}<\tau<T^{-}_{2}$ and $T^{+}_{4}<\tau<T^{-}_{5}$) have the same duration in $\tau$ (see figure 2d). Moreover, both expansion stages are identical in the sense that $a(T_{1})=a(T_{4})=a_{0}=1$, $\Delta p_{1}=\Delta p_{4}$ and $p_{1}=p_{4}$. The analytical solution for $a(\tilde{\tau})$ provided by (62) and (63) is plotted in figure 2(b). The solution is more naturally interpreted when presented in $\tau$ (figure 2c). This requires numerical integration of (14) once the analytical solution for $a(\tilde{\tau})$ has been obtained.

At $\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 $\tau=T^{-}_{4}$. Consequently, the memory effect must be behind the differences in bubble growth observed between these two pressure-wise identical expansion periods (figures 2a,c).

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 3b). Both concentration profiles are bounded by an identical value of $c_{s}(\tau^{*})<0$, and the far field concentration $c(\infty,\tau)=0$. The first concentration profile corresponds to a uniform initial condition $c(\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}<\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.

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 $\tau>T^{+}_{6}$, after the expansion-compression cycles have ended (see figure 2). Even though the pressure returns to the saturation value, $c_{s}(\tau\geq T^{+}_{6})=0$, static equilibrium is not instantaneously achieved. Like 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(\tau=T^{+}_{6})\approx 0.9$ to $a(\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 (19) 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 (15), ([math]) and (22) which takes these effects into account. Details of the numerical treatment of the problem may be found in Appendix B.2.

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 $\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}}=\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

neglecting the dimensionless advection term, i.e. setting $\mbox{{Pe}}=0$ in (72); 2. 2.

neglecting surface tension, i.e. setting $\sigma=0$ in (71) and ([math]); 3. 3.

neglecting both the advection term and surface tension (the same scenario as the one analytically solved, i.e. $\mbox{{Pe}}=0$, $\sigma=0$) and 4. 4.

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 $\mbox{{Pe}}=\sigma=0$). History effects are critical in order to properly describe this particular dissolution (compression) stage ($5<\tau<10$). While growth is well determined by existing solutions (Epstein & Plesset, 1950; Scriven, 1959), 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 at 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 $\tau=5^{+}$. This is corroborated in the inset plot of figure 4, where we plot the Epstein-Plesset solution (28), which assumes a uniform initial concentration field (with $p_{0}=1.1$). To enable proper comparison, the initial radius $a_{0}$ at $\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 $\tau=10$, the bubble continues to grow for some time. As discussed in section 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 though, 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 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 (1969), Equation (43). To do so, we subtract the Epstein-Plesset solution (28) with $\sigma=0$ from Equation (43). For our particular set of initial conditions, $p_{0}=a_{0}=1$, we obtain

[TABLE]

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 (28) with $\sigma=0$ from itself. Provided $p_{0}=a_{0}=1$ and $\sigma\ll 1$, this results in

[TABLE]

We may use these formulas to evaluate the characteristic porcentual contribution to the growth rate of each effect with respect to that given by Epstein-Plesset solution (28) with $\sigma=0$. Halfway through the initial growth stage ($\tau=2.5$, $a=1.18$) we determine that advection enhances growth by about 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 phenomenom 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 (1961), Eller & Flynn (1965), Crum & Hansen (1982) or the more recent work by Zhang & Li (2014a). Here, we will only consider low frequency forcing, which has special relevance in seismological events (Tisato et al., 2015), and in the exposure of marine mammals to low-frequency sonars (Crum & Mao, 1996). More specifically, we shall restrict the present analysis to small amplitude, isothermal oscillations in consistency with Equations (7) and (8).

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 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,

[TABLE]

where $\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 $\tau=D_{m}t/R_{c}^{2}$. The dimensionless pressure $p(\tau)$ is then

[TABLE]

The angular forcing frequency $\Omega$ has been non-dimensionlized with the diffusive timescale. The frequency $\Omega/(2\uppi)$ is in fact equivalent to the frequency-based Péclet number typically used in rectified diffusion problems, $\mbox{{Pe}}_{f}=f_{c}R_{c}^{2}/D_{m}$, as introduced by Fyrillas & Szeri (1994). We propose a solution for the bubble dynamics of the form

[TABLE]

where $\bar{a}=\bar{R}/R_{c}$ is the dimensionless equilibrium radius, while $\delta$ is the oscillation amplitude. Lastly, $\phi=\phi(\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 (37) and the analysis soon to follow make use of three underlying assumptions:

1. (a)

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: $\sigma/\bar{a}\ll 1$. 2. (b)

the strain amplitude of the oscillations is small, $\delta/\bar{a}\ll 1$; 3. (c)

the equilibrium radius $\bar{a}(\tau)$ is assumed to vary sufficiently slowly in time to be treated as constant within an individual oscillation period. This timescale for bubble growth/dissolution must be much larger than the period of oscillation, $\Omega^{-1}$.

The derivation and discussion of the range of frequency $\Omega$, pressure amplitude $\varepsilon$ and saturation level $\Upsilon$ for which assumptions (a)–(c) hold is provided in Appendix C. Here, we choose to highlight that assumption (a) requires that the value of $\Omega$ must satisfy

[TABLE]

where $D_{\mathit{th}}$ is the thermal diffusivity of the gas at constant volume and $\Omega_{\mathit{res}}$ is the dimensionless Minnaert resonance frequency. Note that this range may span several orders of magnitude. Moreover, assumption (b) requires that $\varepsilon\ll 1$.

Assumption (a) allows for the gas pressure in the bubble to be the same as the ambient pressure $p(\tau)$ and consequently $c_{s}(\tau)=p(\tau)-\Upsilon$. The radius $a_{p}$ of the equivalent non-soluble bubble mentioned above (for which $\phi_{p}\equiv 0$) then satisfies $p\>a_{p}^{3}=\bar{a}^{3}$. Linearization results in

[TABLE]

Naturally, (a) 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 (a)–(c) essentially simplify this complex mass transfer problem to the point that the role of the history effect can be smoothly rooted out and analyzed 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, 2006), and in high supersaturation levels with potentially large amplitude oscillations as observed in the capillaries of marine mammals (Crum & Mao, 1996; Ilinskii et al., 2008) in which the effect of the viscoelastic medium is undoubtedly important (Zhang & Li, 2014b).

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 (1994). The smooth problem (diffusion timescale $\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 timescale $\tau\sim\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 timescales of the two processes are well separated: $\Omega\gg 1$. Since we are also considering frequencies $\Omega\sim 1$, splitting the problem is no longer possible. Moreover, the smooth solution, whether following the approach of Eller & Flynn (1965) or Fyrillas & Szeri (1994) is constructed from the mass transport equation in Lagrangian spherical coordinates (as opposed to (40)) 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. (2008) 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 to the unsteady boundary condition. In other words, we are just interested in the solution during an individual period of oscillation about $\bar{a}$. Solving the governing mass transport equation directly in the form

[TABLE]

is suitable to this end. Note that $\mbox{{Pe}}(\tau)=(\mathrm{d}a/\mathrm{d}\tau)\>a\equiv R\dot{R}/D_{m}$ is our $\dot{R}$-based Péclet number (as previously defined in Section 2.2). It follows that under the small amplitude oscillation restriction, the magnitude of the dimensionless advection term in (40) is always smaller than the unsteady and diffusive terms by a factor of $O(\varepsilon)$, where $\varepsilon\ll 1$ (see Appendix D). Note that this condition holds true regardless of the magnitude of Pe (cf. $O(\mbox{{Pe}})\sim\varepsilon\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 provide net growth over many cycles.

Moreover, assumption (b) allows to make the approximation $a(\tau)=\bar{a}(\tau)+O(\delta)\approx\bar{a}$, which is constant over an oscillation period according to (c). While this greatly simplifies the problem, it only induces very small errors of $O(\varepsilon)\ll 1$.

Additionally, we will consider saturation conditions $\Upsilon=1$. This is justified in that the bulk concentration of the liquid comes into play in the smooth solution only (Fyrillas & Szeri, 1994). Similarly, surface tension and time averaged bubble interface motion are mostly relevant in the smooth solution.

Implementing these notions , the solution to (40) for the concentration profile is found to be (see Appendix D for the complete derivation)

[TABLE]

which renders the following concentration gradient at the wall:

[TABLE]

The first term in the RHS of (42) refers to the quasi-steady mass flux imposed by the current value of $c_{s}(\tau)=\varepsilon\sin\Omega\tau$. The last term in the RHS 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 $\varepsilon$. The history effect imposes a phase shift in the sinusoidal interfacial concentration gradient, $0<\phi_{\mathit{grad}}(\Omega)\leq\uppi/4$, with respect to the phase of $-c_{s}(\Omega\tau)$ or the oscillatory component of $a_{p}(\tau)$.

As $\Omega\rightarrow 0$, the history integral term vanishes and $\phi_{\mathit{grad}}\rightarrow 0$. The mass flux is in phase with $a_{p}$ (cf. (39)). Conversely, as $\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: $\phi_{\mathit{grad}}=\uppi/4$. In fact, the frequency dependency of this phase shift is independent of $\varepsilon$ and may be shown to be

[TABLE]

Solution (41) may be described as a damped transverse wave with wavenumber $k=\bar{a}\sqrtsign{\Omega/2}$ and phase velocity $\sqrtsign{2\Omega}/\bar{a}$ attenuating over a penetration depth $l\sim 1/(\bar{a}\sqrt{\Omega})$. Figure 5 shows the close agreement of this (advectionless) analytical solution with the full numerical computation (taking into account advection) involving (40) and the mass conservation equation (25) with no surface tension for $\Omega=1$ and $\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 $\varepsilon=0.01$. On another note, the boundary layer thickness $h$ may be best estimated as the wavelength, $h=2\uppi/k$. Over a distance $h$, (horizontal span in figure 5) the concentration amplitude reduces by a factor of $\mathrm{e}^{2\uppi}\approx 540$.

For $\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, $\phi_{\mathit{grad}}=\uppi/4$. The interfacial gradient is flat at $\Omega\tau=3\uppi/4,\>7\uppi/4$ as opposed to $c_{s}(\tau)=0$ taking place at $\Omega\tau=\uppi,\>2\uppi$. For $\Omega=1$, direct evaluation of (43) yields $\phi_{\mathit{grad}}=\uppi/8$ (as later observed in figure 6a).

It stands to reason that solutions (41) and (42) evaluated in the limit $\Omega\gg 1$ converge to those one may obtain following the approach of Fyrillas & Szeri (1994) (inherently valid for $\Omega\gg 1$) subject to the quasi-static radius approximation, i.e. taking $a(\tau)=\bar{a}$ as constant. The complete derivation may be found in Appendix D.1.

5.2 Strain amplitudes and phase

The frequency-dependent $\phi_{\mathit{grad}}$ induced by the history effect implies that pressure-induced contraction of the bubble radius ($p>1$, $c_{s}>0$, $\mathrm{d}a_{p}/\mathrm{d}\tau<0$) is not in phase with dissolution (negative mass diffusion rate, $\partial c/\partial\xi(1,\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 $\Omega$ of the overall strain amplitude and phase lag with respect to $a_{p}$ $(\phi<0)$ of the isothermally pulsating bubble.

Let us then define the dimensionless pressure corrected radius, $a_{\mathit{corr}}(\tau)$, in an analogous manner to $R_{\mathit{corr}}$ in (3) as follows:

[TABLE]

The pressure corrected radius is associated purely with the mass transfer across the interface. It has the useful property that its rate of change, $\mathrm{d}a_{\mathit{corr}}/\mathrm{d}\tau$, always takes the same sign as the interfacial concentration gradient.Consequently, $\phi_{\mathit{corr}}=\phi_{\mathit{grad}}-\uppi/2$, as shown in Appendix E.

The characteristic amplitude of oscillation for the interfacial concentration gradient may be estimated from (42),

[TABLE]

Equation (44) may be rewritten in terms of the different oscillation amplitudes about the equilibrium radius $\bar{a}$,

[TABLE]

Linearization provides $\delta=\delta_{p}+\delta_{\mathit{corr}}+O(\delta_{p}\delta_{\mathit{corr}})$. Amplitudes $\delta_{p}$ and $\delta_{\mathit{corr}}$ may be reasonably approximated as the product of the characteristic time derivative and the characteristic timescale, $1/\Omega$. We obtain

[TABLE]

which is consistent with (39). Likewise, making use of (104),

[TABLE]

Two limiting cases arise. For small $\Omega\ll 1$, we see that $\delta\approx\delta_{\mathit{corr}}\sim\varepsilon\Lambda/\bar{a}\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 $\Omega$ approaches this limit is:

[TABLE]

Taking $\Omega=1$, $\bar{a}=1$ (figure 6), the dominant contribution in $\delta$ and hence $\phi$ still come from $\delta_{\mathit{corr}}$ and $\phi_{\mathit{corr}}$. The relative contributions of $\delta_{\mathit{corr}}$ and $\delta_{p}$ may be estimated from $\delta_{\mathit{corr}}/\delta_{p}\sim 6\Lambda\approx 5$.

For large $\Omega\gg 1$, we see that $\delta\approx\delta_{p}\sim\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:

[TABLE]

Taking $\Omega=100$, $\bar{a}=1$ (figure 7), $a_{p}$ now provides the main contribution since $\delta_{\mathit{corr}}/\delta_{p}\sim 3\Lambda/10\approx 1/4$.

Ultimately, the difference in the phase lags $\phi$ and $\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, 1965), 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 $\mathrm{d}a/\mathrm{d}\tau>0$.

Let us consider the high frequency limit where $\phi=0$, $\phi_{\mathit{corr}}=\uppi/4$ and $a\approx a_{p}$. The area and shell effects act non-symmetrically in the bubble dissolution ($\mathrm{d}a_{\mathit{corr}}/\mathrm{d}\tau<0$) and growth ($\mathrm{d}a_{\mathit{corr}}/\mathrm{d}\tau>0$) periods, which are of equal duration in time. Figure 8 shows that the area and shell effect increase mass transfer during 3/4 of the growth period and 1/4 of the dissolution period. Consequently, these effects reduce mass transfer during 1/4 of the growth period and 3/4 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. 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{\tau}$, which can be related to the standard time once the EPH equation has been solved for $a(\tilde{\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 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 $\uppi/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.

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(\xi,\tilde{\tau})=\xi\>c(\xi,\tilde{\tau})$, and furthermore neglecting any arising $(\mathrm{d}a/\mathrm{d}\tilde{\tau})/a$ terms, Equation (18) transforms to

[TABLE]

Taking the Laplace transform of (51) subject to the initial condition $f(\xi>1,0)=0$ yields

[TABLE]

A solution compatible with boundary conditions $f(1,\tau)=c_{s}(\tilde{\tau})$ and $f(\infty,\tilde{\tau})=0$ is found to be

[TABLE]

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

[TABLE]

The initial value of $c_{s}$ is required. It is given by

[TABLE]

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 $\Upsilon=p_{0}$ and $c_{s0}=\sigma/a_{0}>0$, corresponding to slow dissolution purely driven by the Laplace pressure. Taking then the inverse Laplace transform of (54) results in

[TABLE]

The term containing the Dirac delta $\delta(\tilde{x})$ may be directly integrated. This finally sets the concentration gradient at the interface to be

[TABLE]

It is not surprising that Jones & Zuber (1978) 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

[TABLE]

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, 1987) reads

[TABLE]

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, 1998), 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$-th plateau (constant pressure segment) after jump $n$. Its value is thus

[TABLE]

Time coefficient $\tilde{T}_{n}^{-}$ describes the instant in time right before the $n$-th pressure jump at $\tilde{\tau}=\tilde{T}_{n}$, while $\tilde{T}_{n}^{+}$ refers to the moment right after.

Solution for the Epstein-Plesset with history term equation

The EPH equation (24) may be integrated in time within the $n$-th segment as

[TABLE]

Evaluating this integral finally yields, for segments $n=1$ to $n=N$,

[TABLE]

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 are evidently identical to the solution in (30) and are straightfowardly given by

[TABLE]

Solution for the concentration field

From (53), the Laplace transformed concentration field $\hat{c}(\xi;s)$ may be conveniently split as the product of two separate functions $\hat{f}(s)$ and $\hat{g}(\xi;s)$ as follows:

[TABLE]

The inverse Laplace transforms of these two functions may be shown to be

[TABLE]

Using these results, the concentration field in the time domain $c(\xi,\tilde{\tau})$ may be then computed by means of the convolution theorem. This gives

[TABLE]

Surface tension is to be neglected once again. This renders $c_{s0}=p_{0}-\Upsilon$ and $\mathrm{d}c_{s}/\mathrm{d}\tilde{\tau}=\mathrm{d}p/\mathrm{d}\tilde{\tau}$, the latter modelled in ([math]). After inserting these expressions into (66), the time evolution of the concentration field at every $n$-th segment is finally obtained:

[TABLE]

B.2 Numerical model

The pressure step functions are analytically modelled by the logistic function

[TABLE]

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

[TABLE]

The governing equations will be numerically integrated in nonlinear time $\tilde{\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 $\tau$. To this end, the physical time $\tau$ was computed at each timestep by integrating $\mathrm{d}\tau=a^{2}\mathrm{d}\tilde{\tau}$. Provided $\tilde{k}$ is large enough, (68) and (69) may be computed from

[TABLE]

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 (22), namely

[TABLE]

subject to the prescribed initial size $a_{0}$. The concentration gradient at the interface is solved for numerically from (15) using a finite-differences scheme. Equation (15) is written here again for the reader’s convenience:

[TABLE]

where $\mbox{{Pe}}(\tilde{\tau})=(\mathrm{d}a/\mathrm{d}\tilde{\tau})/a$. The required boundary and initial conditions are

[TABLE]

Appendix C Conditions for small amplitude, isothermal oscillations

Following the work of Prosperetti (1977) (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: $\Omega_{\mathit{th}}\equiv 2\uppi f_{c}R_{c}^{2}/D_{\mathit{th}}<1$, with $D_{\mathit{th}}=k_{g}/\rho_{g}c_{v,g}$. Symbols $k_{g}$, $\rho_{g}$ and $c_{v,g}$ denote the gas thermal conductivity, density and specific heat at constant volume respectively. Consequently, $\Omega$ does not need to be small for this assumption to hold, as long as $\Omega<D_{\mathit{th}}/D_{m}$. As an example, for air at 300 K, $D_{\mathit{th}}/D_{m}\sim 10^{4}$.

Assumption (a) 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 timescale 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, $\Omega$ must additionally satisfy

[TABLE]

where $\Omega_{\mathit{res}}=\sqrt{3P_{c}/\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, 2006). This results in $\Omega\sim 1$–$100$ for a 100 $\upmu$m bubble. Likewise, a 10 $\upmu$m bubble in the tissue of a marine mammal exposed to 5–5000 Hz sonar (Crum & Mao, 1996) corresponds to $\Omega\sim 30$–$3000$.

Except for small values of $\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), $\delta$ will always be of the same order of magnitude as the pressure-induced isothermal expansion and contraction amplitude. Comparison of (37) and (39) gives $\delta/\bar{a}\sim\varepsilon/3$. Therefore, (b) requires that $\varepsilon\ll 1$.

Turning now to (c), consider a bubble in a gas-saturated liquid $\Upsilon=1$, with negligible Laplace pressure $\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 (a)–(c) was found to be reasonably well approximated by Hsieh & Plesset (1961),

[TABLE]

The inverse of the characteristic timescale for growth is then $\tau_{b}^{-1}\sim\Lambda\varepsilon^{2}/\bar{a}$, which must be much smaller than $\Omega$. Let us now consider low forcing frequencies in addition to a highly supersaturated (say $\Upsilon\sim 2$–3) or undersaturated ($\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., 2008), i.e. assumption (c) no longer holds. The inverse of the characteristic time for bubble growth or dissolution may be estimated from (28) as $\tau_{b}^{-1}\sim|\Upsilon-1-\sigma/\bar{a}|\Lambda/\bar{a}^{2}$. Thus, (c) imposes that the solution saturation level, pressure amplitude and frequency must fulfil the following inequalities:

[TABLE]

Finally note that if surface tension and undersaturation are to be overcome, equating Hsieh-Plesset solution (75) with the Epstein-Pesset solution (28) leads to an approximate, simplistic threshold condition for growth (Safar, 1968): $\varepsilon^{2}=(3/2)\left(1-\Upsilon+\sigma/\bar{a}\right)$. As noted by Safar (1968), this threshold is only valid for large isothermal bubbles in a liquid close to saturation ($\Upsilon\approx 1$) under sufficiently high ambient pressures such that the Laplace pressure is comparatively small: $\sigma/\bar{a}\ll 1$. As an example, $\varepsilon\sim 0.01$ would be required to overcome the surface tension driven dissolution of a CO2 bubble of size $R_{c}\sim 100$ $\upmu$m in saturated water at $P_{c}\sim 20$ 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:

[TABLE]

where $\mbox{{Pe}}(\tau)=(\mathrm{d}a/\mathrm{d}\tau)\>a=R\dot{R}/D_{m}$. Its magnitude is given by $O(\mbox{{Pe}})\sim\varepsilon\Omega\bar{a}^{2}/3$ since $O(a)\sim\bar{a}$ and $O(\mathrm{d}a/\mathrm{d}\tau)=O(\mathrm{d}a_{p}/\mathrm{d}\tau)\sim\varepsilon\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 (77) on $\xi=1$ makes the advection term become identically zero. The diffusive term must therefore be of leading order in a layer bounded by $1<\xi<1+l$ where $O(\partial c)=O(\Delta c_{s})\sim\varepsilon$. Parameter $l$ is the dimensionless penetration depth of diffusion using $R(t)$ as the lengthscale, i.e. $l=L/R(t)$. The relevant lengthscales are $O(\xi)\sim 1$ and $O(\partial\xi)\sim l$, while the characteristic timescale per oscillation is $O(\partial\tau)\sim\Omega^{-1}$. The magnitudes of the unsteady and diffusive terms in (77) are given by

[TABLE]

Balancing these two terms yields $l\sim 1/(\bar{a}\sqrt{\Omega})$. Note that the equivalent dimensional penetration depth is $L\sim\sqrtsign{D_{m}/(2\uppi 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}\Omega)^{-1}\sim 1$ is of the order of the diffusion timescale associated to a lengthscale equal to the current bubble size. Lengthscales are now $O(\xi)=O(\partial\xi)\sim 1$ and it is easy to show that the advection term scales as $\sim\varepsilon^{2}\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}\Omega\gg 1$. A series expansion of the advection term taking $\xi=1+l$ reveals that the magnitude of the dimensionless advection term is independent of $l$ since:

[TABLE]

We conclude that under the small amplitude oscillation restriction, the magnitude of the dimensionless advection term in (77) is always smaller than the unsteady and diffusive terms by a factor of $\varepsilon\ll 1$. Therefore, the dimensionless advection term may be neglected in the oscillatory problem. Moreover, following the discussion presented in section 5.1, we may assume $a(\tau)=\bar{a}(\tau)+O(\delta)\approx\bar{a}$ to remain constant over an oscillation period and additionally take the liquid to be saturated: $\Upsilon=1$.

Hence, letting $f=\xi\>c$, Equation (77) is reduced to a parabolic equation,

[TABLE]

together with boundary conditions

[TABLE]

Resorting to the same treatment given to (51) involving Laplace transforms, we arrive from (80) to an analogous expression for (19), the concentration gradient evaluated at the interface:

[TABLE]

Under our particular conditions, $c_{s0}=0$. The integral term, with $c_{s}(\tau)=\varepsilon\sin\Omega\tau$, may be evaluated from an identity provided by Stepanyants & Yeoh (2009), yielding

[TABLE]

The asymptotic expansion as $\tau\rightarrow\infty$ (since we are interested in the steady periodic state solution) of the last term is

[TABLE]

Finally, inserting this result into (83) we obtain

[TABLE]

The last term is hence the contribution from the history integral. As reiterated by Fyrillas & Szeri (1994), 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, 1987). The concentration profile must be a harmonic function in $\tau$ with the same frequency as the boundary condition $c_{s}(\tau)=\varepsilon\sin\Omega\tau$. Equation (80) may be easily solved for harmonic motion considering a solution of the form

[TABLE]

compatible with the boundary conditions in ([math]). It then follows that $\hat{f}(1)=\varepsilon$ and $\hat{f}(\infty)=0$. The solution for the concentration profile is

[TABLE]

which renders the following concentration gradient profile:

[TABLE]

Evaluating (88) on $\xi=1$ identically results in (85).

D.1 Treatment following Fyrillas & Szeri (1994)

As demonstrated by Fyrillas & Szeri (1994), when considering a large frequency-based Péclet number $\Omega/(2\uppi)\gg 1$, to obtain the asymptotic solution for the oscillatory concentration field one must first solve

[TABLE]

The Lagrangian coordinate $\eta$, linear time variable $\tau^{\prime}$ and nonlinear time $\hat{\tau}$ are defined as follows:

[TABLE]

The boundary condition at interface associated to the oscillatory problem reads

[TABLE]

Neglecting inertial and viscous effects, $p_{g}(\tau^{\prime})$ is defined as

[TABLE]

while time averaging is computed according to

[TABLE]

where $T$ denotes the dimensionless period of oscillation. Expressing this boundary condition as a Fourier series,

[TABLE]

a compatible solution is found to be

[TABLE]

The interfacial concentration gradient is thus

[TABLE]

Comparing (95) and (96), 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 $\uppi/4$ when $\Omega/(2\uppi)\gg 1$.

Next, we shall prove that the solutions (95) and (96) above will converge with those given in (87) and (88) given the right set of assumptions. We are considering harmonic pressure forcing: $p(\tau^{\prime})=1+\varepsilon\sin 2\uppi\tau^{\prime}$ and $T=1$. Making the quasi-static radius approximation, i.e. taking $a(\tau^{\prime})=\bar{a}+O(\delta)\approx\bar{a}$ as constant implies that $\hat{\tau}=\bar{a}^{4}\tau^{\prime}$. Provided the Laplace pressure is small ($\sigma/\bar{a}\ll 1$), then $p_{g}(\tau^{\prime})\approx 1+\varepsilon\sin 2\uppi\tau^{\prime}-\sigma/\bar{a}$ and the boundary condition (91) becomes

[TABLE]

From (94), we infer that $a_{1}=0$, $b_{1}=\varepsilon$, $\omega_{1}=2\uppi/\bar{a}^{4}$, and $a_{m}=b_{m}=0$ for $m>1$. The solution in (95) then simplifies to

[TABLE]

Since $\Omega/(2\uppi)\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 $\xi\rightarrow 1^{+}$. Applying this limit to the Langrangian coordinate $\eta$, we obtain the following identity:

[TABLE]

Using this result, we may rewrite the solution for $c(\eta,\tau^{\prime})$ in (98) as a function of our original variables, $c(\xi,\tau)$. This gives

[TABLE]

which in turn yields the following gradient profile:

[TABLE]

Applying the limit $\xi\rightarrow 1^{+}$ alongside $\Omega\gg 1$ to (87) and (88), those solutions reduce to expressions (100) and (101) above. Notice that effect of the bubble curvature is lost ($\xi$ in the denominator of (41) vanishes), consequence of the thin boundary layer approximation. Consequently, the steady-state (first) term contributing to the gradient in (85), which dominates for small values of $\Omega$, is also lost.

Appendix E Pressure corrected radius

We begin by rewriting the mass conservation equation (25) without surface tension,

[TABLE]

Integrating this equation in time results in

[TABLE]

where $a_{\mathit{corr}}(\tau)$ has been identified using (44). Its time derivative is then

[TABLE]

Since the prefactor multiplying the interfacial concentration gradient in (104) is positive at all times, we can write

[TABLE]

It follows that the oscillating part of $a_{\mathit{corr}}$, namely $a_{\mathit{corr}}-\bar{a}$, is in phase with

[TABLE]

Thus, equivalently, $\phi_{\mathit{corr}}=\phi_{\mathit{grad}}-\uppi/2$.