2.1 Experimental set-up and procedure

Although the experimental set-up has been described in a previous work (Enríquez et al. Reference Enríquez, Hummelink, Bruggert, Lohse, Prosperetti, van der Meer and Sun2013), a brief description is included here for convenience. The facility is fed with water that is demineralised in a purifier (MilliQ A10) and degassed by making it flow through a filter (MiniModule, Liquicel, Membrana). This water enters into the mixing chamber (see figure 2), that has been previously flushed with $\text{CO}_{2}$ to purge the air from the system. There the water is stirred in the presence of $\text{CO}_{2}$ , kept at the desired saturation pressure for approximately 45 min. Finally, the experimental tank is pressurised with $\text{CO}_{2}$ at this same pressure and then slowly flooded with the carbonated water, so bubbles do not appear during the filling. This preparation procedure ensures that in the experimental tank there are no other gas species present within the liquid or gas phases apart from $\text{CO}_{2}$ (at least in quantifiable amounts).

Figure 2.

Sketch of the experimental set-up. See Enríquez et al. (Reference Enríquez, Hummelink, Bruggert, Lohse, Prosperetti, van der Meer and Sun2013) for a detailed description.

Placed at the centre of the experimental tank there is a silicon chip, treated to become hydrophilic, with a black-silicon hydrophobic pit ( $50~\unicode[STIX]{x03BC}\text{m}$ in radius) at its centre. The role of this pit is to force a single bubble to nucleate at a fixed location in a repeatable way. Furthermore, in order to avoid slight temperature variations to affect the diffusion-driven bubble dynamics, the measurement tank is kept at a constant temperature by means of an external chiller.

Once the measurement tank is filled with the carbonated water, the following experimental procedure is followed:

1. (i) The pressure is lowered below the saturation value until a bubble nucleates at the pit and grows up to the desired size, $R_{i}$ .

2. (ii) The tank pressure is set again to the saturation value. Then, the pressure is finely adjusted manually until the bubble size does not vary in an observable way for approximately five minutes. The pressure at which this occurs will be hereafter the one used in the calculations as the saturation pressure. Notice that this procedure allows us to determine the saturation pressure with more accuracy than that given by the pressure controller during the mixing process.

3. (iii) At time $t_{1}$ (see figure 3), the pressure is lowered by a given amount, $\unicode[STIX]{x0394}p_{1}$ , during a prescribed time $T=t_{2}-t_{1}$ . This turns the liquid supersaturated, which leads to bubble growth.

4. (iv) Subsequently, at time $t_{2}$ , the pressure is increased by an amount $\unicode[STIX]{x0394}p_{2}$ , which causes undersaturation and the bubble to shrink.

5. (v) When the bubble becomes slightly smaller than $R_{i}$ at time $t_{3a}$ , the pressure is gradually set back to the saturation level $P_{0}$ ( $t_{3b}$ ) by means of a pressure drop $\unicode[STIX]{x0394}p_{3}$ . During this short period ( $t_{3a}$ – $t_{3b}$ ) the bubble expands and grows.

6. (vi) During a short time $T_{s}$ after $t_{3b}$ , the pressure remains at saturation but the bubble keeps growing and attains the initial size $R_{i}$ at $t_{4}$ due to the history effect (portrayed in figure 1). At this point ( $t_{4}$ ), growth step (iii) and subsequent dissolution step (iv) are immediately repeated at identical $\unicode[STIX]{x0394}p_{1}$ and $\unicode[STIX]{x0394}p_{2}$ respectively.

Figure 3.

Experimental procedure during which the bubble is exposed to two identical supersaturation–undersaturation cycles. The lower plot shows the pressure that the bubble is subjected to, whilst the upper one illustrates its radius time evolution. In brief, once a bubble is stabilised at a given radius $R_{i}$ , it is forced to a growth cycle for a prescribed time $T$ and then to a dissolution cycle down to a size slightly smaller than $R_{i}$ ( $t=t_{3a}$ ), such that a short time $T_{s}$ after the pressure returns to the initial level $P_{0}$ ( $t=t_{3b}$ ), the combination of previous gas expansion (during $t_{3a}$ – $t_{3b}$ ) and history lead the bubble size to the initial radius ( $t=t_{4}$ ). An identical pressure cycle is immediately imposed, which results in a different time evolution of the bubble radius due to the history effect.

It is imperative to realise that the pressure and bubble size conditions at $t_{4}$ , just before the pressure jump, are identical to the initial conditions at $t_{1}$ . Namely, $R(t_{4})=R(t_{1})$ , $P_{\infty }(t_{4})=P_{\infty }(t_{1})$ and $\text{d}P_{\infty }(t_{4})/\text{d}t=\text{d}P_{\infty }(t_{1})/\text{d}t=0$ . The zero-pressure time derivative condition is extremely important to ensure that the bubble is not under the effect of any previous pressure-induced volumetric expansion or compression at the time when $\unicode[STIX]{x0394}p_{1}$ is suddenly imposed. Note that the sole purpose of the pressure change between $t_{3a}$ and $t_{3b}$ and subsequent stabilisation period ( $t_{3b}$ and $t_{4}$ ) is precisely to enforce this last condition. This complex procedure allows us to directly isolate and quantify the history effect through direct comparison. Any differences between the first growth rate (during $t_{1}$ – $t_{2}$ ) and second growth rate ( $t_{4}$ – $t_{5}$ ) must be purely attributed to the history effect. The differences arise because at $t_{1}$ the bubble is in equilibrium with its surroundings (uniform concentration field) and the contribution of history term is essentially negligible. At $t_{4}$ , however, the concentration field surrounding the bubble has evolved. It is no longer uniform, and the bubble is no longer in equilibrium: thus, the contribution of the history term is now larger.

2.2 Experimental results and discussion

In this subsection we present the results of four experiments that manifest the effects of history in bubble growth and shrinkage. Three of the experiments were carried out as described above, while in the fourth the order of the growth and shrinkage stages was swapped, i.e. first shrinking and then growing.

In the experiments where the bubble is first made to grow (figures 4–6), the most apparent difference between the two growth stages is the somewhat larger bubble size achieved at the end of the second stage (see panels labelled $(a)$ ). This is a consequence of a more important effect, namely the higher growth rate found during the first instants of the second stage, as predicted by the modified Epstein–Plesset equation with history effects provided in Part 1. To illustrate this point, panels $(c)$ show the time derivative of the ambient bubble radius. In all cases the growth rate during the second stage lies above that of the first one, although both curves eventually converge at longer times, when the memory of the previous dissolution stage damps out. As demonstrated in Part 1, the $\text{CO}_{2}$ accumulated around the bubble during a dissolution stage yields a steeper concentration gradient at the interface that, in turn, leads to a faster growth rate at short times once the pressure is reduced and the liquid is supersaturated again. As the growth progresses, the influence of the initial concentration profile becomes weaker and both growth rates converge to the same curve.

Figure 4.

Results for experiment 1 showing the time-histories of (a) the measured bubble radius $R$ in response to (b) the imposed pressure $P_{\infty }(t)$ . In (c), the rate of growth of the ambient radius $R_{0}$ , defined in (1.1), is plotted for the two growth cycles. The time axis is initialised at $t_{1}$ or $t_{4}$ accordingly.

Figure 5.

Results for experiment 2 (see caption of figure 4). The range of pressures is slightly different to the ones exposed in figure 4. However, the history effect is repeatable.

Figure 6.

Results for experiment 3 (see caption of figure 4).

Figure 7.

Results for experiment 4 (see caption of figure 4). This time, the direction of the pressure jumps is inverted thereby replacing the two growth stages observed in experiments 1–3 with two dissolution stages.

As demonstrated in Part 1, the $\text{CO}_{2}$ accumulated around the bubble during a dissolution stage yields a steeper concentration gradient at the interface that, in turn, leads to a faster growth rate after a short transient time once the pressure is dropped and the liquid is supersaturated again. As the growth progresses, the influence of the initial concentration profile becomes weaker and both growth rates converge to the same curve. During the very early times after the pressure drop (up to approximately ten seconds later), the contribution of the history effect on mass transfer is masked by the large growth rates induced by the sudden decrease of the interfacial concentration (induced by this pressure drop via Henry’s law) that leads to a steep interfacial concentration gradient. The change in growth rates between the first and second cycles is experimentally indiscernible. This does however agree with theory, as one may observe from the numerically computed rates in figure 3 of Part 1 (Peñas-López et al. Reference Peñas-López, Parrales, Rodríguez-Rodríguez and van der Meer2016b ).

It is interesting to compare this behaviour with that found when the bubble is forced to first dissolve and then to grow (figure 7). Although unavoidable experimental limitations of the control of the pressure in the facility in this case result in a somewhat noisier time derivative of the ambient radius, the same qualitative behaviour is found. Namely, the magnitude of the rate of change of the radius is larger in the second dissolution stage, thus leading to a smaller radius at the end of this stage. Analogously to what occurred in experiments 1–3, this is a consequence of the local depletion of $\text{CO}_{2}$ near the bubble caused by the intermediate growth stage.

Besides illustrating the history effect in the growth and dissolution of bubbles, these experiments will serve as benchmark cases for the numerical simulations described in the following sections. These numerical analyses will allow us to quantify the relative importance of the different physical effects that play a role in the processes illustrated in figures 4–7 which, besides diffusion, include surface tension, boundary-induced advection and natural convection.

In the theory that follows, we will assume that the bubble remains strictly spherical at all times. Two experimental snapshots depicting the upper and lower extremes in bubble size are provided in figure 8. The bubble is actually attached to a cylindrical pit of $50~\unicode[STIX]{x03BC}\text{m}$ diameter and $30~\unicode[STIX]{x03BC}\text{m}$ depth. The gas volume contained inside the pit can be neglected compared to the total volume of the gas bubble. In the experiments in which the cycles start with a growth phase, where $R>200~\unicode[STIX]{x03BC}\text{m}$ , the bubble remains spherical, as observed in figure 8(a). Only at the smallest radii during the experiments starting with a dissolution phase, we observe a spherical cap, figure 8(b). However, the assumption of perfectly spherical bubble at all time yields a relative error of less than 3 % as compared to the actual gas volume of the spherical cap and the pit. Therefore, the assumption of strictly spherical gas bubble for the analysis is more than justified.

Figure 8.

Bubble snapshots at both extremes of the bubble size range measured during our experiments. The largest radius is (a) $R=358~\unicode[STIX]{x03BC}\text{m}$ , corresponding to the maximum radius attained during experiment 1 (see later figure 4), whereas (b) $R=92~\unicode[STIX]{x03BC}\text{m}$ is the smallest radius, obtained during the dissolution experiment 4 (see figure 7). The radius is computed by means of the light-blue circumference fitted to the bubble contour. The horizontal red line marks the height of the bubble–substrate contact line, below which there is the reflection of the bubble on the substrate surface.

3.1 Mass transfer problem

The transport of dissolved gas species in the liquid is governed by the following mass transport equation, usually referred to as the advection–diffusion equation,

(3.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}t}+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}C=D_{m}\hspace{2.22198pt}\unicode[STIX]{x1D6FB}^{2}C,\end{eqnarray}$$

where $C(\boldsymbol{x},t)$ is the molar concentration field, $\boldsymbol{U}(\boldsymbol{x},t)$ is the velocity vector field and $D_{m}$ is the mass diffusion coefficient. The initial concentration of dissolved gas is assumed to be uniform throughout the liquid and equal to $C_{\infty }$ , equal to the gas concentration in the far field. The boundary condition of zero-mass flux holds across the impermeable wall. The concentration boundary condition at the bubble surface, $C_{i}(t)$ , is given by Henry’s law,

(3.2) $$\begin{eqnarray}C_{i}(t)=k_{H}P_{g}(t),\end{eqnarray}$$

where $k_{H}$ is Henry’s (molar-based) solubility constant and $P_{g}(t)$ is the total gas pressure inside the bubble. A constant temperature environment $T_{\infty }$ is assumed, i.e. $k_{H}$ remains constant, while the ambient pressure $P_{\infty }(t)$ is set to vary with time $t$ . The bubble gas volume is related to the gas content and pressure via the equation of state for an ideal gas,

(3.3) $$\begin{eqnarray}{\textstyle \frac{4}{3}}\unicode[STIX]{x03C0}R^{3}P_{g}=nR_{u}T_{\infty },\end{eqnarray}$$

where $n(t)$ is the number of gas moles inside the bubble and $R_{u}$ denotes the universal gas constant. The total gas pressure inside the bubble, $P_{g}$ , considering liquid–gas surface tension $\unicode[STIX]{x1D6FE}_{lg}$ , but neglecting inertial and viscous effects inside the gas phase, is given by

(3.4) $$\begin{eqnarray}P_{g}=P_{\infty }+2\unicode[STIX]{x1D6FE}_{lg}/R.\end{eqnarray}$$

The mass transfer problem is closed with Fick’s first law, which sets the molar flow rate of gas across the bubble surface $S$ to be

(3.5) $$\begin{eqnarray}{\dot{n}}=D\int _{S}\unicode[STIX]{x1D735}C\boldsymbol{\cdot }\hat{\boldsymbol{n}}\hspace{2.22198pt}\text{d}S,\end{eqnarray}$$

where $\text{d}S$ is an infinitesimal area element of the bubble surface, and $\hat{\boldsymbol{n}}$ is the outward-pointing unit normal from the bubble surface.

Equations (3.1)–(3.5) represent the mass transfer problem equations. These must be coupled with the equations of motion from which the velocity field $\boldsymbol{U}(\boldsymbol{x},t)$ may be computed.

3.2 Modelling density-driven natural convection

The dissolved gas concentration profile around the bubble implies a non-uniform density field of the surrounding liquid–gas solution which may trigger the onset of ‘density-induced natural convection’ (Takemura, Liu & Yabe Reference Takemura, Liu and Yabe1996). The change in solution density may be quantified through the concentration expansion coefficient, $\unicode[STIX]{x1D706}$ , usually defined as (Bataller et al. Reference Bataller, Miqueu, Plantier, Daridon, Jaber, Abbasi, Saghir and Bou-Ali2009)

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D706}=\frac{1}{\unicode[STIX]{x1D70C}_{l}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}C},\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}_{l}$ is the density of the pure solvent. Any change in the solution density is therefore assumed proportional to the change in dissolved gas concentration. For dilute, monosolute solutions, the concentration expansion coefficient is approximately given by (see appendix A)

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D706}\approx \frac{M_{g}}{\unicode[STIX]{x1D70C}_{l}}-\overline{V}_{g}^{\infty },\end{eqnarray}$$

where $M_{g}$ is the gas molar mass and $\overline{V}_{g}^{\infty }$ is the (temperature dependent) partial molar volume of the solute in the solvent at infinite dilution. For $\text{CO}_{2}$ gas in pure water, $\overline{V}_{g}^{\infty }\approx 34.2~\text{cm}^{-3}~\text{mol}^{-1}$ (Harvey, Kaplan & Burnett Reference Harvey, Kaplan and Burnett2005), which results in $\unicode[STIX]{x1D706}\approx 9.8~\text{cm}^{-3}~\text{mol}^{-1}$ .

The variations in density considered here are small, of the order of 0.1 %. However, these variations are sufficiently large to have a non-negligible effect on the motion of the flow. Consequently, it was deemed appropriate to take this effect into account via the Boussinesq approximation. This essentially results in the inclusion of a non-uniform buoyancy term imposed by the local dissolved gas concentration into the Navier–Stokes equation (3.11). The Boussinesq approximation allows for the flow to be regarded as incompressible when treating the continuity equation. Therefore, the incompressibility condition,

(3.8) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{U}=0,\end{eqnarray}$$

always holds. Moreover, the density $\unicode[STIX]{x1D70C}(\boldsymbol{x},t)$ of the liquid–gas solution shall be approximated as constant, $\unicode[STIX]{x1D70C}(\boldsymbol{x},t)=\unicode[STIX]{x1D70C}_{l}$ , in both the inertial and viscous terms of the momentum equation. However, we shall allow small variations in density in the body force (gravity) term. The density field $\unicode[STIX]{x1D70C}(\boldsymbol{x},t)$ may then be expressed as

(3.9) $$\begin{eqnarray}\unicode[STIX]{x1D70C}(\boldsymbol{x},t)=\unicode[STIX]{x1D70C}_{l}+\unicode[STIX]{x1D70C}^{\ast }(\boldsymbol{x},t),\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}^{\ast }(\boldsymbol{x},t)$ is the density perturbation field arising from the non-uniform concentration field and, evidently, $|\unicode[STIX]{x1D70C}^{\ast }|\ll \unicode[STIX]{x1D70C}_{l}$ . Similarly, the pressure field in the solution may be decomposed into

(3.10) $$\begin{eqnarray}P(\boldsymbol{x},t)=P_{\infty }(t)+P_{h}(\boldsymbol{x})+P^{\ast }(\boldsymbol{x},t).\end{eqnarray}$$

Here $P_{h}$ is the background hydrostatic pressure, $\unicode[STIX]{x1D735}P_{h}=\unicode[STIX]{x1D70C}_{l}\boldsymbol{g}$ , where $\boldsymbol{g}$ denotes the gravitational acceleration, and $P^{\ast }(\boldsymbol{x},t)$ is the pressure perturbation arising from the fluid motion. It likewise follows that for our experiments, $P_{h}\ll P_{\infty }$ and $|P^{\ast }|\ll P_{\infty }$ .

3.3 Equations of motion in terms of the streamfunction and vorticity

Making use of (3.9) and (3.10), the Navier–Stokes momentum equation may be written as

(3.11) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{U}}{\unicode[STIX]{x2202}t}+(\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{U}=-\frac{\unicode[STIX]{x1D735}P^{\ast }}{\unicode[STIX]{x1D70C}_{l}}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{U}+\frac{\unicode[STIX]{x1D70C}^{\ast }\boldsymbol{g}}{\unicode[STIX]{x1D70C}_{l}},\end{eqnarray}$$

where $\unicode[STIX]{x1D708}$ is the kinematic viscosity of the liquid. Since the flow is axisymmetric around the vertical ( $z$ ) axis, if we are able to define an orthogonal set of coordinates $\unicode[STIX]{x1D702},\unicode[STIX]{x1D709},\unicode[STIX]{x1D719}$ where $\unicode[STIX]{x1D719}$ is the angle of rotation around the vertical axis, then the velocity field has only two components, $\boldsymbol{U}=U_{\unicode[STIX]{x1D702}}(\unicode[STIX]{x1D702},\unicode[STIX]{x1D709})\hspace{2.22198pt}\hat{\boldsymbol{e}}_{\boldsymbol{\unicode[STIX]{x1D702}}}+U_{\unicode[STIX]{x1D709}}(\unicode[STIX]{x1D702},\unicode[STIX]{x1D709})\hspace{2.22198pt}\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D743}}$ , and the whole problem may be treated as two-dimensional. The vorticity field $\unicode[STIX]{x1D734}$ is then also unidirectional and a vorticity scalar, $\unicode[STIX]{x1D6FA}$ , exists:

(3.12) $$\begin{eqnarray}\unicode[STIX]{x1D734}=\unicode[STIX]{x1D6FA}\hspace{2.22198pt}\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D753}}=\unicode[STIX]{x1D735}\times \boldsymbol{U}.\end{eqnarray}$$

Taking the curl of (3.11) eliminates the pressure term and the vorticity scalar transport equation is obtained:

(3.13) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}t}+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D6FA}\hspace{2.22198pt}\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D753}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{U}-\unicode[STIX]{x1D708}\unicode[STIX]{x1D735}\times \left(\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D6FA}\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D753}})\right)\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D753}}+\frac{1}{\unicode[STIX]{x1D70C}_{l}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{\ast }\times \boldsymbol{g})\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D753}}.\end{eqnarray}$$

It follows from (3.6) that $\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D706}\unicode[STIX]{x1D70C}_{l}\unicode[STIX]{x1D735}C$ , so the vorticity transport equation becomes

(3.14) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}t}+h_{\unicode[STIX]{x1D719}}\,\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\unicode[STIX]{x1D6FA}/h_{\unicode[STIX]{x1D719}}\right)={\mathcal{L}}^{2}\unicode[STIX]{x1D6FA}+\unicode[STIX]{x1D706}(\unicode[STIX]{x1D735}C\times \boldsymbol{g})\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D753}}.\end{eqnarray}$$

Here, we have made use of the following linear operator:

(3.15) $$\begin{eqnarray}{\mathcal{L}}^{2}=\frac{1}{h^{2}}\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\left(\frac{1}{h_{\unicode[STIX]{x1D719}}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\left(h_{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D6FA}\right)\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\left(\frac{1}{h_{\unicode[STIX]{x1D719}}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\left(h_{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D6FA}\right)\right)\right],\end{eqnarray}$$

where $h_{\unicode[STIX]{x1D719}}$ denotes the scale factor in the $\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D753}}$ direction and $h$ the scale factor in both the $\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D743}}$ and $\hat{\boldsymbol{e}}_{\boldsymbol{\unicode[STIX]{x1D702}}}$ directions. The coordinate system and scale factors will be introduced quantitatively in § 4. Incompressibility, along with the axisymmetric nature of the flow, allow for the velocity field to be expressed in terms of a scalar streamfunction, $\unicode[STIX]{x1D6F9}$ :

(3.16) $$\begin{eqnarray}\boldsymbol{U}=\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D6F9}/h_{\unicode[STIX]{x1D719}}\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D753}}).\end{eqnarray}$$

Combining (3.12) and (3.16) results in the following equation for the streamfunction,

(3.17) $$\begin{eqnarray}{\mathcal{L}}^{2}\left(\unicode[STIX]{x1D6F9}/h_{\unicode[STIX]{x1D719}}\right)=-\unicode[STIX]{x1D6FA}.\end{eqnarray}$$

In § 4, it will be shown how the fluid motion may be obtained by simultaneously solving $\unicode[STIX]{x1D6F9}$ from (3.17) and $\unicode[STIX]{x1D6FA}$ from (3.14) numerically by employing a streamfunction–vorticity method in dynamic tangent-sphere coordinates. It will be seen that the boundary conditions for $\unicode[STIX]{x1D6F9}$ and $\unicode[STIX]{x1D6FA}$ can be determined from those for $\boldsymbol{U}$ through careful analysis. From the physical point of view, the velocity field must satisfy the kinematic and dynamic (zero-shear stress) boundary conditions along the moving bubble boundary, in addition to the no-slip condition at the wall.

3.4 On the parameters and time scales of the problem

This subsection intends to shed light on the physics governing the diffusion-driven growth and dissolution of a bubble attached to a flat plate. More specifically, the goal is to prove that the concentration and velocity fields evolve over very disparate time scales, which will allow for an efficient procedure to numerically solve the problem formulated in previous subsections.

The processes involved in this problem introduce four characteristic time scales: $t_{s}$ for bubble growth and dissolution, $t_{m}$ for mass diffusion of the dissolved gas, $t_{v}$ for viscous diffusion of momentum and $t_{b}$ for the density-induced convection. Let $U$ denote the characteristic flow velocity. When the advection induced by the moving boundary dominates over natural convection, then $U$ is the interface velocity $U_{s}\sim {\dot{R}}$ . When convection overcomes boundary-induced advection, then $U$ becomes the convection velocity $U_{b}$ . The characteristic length scale is the bubble radius $R$ . For mass-diffusion-controlled growth driven by a molar concentration difference $\unicode[STIX]{x0394}C$ between the bubble boundary and the bulk fluid, the flow behaviour may be characterised using three dimensionless parameters. These are the Jakob (Szekely & Martins Reference Szekely and Martins1971) and Grashof numbers for mass transfer, in addition to the Schmidt number, defined as follows:

(3.18a-c ) $$\begin{eqnarray}Ja=\frac{M_{g}|\unicode[STIX]{x0394}C|}{\unicode[STIX]{x1D70C}_{g}},\quad Gr=\frac{\unicode[STIX]{x1D706}|\unicode[STIX]{x0394}C|gR^{3}}{\unicode[STIX]{x1D708}^{2}},\quad Sc=\frac{\unicode[STIX]{x1D708}}{D_{m}},\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}_{g}$ is the density of the gas bubble and $g$ is the magnitude of the acceleration due to gravity. The Jakob number, Ja, may be regarded as a measure of the driving force for bubble growth induced by the concentration difference and gas solubility. The Grashof number, Gr, represents the ratio of buoyancy (convection) and viscous forces. The Schmidt number, Sc, is the ratio of momentum and mass diffusivities.

Here we shall consider bubble growth or dissolution that is primarily driven by mass diffusion. We may then use the approximate result obtained by Epstein & Plesset (Reference Epstein and Plesset1950) or Scriven (Reference Scriven1959) to estimate bubble growth as

(3.19) $$\begin{eqnarray}R\sim Ja\sqrt{D_{m}t}.\end{eqnarray}$$

It then follows that the bubble growth time scale and boundary-induced advection velocity scales are

(3.20a,b ) $$\begin{eqnarray}U_{s}={\dot{R}}\sim \frac{Ja^{2}D_{m}}{R},\quad t_{s}=\frac{R}{U_{s}}\sim \frac{R^{2}}{Ja^{2}D_{m}}.\end{eqnarray}$$

The magnitudes of the terms in the mass transport equation (3.1) are

(3.21a-c ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}t}\sim \frac{|\unicode[STIX]{x0394}C|}{t_{m}},\quad \boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}C\sim \frac{U|\unicode[STIX]{x0394}C|}{R},\quad D\unicode[STIX]{x1D6FB}^{2}C\sim \frac{D_{m}|\unicode[STIX]{x0394}C|}{R^{2}},\end{eqnarray}$$

where $t_{m}$ is the characteristic time required for a significant concentration change over characteristic length scale $R$ . Similarly, taking $\unicode[STIX]{x1D6FA}\sim U/R$ , the magnitudes of the terms in the vorticity transport equation (3.14) are

(3.22a-d ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}t}\sim \frac{U}{R\hspace{2.22198pt}t_{v}},\quad \boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FA}\sim \frac{U^{2}}{R^{2}},\quad \unicode[STIX]{x1D708}{\mathcal{L}}^{2}\unicode[STIX]{x1D6FA}\sim \frac{\unicode[STIX]{x1D708}U}{R^{3}},\quad \unicode[STIX]{x1D706}\unicode[STIX]{x1D735}C\times \boldsymbol{g}\sim \frac{\unicode[STIX]{x1D706}|\unicode[STIX]{x0394}C|g}{R},\end{eqnarray}$$

where similarly $t_{v}$ refers to the time required for a significant vorticity change over the same characteristic length scale $R$ . For $Sc\sim 1$ or $Sc\gg 1$ , the characteristic convection velocity and time scale may be obtained from a balance between the viscous term (3.22c ) and the buoyancy term (3.22d ) in the vorticity transport equation,

(3.23a,b ) $$\begin{eqnarray}U_{b}\sim \frac{\unicode[STIX]{x1D706}|\unicode[STIX]{x0394}C|gR^{2}}{\unicode[STIX]{x1D708}},\quad t_{b}=\frac{R}{U_{b}}\sim \frac{\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D706}|\unicode[STIX]{x0394}C|gR}.\end{eqnarray}$$

The ratio of velocities is given by $U_{b}/U_{s}=Gr\hspace{2.22198pt}Sc/Ja^{2}$ . The ratio of the advection term (3.21b ) and the diffusive term (3.21c ) in the mass transport equation yields a Péclet number, $Pe=UR/D$ . The ratio of the advection term (3.22b ) over the diffusive term (3.22c ) in the vorticity transport equation similarly yields a Reynolds number, $Re=UR/\unicode[STIX]{x1D708}$ . Neglecting natural convection, setting $U=U_{s}$ gives $Pe=Ja^{2}$ and $Re=Ja^{2}/Sc$ . Likewise, natural convection dominating over boundary-induced advection, $U=U_{b}$ , results in $Pe=Gr\hspace{2.22198pt}Sc$ and $Re=Gr$ .

From the above analysis, we may conclude that mass and momentum diffusion will clearly dominate over advection and natural convection provided $Gr\hspace{2.22198pt}Sc<1$ and $Ja<1$ (i.e. $Pe$ and $Re$ are small). In such a case, the mass diffusion and viscous time scales are the leading time scales in the mass transport and vorticity transport equations respectively. The unsteady term in each transport equation may then be balanced by the corresponding diffusive term, yielding

(3.24a,b ) $$\begin{eqnarray}t_{m}\sim \frac{R^{2}}{D_{m}},\quad t_{v}\sim \frac{R^{2}}{\unicode[STIX]{x1D708}}.\end{eqnarray}$$

The ratios between the mass diffusion time scale and the other time scales are

(3.25a-c ) $$\begin{eqnarray}\frac{t_{m}}{t_{s}}=Ja^{2},\quad \frac{t_{m}}{t_{b}}=Gr\hspace{2.22198pt}Sc,\quad \frac{t_{m}}{t_{v}}=Sc.\end{eqnarray}$$

As reference for the conditions explored in this work, a $\text{CO}_{2}$ gas bubble with $R=0.25$  mm growing in a 15 % supersaturated $\text{CO}_{2}$ -water solution at 5 bar and 293 K, with $\unicode[STIX]{x1D706}=9.8~\text{cm}^{-3}~\text{mol}^{-1}$ , results in $Ja=0.12$ , $Gr=0.038$ and $Sc=523$ . The Rayleigh number is $Gr\hspace{2.22198pt}Sc=19.6$ . Under these conditions, intentionally similar to those of our experiments, equation (3.25) translates to

(3.26) $$\begin{eqnarray}t_{s}>t_{m}\sim t_{b}\gg t_{v}.\end{eqnarray}$$

The vorticity/velocity field around a bubble evolves at a time scale $t_{v}$ provided by the viscous diffusion of momentum. This time scale is much faster than the time scale $t_{m}$ of mass transfer, i.e. the time required to observe a significant change in the concentration field surrounding the bubble. Likewise, the time scale $t_{s}$ in which a substantial change in the bubble radius may be observed is significantly larger than $t_{m}$ . This means that the thin boundary layer approximation (valid when $t_{s}\ll t_{m}$ , i.e. $Ja\gg 1$ ), while suitable for treating the fast growth of bubbles in highly supersaturated liquids (Rosner & Epstein Reference Rosner and Epstein1972), is clearly not applicable here.

The time scale of interest is of course $t_{m}$ . Let us neglect density-driven convection for the moment. At every time step of this slow time scale $t_{m}$ , provided $t_{m}\gg t_{v}$ ( $Sc\gg 1$ ), viscous action ensures that the flow always reaches (over a much faster time scale $t_{v}$ ) a steady-state solution. In other words, at every time step of $t_{m}$ , the advection term in the mass transport equation may then be computed from the steady-state vorticity (hence velocity) solution imposed by the instantaneous concentration field and interface velocity. We shall refer to it as the quasi-steady advection approximation. It is worth pointing out that bubbles of other gases with solubility parameter, $\unicode[STIX]{x1D6EC}$ , smaller than that of $\text{CO}_{2}$ (such as nitrogen or oxygen) can be described as well with this approximation, as $t_{m}$ will be much smaller than $t_{s}$ . In these cases, the history effect – which is a diffusive effect – will be even more apparent since boundary-driven advection will have a smaller influence.

Considering now density-driven convection, provided $t_{b}\gg t_{v}$ ( $Gr\ll 1$ ), then viscosity is able to establish a quasi-steady velocity field in a time much shorter than that taken by buoyancy to induce changes in the flow. In other words, although buoyancy must be taken into account to properly compute the velocity field around the bubble, it does not affect the validity of the quasi-steady advection approximation. This knowledge will now be used in the next section when implementing the equations into a numerical model.

4.1 Non-dimensionalisation

We begin by introducing the dimensionless time, radius and Cartesian coordinates, ambient pressure and mole number:

(4.1a-e ) $$\begin{eqnarray}\unicode[STIX]{x1D70F}=\frac{D_{m}}{R_{i}^{2}}t,\quad a=\frac{R}{R_{i}},\quad \tilde{\boldsymbol{x}}=\frac{\boldsymbol{x}}{R_{i}},\quad p=\frac{P_{\infty }}{P_{0}},\quad \unicode[STIX]{x1D707}=\frac{R_{u}T_{\infty }}{4/3\unicode[STIX]{x03C0}R_{i}^{3}P_{0}}n.\end{eqnarray}$$

In this work we have chosen the characteristic radius $R_{i}$ to be the initial radius $R(t=0)$ . Similarly, the characteristic ambient pressure $P_{0}$ corresponds to the initial liquid pressure, $P_{\infty }(0)$ , whereas the mole number $n$ is made dimensionless with that contained in a bubble of radius $R_{i}$ , immersed in a liquid at pressure $P_{0}$ and in the absence of surface tension. Note that the time scale of choice has been that of mass diffusion, $t_{m}$ , presented in (3.24a ). Additionally, the molar concentration field $C$ and the interfacial molar concentration $C_{i}$ may be non-dimensionalised through

(4.2a,b ) $$\begin{eqnarray}c=\frac{C-C_{\infty }}{k_{H}P_{0}},\quad c_{i}=\frac{C_{i}-C_{\infty }}{k_{H}P_{0}}.\end{eqnarray}$$

The dimensionless counterparts of the vorticity scalar, velocity and streamfunction are

(4.3a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D714}=\frac{R_{i}^{2}}{D_{m}}\unicode[STIX]{x1D6FA},\quad \boldsymbol{u}=\frac{R_{i}}{D_{m}}\boldsymbol{U},\quad \unicode[STIX]{x1D713}=\frac{1}{R_{i}D_{m}}\unicode[STIX]{x1D6F9}.\end{eqnarray}$$

Lastly, let us present the following dimensionless parameters and dimensionless numbers:

(4.4a-d ) $$\begin{eqnarray}\unicode[STIX]{x1D6F6}=\frac{C_{\infty }}{k_{H}P_{0}},\quad \unicode[STIX]{x1D6EC}=k_{H}R_{u}T_{\infty },\quad \unicode[STIX]{x1D70E}=\frac{2\unicode[STIX]{x1D6FE}_{lg}}{R_{i}P_{0}},\quad Gr_{0}=\frac{\unicode[STIX]{x1D706}k_{H}P_{0}gR_{i}^{3}}{\unicode[STIX]{x1D708}^{2}}.\end{eqnarray}$$

The parameter $\unicode[STIX]{x1D6F6}$ refers to the initial saturation level of the solution, $\unicode[STIX]{x1D6EC}$ is a solubility parameter, $\unicode[STIX]{x1D70E}$ is the surface tension parameter, while $Gr_{0}$ is a reference Grashof number, $Gr_{0}=Gr\,k_{H}P_{0}/|\unicode[STIX]{x0394}C|$ .

4.2 The tangent-sphere coordinate system

The problem can be conveniently recast in dimensionless tangent-sphere spatial coordinates $(\unicode[STIX]{x1D702},\unicode[STIX]{x1D709},\unicode[STIX]{x1D719})$ , where

(4.5a-c ) $$\begin{eqnarray}\tilde{x}=2a\frac{\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D702}^{2}+\unicode[STIX]{x1D709}^{2}}\cos \unicode[STIX]{x1D719},\quad {\tilde{y}}=2a\frac{\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D702}^{2}+\unicode[STIX]{x1D709}^{2}}\sin \unicode[STIX]{x1D719},\quad \tilde{z}=2a\frac{\unicode[STIX]{x1D709}}{\unicode[STIX]{x1D702}^{2}+\unicode[STIX]{x1D709}^{2}}.\end{eqnarray}$$

The contours of $\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D709}$ satisfy the following inverse relations (Moon & Spencer Reference Moon and Spencer1988; Batchelor Reference Batchelor1979):

(4.6a,b ) $$\begin{eqnarray}\tilde{x}^{2}+{\tilde{y}}^{2}+\tilde{z}^{2}=(2a/\unicode[STIX]{x1D702})\sqrt{\tilde{x}^{2}+{\tilde{y}}^{2}},\quad \tilde{x}^{2}+{\tilde{y}}^{2}+\tilde{z}^{2}=\left(2a/\unicode[STIX]{x1D709}\right)\tilde{z}.\end{eqnarray}$$

These coordinates, represented in figure 10, scale with the dimensionless radius of the bubble, $a(\unicode[STIX]{x1D70F})$ . The scale factors are defined as

(4.7a,b ) $$\begin{eqnarray}\tilde{h}=\frac{h_{\unicode[STIX]{x1D702}}}{R_{i}}=\frac{h_{\unicode[STIX]{x1D709}}}{R_{i}}=\frac{2a}{\unicode[STIX]{x1D702}^{2}+\unicode[STIX]{x1D709}^{2}},\quad \displaystyle \tilde{h}_{\unicode[STIX]{x1D719}}=\frac{h_{\unicode[STIX]{x1D719}}}{R_{i}}=\frac{2a\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D702}^{2}+\unicode[STIX]{x1D709}^{2}}.\end{eqnarray}$$

Figure 10.

Contour lines of the tangent-sphere $\unicode[STIX]{x1D702}$ , $\unicode[STIX]{x1D709}$ coordinates, plotted in the $y=0$ ( $\unicode[STIX]{x1D719}=0$ ) Cartesian plane. $\unicode[STIX]{x1D702}=0$ lies on the $z$ -axis,  $\unicode[STIX]{x1D702}\rightarrow \infty$ at the contact point. The horizontal wall lies on the $\unicode[STIX]{x1D709}=0$ isosurface, while the bubble surface is always mapped by $\unicode[STIX]{x1D709}=1$ . The separation of the plotted contours is uniform ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D702}=\unicode[STIX]{x0394}\unicode[STIX]{x1D709}=0.1$ ).

The partial time derivative of any scalar function $f$ described by fixed Cartesian coordinates $(x,y,z)$ expands as the material derivative when described by our $R(t)$ -scaling spatial coordinates ( $\unicode[STIX]{x1D702},\unicode[STIX]{x1D709}$ ). Taking the partial derivative of $f$ with respect to time $\unicode[STIX]{x1D70F}$ , we find

(4.8) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}f(\tilde{x},{\tilde{y}},\tilde{z},\unicode[STIX]{x1D70F})=\frac{\text{D}}{\text{D}\unicode[STIX]{x1D70F}}f(\unicode[STIX]{x1D702}(\unicode[STIX]{x1D70F}),\unicode[STIX]{x1D709}(\unicode[STIX]{x1D70F}),\unicode[STIX]{x1D70F})=\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}+\unicode[STIX]{x1D702}^{\prime }\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\unicode[STIX]{x1D709}^{\prime }\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}.\end{eqnarray}$$

The prime notation ( $^{\prime }$ ) denotes $\text{d}/\text{d}\unicode[STIX]{x1D70F}$ . The terms containing $\unicode[STIX]{x1D702}^{\prime }$ and $\unicode[STIX]{x1D709}^{\prime }$ represent the apparent advection of a quiescent fluid relative to our time-varying coordinate system. Let us define the a priori unknown corresponding (dimensionless) apparent velocity field as $\boldsymbol{u}_{rel}(\unicode[STIX]{x1D702},\unicode[STIX]{x1D709},\unicode[STIX]{x1D70F})=u_{rel,\unicode[STIX]{x1D702}}\hspace{2.22198pt}\hat{\boldsymbol{e}}_{\boldsymbol{\unicode[STIX]{x1D702}}}+u_{rel,\unicode[STIX]{x1D709}}\hspace{2.22198pt}\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D743}}$ . The advection term on $f$ would then be

(4.9) $$\begin{eqnarray}(\boldsymbol{u}_{rel}\boldsymbol{\cdot }\tilde{\unicode[STIX]{x1D735}})f=\frac{u_{rel,\unicode[STIX]{x1D702}}}{\tilde{h}}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\frac{u_{rel,\unicode[STIX]{x1D709}}}{\tilde{h}}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}},\end{eqnarray}$$

where the operator $\tilde{\unicode[STIX]{x1D735}}=R_{i}\unicode[STIX]{x1D735}$ is dimensionless. Comparing (4.8) and (4.9) immediately reveals that $u_{rel,\unicode[STIX]{x1D702}}=\tilde{h}\unicode[STIX]{x1D702}^{\prime }$ and $u_{rel,\unicode[STIX]{x1D709}}=\tilde{h}\unicode[STIX]{x1D709}^{\prime }$ . Thus, the dimensionless velocity field of our scaling coordinate system (relative to any fixed point in the physical domain) is just equal to $-\boldsymbol{u}_{rel}=-(\tilde{h}\unicode[STIX]{x1D702}^{\prime }\hspace{2.22198pt}\hat{\boldsymbol{e}}_{\boldsymbol{\unicode[STIX]{x1D702}}}+\tilde{h}\unicode[STIX]{x1D709}^{\prime }\hspace{2.22198pt}\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D743}})$ . Differentiating (4.6a ) and (4.6b ) independently with respect to $\unicode[STIX]{x1D70F}$ , one finds that

(4.10a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D702}^{\prime }=a^{\prime }\unicode[STIX]{x1D702}/a,\quad \unicode[STIX]{x1D709}^{\prime }=a^{\prime }\unicode[STIX]{x1D709}/a.\end{eqnarray}$$

4.3 Streamfunction–vorticity formulation

The streamfunction satisfies the following equation,

(4.11) $$\begin{eqnarray}\tilde{{\mathcal{L}}}^{2}(\unicode[STIX]{x1D713}/\tilde{h}_{\unicode[STIX]{x1D719}})=-\unicode[STIX]{x1D714},\end{eqnarray}$$

which makes use of the dimensionless operator $\tilde{{\mathcal{L}}}^{2}=R_{0}^{2}{\mathcal{L}}^{2}$ , and (4.11) may be rewritten in terms of the coordinates  $\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D709}$ as

(4.12) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D702}^{2}-\unicode[STIX]{x1D709}^{2}}{2a\tilde{h}_{\unicode[STIX]{x1D719}}^{2}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\frac{\unicode[STIX]{x1D702}\unicode[STIX]{x1D709}}{a\tilde{h}_{\unicode[STIX]{x1D719}}^{2}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}+\frac{1}{\tilde{h}^{2}\tilde{h}_{\unicode[STIX]{x1D719}}}\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}^{2}}\right)=-\unicode[STIX]{x1D714}.\end{eqnarray}$$

Boundary conditions are derived and explained in appendix B. These are

(4.13a-d ) $$\begin{eqnarray}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702},0)=0,\quad \unicode[STIX]{x1D713}(\unicode[STIX]{x1D702},1)=-\frac{2a^{2}a^{\prime }}{\unicode[STIX]{x1D702}^{2}+1},\quad \unicode[STIX]{x1D713}(0,\unicode[STIX]{x1D709})=-2a^{2}a^{\prime },\quad \unicode[STIX]{x1D713}(\infty ,\unicode[STIX]{x1D709})=0.\end{eqnarray}$$

Once $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702},\unicode[STIX]{x1D709})$ is known, from the definition in (4.12), the dimensionless velocity components may be computed by numerically differentiating $\unicode[STIX]{x1D713}$ according to

(4.14a,b ) $$\begin{eqnarray}u_{\unicode[STIX]{x1D702}}=\frac{1}{\tilde{h}_{\unicode[STIX]{x1D719}}\tilde{h}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}},\quad u_{\unicode[STIX]{x1D709}}=-\frac{1}{\tilde{h}_{\unicode[STIX]{x1D719}}\tilde{h}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}.\end{eqnarray}$$

Expressions for the velocity field in the Cartesian reference frame are included in appendix C. The dimensionless vorticity transport equation reads

(4.15) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}+\tilde{h}_{\unicode[STIX]{x1D719}}\,\boldsymbol{u}\boldsymbol{\cdot }\tilde{\unicode[STIX]{x1D735}}\left(\frac{\unicode[STIX]{x1D714}}{\tilde{h}_{\unicode[STIX]{x1D719}}}\right)=Sc\hspace{2.22198pt}\tilde{{\mathcal{L}}}^{2}\unicode[STIX]{x1D714}+Gr_{0}\hspace{2.22198pt}Sc^{2}\hspace{2.22198pt}(\tilde{\unicode[STIX]{x1D735}}c\times \hat{\boldsymbol{g}})\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D753}},\end{eqnarray}$$

where $\hat{\boldsymbol{g}}=\boldsymbol{g}/g$ . The Rayleigh number $Gr_{0}\hspace{2.22198pt}Sc$ represents the ratio between the reference buoyant transport of momentum and diffusive transport of mass. In ( $\unicode[STIX]{x1D702}$ , $\unicode[STIX]{x1D709}$ ) coordinates, the vorticity transport equation becomes

(4.16a ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}=L\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+M\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}+P\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}^{2}}\right)+Q+S\unicode[STIX]{x1D714},\end{eqnarray}$$

where $L$ , $M$ , $P$ , $Q$ and $S$ are time- and space-dependent coefficients given by

(4.16b,c ) $$\begin{eqnarray}L=-Sc\hspace{2.22198pt}\frac{\unicode[STIX]{x1D702}^{2}-\unicode[STIX]{x1D709}^{2}}{2a\tilde{h}_{\unicode[STIX]{x1D719}}}-\frac{a^{\prime }\unicode[STIX]{x1D702}}{a}-\frac{u_{\unicode[STIX]{x1D702}}}{\tilde{h}},\quad M=-Sc\hspace{2.22198pt}\frac{\unicode[STIX]{x1D709}\unicode[STIX]{x1D702}}{a\tilde{h}_{\unicode[STIX]{x1D719}}}-\frac{a^{\prime }\unicode[STIX]{x1D709}}{a}-\frac{u_{\unicode[STIX]{x1D709}}}{\tilde{h}},\end{eqnarray}$$

(4.16d,e ) $$\begin{eqnarray}P=\frac{Sc}{\tilde{h}^{2}},\quad Q=-Gr_{0}\hspace{2.22198pt}Sc^{2}\hspace{2.22198pt}\left(\frac{\unicode[STIX]{x1D702}^{2}-\unicode[STIX]{x1D709}^{2}}{2a}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\frac{\unicode[STIX]{x1D702}\unicode[STIX]{x1D709}}{a}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\right),\end{eqnarray}$$

(4.16f ) $$\begin{eqnarray}S=-\frac{Sc}{\tilde{h}}\,\frac{\unicode[STIX]{x1D702}^{2}+2\unicode[STIX]{x1D709}^{2}+\unicode[STIX]{x1D709}^{4}/\unicode[STIX]{x1D702}^{2}}{\left(\unicode[STIX]{x1D702}^{2}+\unicode[STIX]{x1D709}^{2}\right)^{2}}-\frac{1}{\tilde{h}\unicode[STIX]{x1D702}}\,\frac{u_{\unicode[STIX]{x1D702}}\left(\unicode[STIX]{x1D702}^{2}-\unicode[STIX]{x1D709}^{2}\right)+2u_{\unicode[STIX]{x1D709}}\,\unicode[STIX]{x1D702}\unicode[STIX]{x1D709}}{\unicode[STIX]{x1D702}^{2}+\unicode[STIX]{x1D709}^{2}}.\end{eqnarray}$$

Coefficient $Q$ contains the density-driven convection term. The cross product was evaluated by first expressing $\hat{\boldsymbol{g}}$ in its ( $\hat{\boldsymbol{e}}_{\boldsymbol{\unicode[STIX]{x1D702}}}$ , $\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D743}}$ ) components via the transformation matrix described in appendix C. Boundary conditions are derived and explained in appendix B.3. These are

(4.17a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}(\unicode[STIX]{x1D702},0)=\frac{\tilde{h}u_{n}}{Sc\hspace{2.22198pt}\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}_{v}},\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}(\unicode[STIX]{x1D702},1)=\frac{4a^{\prime }\unicode[STIX]{x1D702}}{a(\unicode[STIX]{x1D702}^{2}+1)}+\frac{2}{a}u_{\unicode[STIX]{x1D702}}(\unicode[STIX]{x1D702},1),\end{eqnarray}$$

(4.17c,d ) $$\begin{eqnarray}\unicode[STIX]{x1D714}(0,\unicode[STIX]{x1D709})=0,\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}(\infty ,\unicode[STIX]{x1D709})=0.\end{eqnarray}$$

Here, $\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}_{v}$ is the computational time step for the viscous transport of momentum. We shall report its meaning in § 4.5.

4.4 Formulation for the mass transfer problem

In dimensionless form, the advection–diffusion equation (3.1) becomes

(4.18) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}+\boldsymbol{u}\boldsymbol{\cdot }\tilde{\unicode[STIX]{x1D735}}c=\tilde{\unicode[STIX]{x1D6FB}}^{2}c,\end{eqnarray}$$

or equivalently,

(4.19a ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}=F\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+G\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}+H\left(\frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}+\frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}^{2}}\right),\end{eqnarray}$$

with

(4.19b-d ) $$\begin{eqnarray}F=-\frac{\unicode[STIX]{x1D702}^{2}-\unicode[STIX]{x1D709}^{2}}{2a\tilde{h}\unicode[STIX]{x1D702}}-\frac{a^{\prime }\unicode[STIX]{x1D702}}{a}-\frac{u_{\unicode[STIX]{x1D702}}}{\tilde{h}},\quad G=-\frac{\unicode[STIX]{x1D709}}{a\tilde{h}}-\frac{a^{\prime }\unicode[STIX]{x1D709}}{a}-\frac{u_{\unicode[STIX]{x1D709}}}{\tilde{h}},\quad H=\frac{1}{\tilde{h}^{2}}.\end{eqnarray}$$

The velocity field components $u_{\unicode[STIX]{x1D702}}(\unicode[STIX]{x1D702},\unicode[STIX]{x1D709},\unicode[STIX]{x1D70F})$ and $u_{\unicode[STIX]{x1D709}}(\unicode[STIX]{x1D702},\unicode[STIX]{x1D709},\unicode[STIX]{x1D70F})$ must of course be computed beforehand as detailed in § 4.3. Boundary conditions for the concentration are

(4.20a-d ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}(\unicode[STIX]{x1D702},0)=0,\quad c(\unicode[STIX]{x1D702},1)=c_{i},\quad \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}(0,\unicode[STIX]{x1D709})=0,\quad \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}(\infty ,\unicode[STIX]{x1D709})=0.\end{eqnarray}$$

The interfacial concentration $c_{i}(\unicode[STIX]{x1D70F})$ , appealing to Henry’s law, is given by

(4.21) $$\begin{eqnarray}c_{i}=\left(p+\frac{\unicode[STIX]{x1D70E}}{a}\right)-\unicode[STIX]{x1D6F6}.\end{eqnarray}$$

The dimensionless molar flow rate $\unicode[STIX]{x1D707}^{\prime }$ across the interface may be computed non-dimensionalising (3.5), resulting in

(4.22) $$\begin{eqnarray}\unicode[STIX]{x1D707}^{\prime }=-3\unicode[STIX]{x1D6EC}a\int _{0}^{\infty }\frac{\unicode[STIX]{x1D702}}{1+\unicode[STIX]{x1D702}^{2}}\left.\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\right|_{\unicode[STIX]{x1D709}=1}\text{d}\unicode[STIX]{x1D702}.\end{eqnarray}$$

Finally, the last equation remaining is the dimensionless ideal gas law,

(4.23) $$\begin{eqnarray}\left(p+\frac{\unicode[STIX]{x1D70E}}{a}\right)a^{3}=\unicode[STIX]{x1D707}.\end{eqnarray}$$

4.5 Numerical procedure

The experiments were simulated by numerically solving the governing equations presented in the previous section. To do so, we used a second-order finite-difference discretisation in space and an implicit Euler method in time. The latter was chosen in search for unconditional stability, bearing in mind that the grid spacing becomes infinitesimally small as $\unicode[STIX]{x1D702}\rightarrow \infty$ .

We have seen in § 3.4 that this problem involves multiple scales. The governing equations have been non-dimensionalised in time with the time scale $t_{m}=R_{i}^{2}/D_{m}$ for mass diffusion, $\unicode[STIX]{x1D70F}=t/t_{m}$ . The mass transport equation will therefore require a time step $\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}_{m}$ that suitably advances within $t_{m}$ . However, the momentum transport equation requires a much smaller time scale $t_{v}$ and consequently, $\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}_{v}\sim \unicode[STIX]{x0394}\unicode[STIX]{x1D70F}_{m}/Sc$ . Stability of the scheme requires to advance in $\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}_{v}$ . Doing so, however, the overall number of time iterations and computational cost required to span the whole duration of the experiments would be exceedingly high.

This issue could be solved making use of the quasi-steady advection approximation, validated in § 3.4. Essentially, this approximation allows to advance the simulation in $\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}_{m}$ rather than in $\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}_{v}$ . After time step $\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}_{m}$ , time advances from $\unicode[STIX]{x1D70F}_{n}$ to $\unicode[STIX]{x1D70F}_{n+1}$ , where subscript $n$ refers to the $n$ th time iteration. Given the actual concentration field $c_{n}$ , suppose mass transfer across the bubble interface results in a change of radius from $a_{n}$ to $a_{n+1}$ with corresponding rate $a^{\prime }$ . What is the vorticity field $\unicode[STIX]{x1D714}_{n+1}$ for this new configuration? We first make the initial guess: $\unicode[STIX]{x1D714}_{n+1}=\unicode[STIX]{x1D714}_{n}$ . The vorticity field is then allowed to independently evolve through the vorticity transport equation, advancing with time step $\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}_{v}$ . By virtue of the quasi-steady advection approximation, the concentration field, bubble radius and velocity are treated as invariants. After $k$ iterations, the vorticity converges to the steady-state solution. We used the following criterion for convergence:

(4.24) $$\begin{eqnarray}\unicode[STIX]{x1D700}_{max}^{2}=\text{max}\left\{\frac{(\unicode[STIX]{x1D714}_{k+1}-\unicode[STIX]{x1D714}_{k})^{2}}{\text{max}(\unicode[STIX]{x1D714}_{k}^{2})}\right\}<10^{-6}.\end{eqnarray}$$

The resulting vorticity and velocities are then used to solve the advection–diffusion equation to find $c_{n+1}$ . The process is then repeated for the following time step. This way, the overall number of iterations and computational time are greatly reduced.

The numerical procedure followed at every time iteration $n$ consists of the following steps:

1. (i) Update time $\unicode[STIX]{x1D70F}_{n}$ and primary variables: vorticity $\unicode[STIX]{x1D714}_{n}$ , concentration $c_{n}$ and mole number $\unicode[STIX]{x1D707}_{n}$ . Define the time step $\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}_{m}$ .

2. (ii) Obtain value of the externally applied pressures $p_{n}$ and $p_{n+1}$ through linear interpolation of the experimental pressure data.

3. (iii) (4.23) yields a cubic equation, $p_{n}a_{n}^{3}+\unicode[STIX]{x1D70E}a_{n}^{2}-\unicode[STIX]{x1D707}_{n}=0$ . Solve to obtain the radius $a_{n}$ .

4. (iv) Compute $\unicode[STIX]{x1D707}_{n}^{\prime }$ from (4.22) and integrate over $\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}_{m}$ to obtain $\unicode[STIX]{x1D707}_{n+1}$ .

5. (v) Obtain the radius $a_{n+1}$ in the same way as in step (iii) and then compute $a^{\prime }$ .

6. (vi) Compute $\unicode[STIX]{x1D714}_{n+1}$ and velocity components $u_{\unicode[STIX]{x1D709}}$ and $u_{\unicode[STIX]{x1D702}}$ with the streamfunction–vorticity method. First set $\unicode[STIX]{x1D714}_{k=0}=\unicode[STIX]{x1D714}_{n}$ . Define the secondary time step $\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}_{v}$ .

   1. (a) Update $\unicode[STIX]{x1D714}_{k}$ .

   2. (b) Compute $\unicode[STIX]{x1D713}$ from (4.12) using $\unicode[STIX]{x1D714}_{k}$ .

   3. (c) Compute $u_{\unicode[STIX]{x1D709}}$ and $u_{\unicode[STIX]{x1D702}}$ from (4.14).

   4. (d) Compute $\unicode[STIX]{x1D714}_{k+1}$ from (4.16).

   5. (e) Check for convergence of $\unicode[STIX]{x1D714}_{k+1}$ . If tolerance in (4.24) is met, $\unicode[STIX]{x1D714}_{n+1}=\unicode[STIX]{x1D714}_{k+1}$ ; otherwise, update $k$ and return to step (a).

7. (vii) Compute $c_{n+1}$ from (4.19).

8. (viii) Update $n$ and return to step (i).