Growth of a bubble cloud in CO2-saturated water under microgravity

Patricia Vega-Martínez a, Javier Rodríguez-Rodríguez *a and Devaraj van der Meer b
aFluid Mechanics Group, Carlos III University of Madrid, Leganés, Spain. E-mail:
bPhysics of Fluids Group, University of Twente, Enschede, The Netherlands

Received 3rd January 2020 , Accepted 26th March 2020

First published on 30th March 2020

The diffusion-driven growth of a dense cloud of bubbles immersed in a gas-supersaturated liquid is a problem that finds applications in several modern technologies such as solvent-exchange micro-reactors, nanotechnology or the manufacturing of foamy materials. However, under Earth's gravity conditions, these dynamics can only be observed for a very limited time if the cloud is not attached to a surface, due to the action of buoyancy, i.e. of gravity effects. Here, we present experimental observations of the time evolution of dense bubble clouds growing in CO2-supersaturated water under microgravity conditions. We report the existence of three regimes where the bubble cloud exhibits different growth rates. At short times, each bubble grows independently following the Epstein–Plesset equation. Later on, bubbles start to interact with each other and their growth rate diminishes as they compete for the available CO2. When this happens, the growth rate slows down. This occurs earlier the deeper the bubble is in the cloud. Finally, at long times, only those bubbles on the husk continue growing. These regimes may be qualitatively described by a mathematical model where each individual bubble grows in the presence of a constellation of point mass sinks. Despite the model being only valid for dilute bubble clouds, its predictions are consistent with the experimental observations, even though the bubble clouds we observe are rather dense.

1 Introduction

The diffusion-driven dynamics of a dense cloud of bubbles (or drops) have received renewed attention lately due to their applications in modern chemical processes such as solvent extraction1 and in nanoscience.2 Besides these immediate applications, understanding these diffusive dynamics also sheds light on certain processes taking place in foams, such as coarsening or coalescence. Indeed, although many other physical effects are involved in the evolution of a wet foam, still diffusion-related phenomena play a central role.3 However, despite the diffusive interaction between nearby bubbles in bubbly liquids or foams having been studied for a long time, still there is a need for quantitative predictive models.4 In part, this is a consequence of the difficulties that arise in isolating purely diffusive effects in experiments. In fact, the large difference in density between the bubbles and the liquid makes gravity effects manifest rather early in the evolution of foams or bubble clouds.3

In foams, it is hard to disentangle the role of liquid drainage (a gravity-related effect) from other effects associated with the transport of gas5 or anti-foaming agents6 in the rate of bubble coalescence. Moreover, drainage turns wet foams into dry ones very fast, so the former are hard to characterize experimentally.3 Even in a situation where diffusion is the leading effect, such as the growth of a bubble cloud in a gas-supersaturated liquid, gravity effects become of comparable importance very quickly. For example, a free cloud of sub-millimetric bubbles exhibits a purely diffusion-driven growth only for a very limited time (hundreds of milliseconds), even for high gas concentrations such as those found, for instance, in beer or champagne, as shown by Rodríguez-Rodríguez et al.7 In this study, a dense bubble cloud was produced in the bulk of a CO2-saturated liquid (beer) and its evolution was studied using high-speed imaging. This evolution could be divided into three stages: in the first stage, which lasts for only about 10 ms, the size of the bubbles increases with the square root of time, as predicted by the Epstein–Plesset equation for the case of a single isolated bubble.8 Thus, it is reasonable to assume that in this stage, bubbles grow without interacting much with each other. Then, another stage follows where bubbles moderate their growth, due to the depletion of the dissolved gas in the space between bubbles. These two stages, where the dynamics are driven by diffusive bubble growth, come to an end when gravity effects become important and make the bubble cloud rise, which occurs as early as about 100 ms even for a bubble cloud as small as about 1 mm in diameter. Of course, bubbles could be prevented from rising by keeping them pinned to a substrate. However, the pinning introduces additional complexities that obscure, to some extent, diffusive effects. Moreover, in some situations, gravity plays a role even if the bubbles remain pinned to a substrate. In the case of sessile bubbles in CO2-supersaturated water, Enríquez et al.9 and Soto et al.10 have shown that although bubbles are not allowed to move, the depletion of CO2 that takes place around them induces density gradients in the liquid that ultimately trigger a convection plume, which in turn dominates the flow and overcomes the effect of diffusion.

Therefore, avoiding or minimizing gravity effects is essential to study the purely diffusion-driven dynamics of wet foams or bubble clouds. Moreover, besides this purely fundamental interest, there are situations in which bubble clouds or wet foams evolve under microgravity conditions or, at least, under reduced gravity conditions. For example, in the field of planetary physics, the diffusive growth of bubbles is responsible for the loss of helium in meteorites.11 Thus, modeling the evolution of helium content in these objects is important to understand the history of the formation of our solar system. In the Moon, rocks with a highly vesicular structure, such as the “Seatbelt Rock”, have been found (see Fig. 4 in the report by Jolliff and Robinson12). This structure suggests that the rock has formed by solidification of a dense cloud of bubbles containing volatile gases, as a result of ancient volcanism. From a more technological point of view, the behavior of foams under reduced gravity conditions is of interest for the development of advanced manufacturing processes in space, as pointed out by Koursari et al.13 These authors studied experimentally the drying of a wet foam in the absence of gravity. Since on Earth, gravity is essential to drain out the water filling the interstices between bubbles, in space, they use a porous medium to soak up the liquid. Note that if bubbles were allowed to grow by supersaturating the liquid with a gas, their growth could push out the interstitial water without the need of an external absorbing body.

In summary, studying experimentally the diffusive interaction of free growing bubbles in dense bubble clouds or in wet foams at long times requires avoiding the effect of gravity. This has been done in several studies in the case of foams5,13–15 but, to the best of our knowledge, there are no studies focusing on dense bubble clouds growing diffusively in gas-supersaturated liquids under microgravity conditions. With this motivation in mind, we have carried out experiments in which dense bubble clouds grow in CO2-supersaturated water under microgravity conditions. These experiments were performed in the drop tower of the German Center of Applied Space Technology and Microgravity (ZARM) using an improved version of the facility described by Vega-Martínez et al.16 Besides, inspired by these experiments, we have put together a mathematical model for the growth of a bubble cloud in a supersaturated liquid that is able to explain, albeit in a qualitative fashion, the features of the time evolution of bubble clouds observed in the experiments.

2 Experimental setup and procedures

The experiment aims at generating a dense bubble cloud in a supersaturated liquid, in our case CO2-saturated water, and then observing its evolution using high-speed imaging under microgravity conditions. Fig. 1 shows a sketch of the experimental setup, which is based on that described by Vega-Martínez et al.16 This experimental setup is fitted into a capsule that is then dropped inside the drop tower of the German Center of Applied Space Technology and Microgravity (ZARM). This tower has a total height of 146 meters and contains in its interior a 120 m-tube where the capsule falls in near vacuum for about 4 s. This generates a microgravity environment for the duration of the drop. Then, the capsule decelerates and comes to a stop in a pool filled with polystyrene pellets.
image file: d0sm00015a-f1.tif
Fig. 1 Layout of the experimental setup and the drop tower facility. The two high-speed cameras form an angle of 90°. The chamber has a cylindrical shape, 24.4 mm in diameter and 200 mm in height. The bubble cloud is generated at the center of the chamber to avoid the effect of walls. The right panel is a sketch of the capsule inside the drop tower (not to scale).

The measurement chamber, where the bubble cloud will form and grow, is filled with carbonated water. This chamber is immersed in a rectangular tank filled with degassed water to minimize optical distortion. The top wall of the chamber is connected to a pressurization and depressurization system. At the bottom of the measurement chamber, there are two electrodes connected to the spark generator. The electrodes are two thin (100 μm) copper wires that touch each other near their tips. The spark generator is a discharge circuit similar to that described by Willert et al.17 but replacing the LEDs with electrodes, as suggested by Goh et al.18 In our experiments, the discharge has a peak between 50 and 100 A and lasts for about 10 μs. The discharge induces cavitation, and the collapse of the imploding cavitation bubble generates the bubble cloud, the growth of which is the target of the experiment. To measure the time evolution of the shape and size of the cloud, two high-speed cameras, suitable for use in the drop tower's capsule, image the bubble cloud from two perpendicular directions at 2000 fps. They produce images of 512 × 512 pixels and a spatial resolution of 25 μm per pixel. As can be seen in the appendix, the information provided by the two cameras is qualitatively similar for all drops. For that reason, in the discussion of the experimental results, we will only use images taken from camera 1.

The supersaturated water is prepared in a pressurized facility built to make carbonated water on-the-site, which is subsequently transferred to the measurement chamber at a pressure higher than the saturation one, such that no gas is lost during the process. This procedure is similar to that followed by Enríquez et al.19

Fig. 2 follows the evolution of a bubble cloud under microgravity conditions (10−6 g and 1 g) from its generation (the spark) to the end of the experiment, approximately 3.4 seconds later. During the first 100 ms, both tests show a similar diffusion-driven growth. Subsequently, in the 1 g test, the effects of gravity become abundantly clear: the cloud starts a buoyancy-induced rising motion, as can be seen in the first row of the figure. In the microgravity test, however, the bubble cloud continues growing by diffusion and no net motion is observed. For the model developed in the next section, it is important to point out that despite the high bubble density of the cloud, we observe that very few bubbles coalesce during the experiments. Thus, in a first approximation, the number of bubbles may be regarded as constant.

image file: d0sm00015a-f2.tif
Fig. 2 Comparison of experiments where the bubble cloud develops with (upper row) and without gravity (lower row). See Movies in the ESI.

3 Mathematical model of transient mass transfer in a bubble cloud

Let us consider a spherical bubble of radius Ri surrounded by a cloud of other bubbles of radii Rj, with j = 1…N and ji. We introduce here three hypotheses. First, from the point of view of mass transfer, one bubble sees the others in the cloud as mass sinks of intensity j. Second, we do not consider surface tension effects. Third, advective effects caused by bubble growth are neglected in the gas transport equation. This means that the concentration of dissolved gas obeys the non-convective heat equation. Due to the linear nature of this equation, the interaction of the i-th bubble with the surrounding ones can be treated by superposition. For this reason, we consider the growth of a bubble, labelled i, in the presence of a point mass sink of intensity j a distance dij apart from its center (see Fig. 3). The problem we will solve now is to determine the gas mass flux entering bubble i under these conditions. This problem exhibits cylindrical symmetry around the line connecting the centers of both bubbles, so the concentration field obeys
image file: d0sm00015a-t1.tif(1)

image file: d0sm00015a-f3.tif
Fig. 3 Bubble of radius Ri in the presence of a mass sink of intensity j located at a distance dij.

This equation must be subject to boundary conditions at the bubble surface, r = Ri, and in the far field r → ∞. At the bubble surface, the concentration is given by Henry's law, in other words C(R,t) = Cs = kHPs, where kH is Henry's constant and Ps is the ambient pressure, whereas in the far field, C([x with combining right harpoon above (vector)],t) = C. Note that although eqn (1) is linear, the fact that the radius of the i-th bubble, Ri, depends on time effectively couples the computation of the mass fluxes of the different bubbles, as every bubble j affects the concentration field induced by the others through the boundary condition at r = Ri. Nonetheless, as will be explained below, this effect can be neglected with the assumptions adopted here.

The solution of the problem will be sought by exploiting its linear nature to decompose it into a spherically-symmetric one that satisfies these boundary conditions plus another one that accounts for the effect of the point sink but that has homogeneous boundary conditions. The former solution recovers the mass flux given by the Epstein-Plesset equation,8i0. To calculate the latter, denoted by ij, we transform eqn (1) into the Laplace plane. In what folows, Ĉ(r,s) and [m with combining circumflex]j(s) are the Laplace transforms of C(r,t) and j(t), respectively. With this notation, eqn (1) turns into

image file: d0sm00015a-t2.tif(2)
To transform the time derivative of the concentration field, we take into account that with the decomposition explained above, the concentration field induced by an individual mass point sink has homogeneous initial conditions, thus ∂C/∂t transforms to .

The general solution to this equation is:20

image file: d0sm00015a-t3.tif(3)
where h(1)l and h(2)l are Hankel's functions of the first and second kind, respectively. Note that since they contain exponential functions of imaginary exponents, the minus sign inside the square root of their arguments guarantees that they take real values. Although to compute the concentration field, it is necessary to determine the full set of coefficients (al,bl), the calculation of the mass flux [m with combining circumflex]ij that sink j induces on the bubble only requires knowledge of the coefficients a0 and b0. Indeed, the mass transfer due to higher-order harmonics is zero, as
image file: d0sm00015a-t4.tif(4)
with the last expression being zero for l ≥ 1. Moreover, coefficient b0 must vanish, as the function image file: d0sm00015a-t5.tif is unbounded as r → ∞.

To obtain a0, we introduce (3) into (2) and project onto P0(cos[thin space (1/6-em)]θ) to find

image file: d0sm00015a-t6.tif(5)
where the variable change x = cos[thin space (1/6-em)]θ has been introduced. Here, it is good to note that the expression of eqn (5) is a hybrid, as it mixes expressions in the Laplace space with the radius of bubble i, which generally is a function of time, Ri(t). We will, however, ignore this fact and treat Ri as a constant in the Laplace space to make analytical progress. Physically, this amounts to saying that we assume that the time scale with which the concentration field responds to changes in the configuration is much smaller than that at which the bubble grows. This is an assumption that is quite similar to the “frozen bubble” assumption that was employed in solving the single bubble problem in the original paper by Epstein and Plesset.8 Evaluating the integrals, and taking into account image file: d0sm00015a-t7.tif,
image file: d0sm00015a-t8.tif(6)
image file: d0sm00015a-t9.tif(7)

Introducing this value of a0 together with b0 = 0 into (3) and differentiating, we get

image file: d0sm00015a-t10.tif(8)

By transforming this expression back into the time domain and introducing it in the definition of the mass flux (4),

image file: d0sm00015a-t11.tif(9)

To obtain the total mass flux into bubble i, this expression must be summed over all bubbles in the cloud (with the exception of i itself) and then added to the mass flux provided by the Epstein–Plesset equation, i0, to yield

image file: d0sm00015a-t12.tif(10)

To make the above expressions dimensionless, we introduce the following variables (suppressing indices), a = R/Rc, Δ = d/Rc, τ = Dt/Rc2 and = dM/dτ = /(ρgDRc), where Rc is a fixed characteristic bubble diameter, and ρg is the gas density at the ambient pressure Ps. Besides, the saturation level, which indicates the amount of CO2 dissolved in the liquid, is defined by ζ = C/Cs. Let us define also the Ostwald coefficient as Λ = Cs/ρg. With these definitions, eqn (10) results in

image file: d0sm00015a-t13.tif(11)

This mass flux must be completed with the mass conservation equation for the i-th bubble, namely

image file: d0sm00015a-t14.tif(12)

The changes in the volume of the bubble induce a velocity field, which, sufficiently far away from its center, can be modelled as that of a volume point source. Thus, considering the whole cloud, each individual bubble moves as a result of the superposition of these flows. In the Stokes limit, the resulting bubble velocity can be computed using the solution provided by Michelin et al.4

image file: d0sm00015a-t15.tif(13)
In this equation, ēij is the unit vector pointing from the center of the i-th bubble to that of the j-th bubble. In this work, we are going to let bubbles translate with only the first term of this expression since we expect (aj/Δjl)3 ≪ 1.

3.1 Numerical method

The computation of the time evolution of the bubble cloud is divided into two stages. First, we integrate eqn (11) and (12) together assuming that the bubbles do not displace and then, in a second stage, we compute the motion of the bubbles using the velocity field given by (13).

The system of integro-differential eqns (11) and (12) is solved using an explicit Euler method for the temporal derivatives and a trapezoid method for the integral. Note that the kernel of the integral equation tends smoothly to zero at the integration limit τ′ = τ, which makes the usage of more sophisticated integration methods unnecessary. We can distinguish two terms in eqn (11): the first one is the mass flux provided by the Epstein–Plesset equation,8 which corresponds to the evolution of an isolated bubble. The second term, the integral, models the interaction between each individual bubble and the rest of the cloud. For the calculation of this second term, we will adopt the simplifying assumption that the radius of the bubble corresponds to the initial one (frozen bubble). Note that this hypothesis is also adopted in the derivation of the Epstein–Plesset equation.8 The time step chosen for the simulations presented in this paper is Δτ = 10−3, which has been selected after verifying that taking a smaller time step does not affect the results.

Once the time evolution of the bubble radii is computed, the bubbles are displaced using the velocity field given by (13) using also, for consistency with the previous stage, an Euler method with the same time step.

The numerical computation is started by integrating solely the Epstein–Plesset equation up to the time step used for the rest of the calculation, Δτ. Note that this is possible since at very short times, the concentration boundary layer is confined to a thin shell around the bubble of thickness image file: d0sm00015a-t16.tif, which is much smaller than the typical bubble–bubble distance. This makes it possible to neglect the integral term in this first step. This strategy is also used in the computation of the bubble velocities.

4 Results and discussion

We split this section into two parts. First, we explore the predictions of the theoretical model just described for conditions representative of our experiments. Then, we compare these predictions, in a qualitative way, with those features of the evolution of the bubble clouds that could be determined from the experiments.

4.1 Results of the model

In this section, we show the predictions of the model for bubble clouds representative of those observed in the experiments. Since the number of dimensionless parameters of the problem is large, we restrict our results to initially spherical bubble clouds with a number of bubbles large enough to exhibit collective effects, namely N = 100, and evolving in liquids with different saturation levels, ζ. Furthermore, all the bubbles have the same initial size, R0 = 30 μm. This radius is consistent with the estimated size of the bubble fragments that result in a similar cavitation experiment.7 Note that for bubbles of this size, it is reasonable to neglect the effect of surface tension in their dynamics, as the Laplace pressure, 2σ/R0, with σ ≈ 0.7 N m−1 as the surface tension coefficient, only becomes of the order of ambient pressure for bubbles of size 2σ/Ps ≈ 1.4 μm. In all the calculations, we have adopted Λ = 0.842 and D = 1.92 × 10−5 m2 s−1, corresponding to the values of the system CO2–water at 25 °C.21 The initial positions of the bubbles are chosen randomly within a sphere, the cloud, of radius Rbc, such that the average bubble density is approximately uniform. To clearly identify the moment when bubbles start to interact with each other, it is convenient to seed the bubbles in the cloud with an average inter-bubble nearest neighbor distance, d, larger than, but still of the order of, one bubble diameter. Assuming that the volume of cloud available per bubble is Vb ∼ πd3/6, setting Rbc ≈ 10R0 yields d ≈ 4.3R0, which fulfills the above criterion. In what follows, all length scales are made dimensionless with R0 and time scales with R02/D, using the same notation as in the previous section.

Fig. 4 shows qualitatively the evolution of a modeled bubble cloud with a saturation ζ = 2. The top left side is the initial state, τ = 0, whereas the top right one corresponds to the cloud at a later time, τ = 3.2. The color indicates the initial distance from the bubble to the center of the cloud, growing from light to dark hue. To illustrate this color coding, the bottom image represents a cut of the cloud at τ = 3.2. To average results for bubbles located at a given initial distance from the center, the cloud radius in the initial state has been divided into B = 10 layers equally spaced along radius Rbc, with the color indicating the layer number. Using this averaging strategy, we can show how bubbles grow at different depths inside the cloud.

image file: d0sm00015a-f4.tif
Fig. 4 Three-dimensional view of a bubble cloud generated with the model described in Section 3, (a) at time τ = 0 and (b) at a later time, τend = 3.2, in dimensionless units. Panel (c) represents half the cloud, consisting of only those bubbles with centers obeying x < 0, at the same time. The saturation of the liquid is ζ = 2. The color represents the initial distance of each bubble to the center of the cloud, growing from lighter to darker colors.

Panels (a, c and e) of Fig. 5 show the time evolution of the radii of the bubbles in three clouds with saturation levels ζ = 1.1, 2 and 4, averaged over each radial layer. Correspondingly, panels (b, d and f) show the average bubble radius of each layer at the end of the simulation (τ ≈ 6.5). The duration of the simulation has been chosen to allow for the observation of the different growth regimes. At early stages, all the bubbles grow freely, that is, as isolated bubbles. In this situation, their growth is well described by the Epstein–Plesset equation,8 corresponding to the red dashed line in the figure. Using our notation, this equation reads

image file: d0sm00015a-t17.tif(14)
The departure from the growth predicted by this equation is a consequence of the fact that the thickness of the diffusive boundary layer around the bubble grows as image file: d0sm00015a-t18.tif. Therefore, until times of order τ ∼ (Δ/2 − 1)2, the boundary layers of the different bubbles do not interact and bubbles do not feel each other. However, from that time onwards, the diffusive boundary layers overlap and bubbles compete for the available CO2. As a consequence, their growth rates diminish quickly, even becoming zero or slightly negative for bubbles in the innermost layers, where the access to the CO2 present in the bulk outside the cloud is limited. In contrast, the growth of the bubbles in the husk of the cloud continues, albeit at a much smaller pace. Note that the freezing of the growth of the innermost bubbles, which is also clearly visible in Fig. 4c, is consistent with the experimental observations (see Fig. 9 and the associated discussion). Regarding the slight reduction of the bubble sizes for bubbles deep inside the cloud for τ ≳ 2, observed in Fig. 5, although it seems reminiscent of Ostwald-rippening, where surface tension causes large bubbles to grow at the expense of small ones,22 surface tension is not considered here. Consequently, this bubble size reduction cannot be attributed to this effect. In fact, in the absence of surface tension, the gas concentration at the bubble surface could never be smaller than the saturation one, Cs, which means that bubbles should never shrink. The cause of this anomalous effect is the fact that we model nearby bubbles as point sinks. For a point sink, imposing a given mass flux means that the concentration is minus infinity at the sink. So, if the sink is very close to the bubble, it may occur that locally, the concentration becomes smaller than the saturation one, which is a non-physical effect. Naturally, this would not occur where the sink is replaced by a physical bubble, as in its surroundings, the concentration follows CCs at all times. Nonetheless, the bubble size reduction created by this simplification is at most about 15% in the case ζ = 4 and smaller in other cases.

image file: d0sm00015a-f5.tif
Fig. 5 (a, c and e) Dimensionless time evolution of the bubble radii for saturation levels ζ = 1.1, 2 and 4, respectively. The solid-line shows the mean radius for each cloud layer, with the color evolving from lighter to darker according to the distance of the layer to the cloud center. The error bars mark the maximum and minimum radius of the bubbles for each layer at that time instant. The red-dashed line represents the growth of an isolated bubble predicted by the Epstein–Plesset equation.8 (b, d and f) Final layer-averaged radius of the bubbles in the cloud as a function of its distance to the center. Markers represent the mean final radius, whereas the vertical bars show the maximum and minimum value of the radius. The horizontal error bars correspond to the maximum and minimum final distance to the cloud center of the bubbles contained in each layer. The black horizontal dashed line gives the value of af predicted by eqn (15).

The final radius of bubbles belonging to the innermost layers can be estimated using simple mass conservation arguments. Since these bubbles lose access to the CO2 from the bulk liquid, they almost stop growing when they have absorbed all the excess CO2 dissolved in the liquid volume of the cloud available per bubble, which we can regard as that of a sphere of diameter d. Thus, the final dimensional mass, 4π/3Rf3ρg, must be the sum of the initial one, 4π/3R03ρg, plus the excess mass of CO2 dissolved in the liquid volume with respect to the saturation one, namely (CCs)π/6d3. In dimensionless terms, this yields for the final radius

af = (1 + Λ(ζ − 1)Δ3/8)1/3.(15)
For the saturation levels corresponding to the simulations in Fig. 5, namely ζ = 1.1, 2 and 4, the values of af = 1.21, 2.22 and 3.14, respectively. We must point out here that when we refer to the final radius or bubble sizes, we mean at the end of the simulation. Note that, even at the center of the bubble cloud, still bubble radii keep changing with time at long times. However, they do so at a much smaller pace than that predicted by the Epstein–Plesset equation, thus we find it reasonable to regard them as quasi-frozen.

As can be observed in Fig. 5b, d and f, the agreement with the numerical results is fairly reasonable, given the assumptions adopted in this estimation. In this figure, markers denote the layer-averaged final radius, with the vertical error bars indicating the minimum and maximum bubble size for each layer. In an analogous way, the horizontal position of each marker represents the mean radial position of the bubbles in the layer at the end of the simulation, which are all contained in the segment marked by the horizontal error bars.

Once bubbles deep into the cloud reach this maximum size, it seems reasonable to assume that only those in the outer shell continue growing. Because there is almost no CO2 inside the cloud, all the mass flux comes across the outer boundary of the cloud. Moreover, this mass flux must be shared between all the bubbles in the outer shell. In a first approximation, we can model this effect by assuming that each bubble no longer absorbs gas from a solid angle 4π (that is, from all directions) but from a solid angle 4π/Nos, where Nos is the number of bubbles in the outer shell. This means that the Epstein–Plesset equation can still be used to model bubble growth in the husk, by dividing its right hand side by the number of bubbles in the outer shell. To check this hypothesis, we plot in Fig. 6 the time derivative of the radius of the bubbles in the outermost layer at the end of the simulation divided by the one predicted by the Epstein–Plesset equation as a function of the inverse of the number of bubbles in this layer. These results correspond to the same saturation level, ζ = 2, but to four different realizations. Because in each realization, the bubbles are distributed randomly, their number in the outermost layer varies. The average value of this number can be estimated by assuming that the bubble density in the cloud is constant and that the outermost layer has a thickness of about d/2. Using simple geometrical analysis, this yields Nos ≈ 3Nd/2Rbc.

image file: d0sm00015a-f6.tif
Fig. 6 Rate of change of the radius of the bubbles in the outermost layer divided by that predicted by the Epstein–Plesset equation, as a function of the inverse of the number of bubbles in this layer. The solid line is 1/Nos.

We conclude this section with a comment about the effect of the saturation level on the growth rate of the cloud. Although in the experiments, the radius of the individual bubbles cannot be determined experimentally, it is possible to measure the time evolution of the projected area of the whole cloud by applying image processing to each frame of the high-speed movies. To compare these results against the simulations, we need to compute the time evolution of the projected area of the simulated cloud. To that end, we create synthetic images projecting the three-dimensional bubble cloud onto a plane, as sketched in Fig. 7 and shown in Fig. 4. Using this technique, we are able to take into account the effect of bubbles overlapping on the projected cloud area. Note that, since the cloud is spherical and the bubbles are distributed randomly, the choice of the projection plane does not have an effect on the results.

image file: d0sm00015a-f7.tif
Fig. 7 Schematic of the definition of the projected area.

The physics behind the growth of the bubble cloud becomes more clear if we focus on the time derivative of the projected area rather than on the area itself, as we explain in the next subsection. Fig. 8a shows how, for ζ ≥ 2, the evolution of both the projected area and the total gas volume is qualitatively independent of this value, just affecting the absolute value of the growth rates.

image file: d0sm00015a-f8.tif
Fig. 8 (a) Time derivative of the measured projected area and the measured 2/3 root of the volume for numerically created bubble clouds with identical initial configurations but different saturation levels, ζ. (b) Time derivative of the projected area divided by (ζ − 1). Note that to facilitate the comparison with the experiments, we show here the numerical results in dimensional form, using R0 = 30 μm and R02/D ≈ 0.47 s as length and time scales, respectively.

Indeed, examination of the Epstein–Plesset eqn (14) suggests that, provided the radius of the bubble is sufficiently larger than the initial one, we can assume image file: d0sm00015a-t19.tif. In this case, the growth rate of a2, which is roughly proportional to the projected area, scales with (ζ − 1). To check this scaling, we plot in Fig. 8b the time derivative of the dimensional projected area divided by the saturation level, ζ − 1. As predicted, the curves collapse fairly well at intermediate and long times. At short times, because the bubble radii are still close to the initial ones, the collapse is not so good. Furthermore, the curve corresponding to the smallest saturation level, ζ = 1.1, does not exhibit such a good collapse at any time. This is due to the influence of the initial bubble radius, which lasts longer due to the smaller growth rate.

The model predicts the existence of three different stages in the evolution of the cloud, which, despite their noise and variability, are also observed in the experiments. At short times, the area grows relatively fast, at a nearly constant rate. Later on, the rate of change of the area decays slightly faster than the square root of the time. Finally, at long times, the area grows again at a constant rate, albeit at a much slower pace.

The last two regimes can be easily explained with a simple heuristic model. Let us assume that the volume of the cloud is the sum of that of N bubbles plus their associated water volume (a sphere of radius d/2 for each bubble). Since the water volume in the cloud does not change, the radius of such a cloud would then be

image file: d0sm00015a-t20.tif(16)
where, for simplicity, we take all bubbles to have the same size. The projected area, Abc = πabc2, has a rate of change
image file: d0sm00015a-t21.tif(17)
At short times, the radius of the bubble does not differ much from its initial one, a ≈ 1, whereas its time derivative goes as ȧτ−1/2, as predicted by the Epstein–Plesset equation. Introducing this into eqn (17), we get dAbc/dττ−1/2. This is in fairly good agreement with the behavior observed at intermediate times in both experiments and theory. Actually, the rate of change of the area seems to decay slightly faster than this prediction, which is consistent with the fact that bubbles deep inside the cloud stop growing once they deplete their surroundings of CO2. In other words, the effective value of N in eqn (17) also decays with time, until it only refers to those bubbles located in an outer shell of thickness of order Δ.

At long times, the radius of the bubbles in the outer shell still grows as aτ1/2, albeit with a prefactor smaller than that predicted by Epstein–Plesset, as discussed above. Moreover, ad/2. In this limit, eqn (17) yields dAbc/dτ, which is constant. As we will see in the next subsection, this is indeed what is observed to occur in the simulations and is qualitatively consistent with the experiments as well, despite the noise.

At short times, the model predicts a fast growth rate, not only in absolute value, which is to be expected as a consequence of bubbles growing without competing, but also in terms of how slowly it decays. Indeed, although eqn (17) suggests that dAbc/dττ−1/2 at short times as well, we observe that dAbc/dτ is rather constant instead. This is a consequence of the fact that the void fraction inside the cloud grows fast at short times. Thus, the growth of the projected area is due to the addition of two effects: its boundaries grow, as predicted by eqn (17), and the interior of the projected cloud image becomes more opaque as more bubbles overlap on the projection. Naturally, if we compute the total gas volume inside the cloud, V, we observe a much closer agreement with dV2/3/dtτ−1/2, since the fact that bubbles overlap in the image does not affect the evolution of the volume. This can be observed in Fig. 8, where the evolution of the time derivative of V2/3 is compared with that of the projected area of cloud A for simulations with identical initial bubble distributions but different saturation levels.

4.2 Comparison with experiments

In this section, we analyze the evolution of the size of the bubble clouds observed in our microgravity experiments, comparing it with that predicted by the model. Several key features of the bubble cloud, such as the number of bubbles or their initial size and spatial distributions, cannot be determined experimentally. For this reason, the comparison will be predominantly of qualitative nature and restricted to the area of the bubble cloud projected on an image, which is the magnitude that can be measured experimentally. Nonetheless, despite this limitation, we will show that the different regimes and growth rate behaviors predicted by the model are consistent with those observed experimentally.

An important aspect is to determine the saturation level of CO2 of the water in the experiments. Some isolated bubbles appear in almost all the experiments, and we use them to determine the saturation level, ζ, as described in ref. 16. To this end, we track the evolution of their radii by image processing. Then, we estimate ζ by fitting that evolution to that predicted by the Epstein–Plesset equation,8 obtaining values between 3.5 and 4.3, depending on the experiment.

As stated above, the high bubble density of the cloud precludes the observation of its interior. Consequently, it is not possible to determine experimentally the number of bubbles nor their size distribution, at least not accurately. However, as soon as the capsule breaks upon reaching the bottom of the tower, the sudden appearance of a large apparent gravity allows us to instantly observe those bubbles that were at the center of the cloud (Fig. 9). This observation reveals that bubbles deep inside the cloud are much smaller than those near the outer edge, as predicted by the model (Fig. 5).

image file: d0sm00015a-f9.tif
Fig. 9 Snapshot of a high-speed movie showing the cloud during the deceleration of the capsule (ζ ≈ 4). As the capsule comes to a stop, the large deceleration makes the bubble cloud rise fast. Because large bubbles rise faster, this allows us to see those that were inside the cloud (within the red ellipse region). See also Movies in the ESI.

Fig. 10 shows the time evolution of the projected area for the six drops. In this and the following plot, all the magnitudes are plotted with dimensions, to give an idea of the orders of magnitude of the spatial and temporal scales. The jumps observed in some of the curves are due to the merging or separation of individual bubbles from the main cloud. Besides the fact that all the curves grow monotonically, it is not easy to extract any trend, especially not the existence of different regimes. In part, this is a consequence of the large dispersion in the growth rates of the projected area, which arises from the variability in the initial conditions of the bubble cloud caused by the collapse of the spark-induced cavitation bubble. But also, even if all the clouds were generated identically, the projected area at a given instant is an integral of all the previous growth rates. Thus, this magnitude is not suitable to identify different stages of the cloud growth rate. For this reason, we focus our analysis on the time derivative of the projected area, which is an instantaneous magnitude, i.e. not depending on the bubble cloud history or initial conditions.

image file: d0sm00015a-f10.tif
Fig. 10 Evolution of the measured projected area of six bubble clouds. The level of saturation of the water in these experiments varies between 3.5 and 4.3.

Fig. 11 shows on the logarithmic scale the rate of change of the area of the cloud for the six drops studied in the campaign (see the appendix for details on how the time derivative of the area is calculated). The curves in the figure have been scaled with their value at the center of the time range where they follow the law dA/dtt−1/2 to be able to compare the evolution of clouds with very different initial sizes and thus growth rates. We plot also the rate of change of the area computed from a numerical simulation with ζ = 4.

image file: d0sm00015a-f11.tif
Fig. 11 Rate of change of the projected area for the six drops, compared with that predicted by a simulation with ζ = 4. The curves have been scaled with their corresponding value at the center of the time range where they follow the scaling law dA/dt = Ct−1/2, that is, tr = 0.4 s. Here, C is the fitting constant of each experiment.

In the experiments, the behavior at short times exhibits a larger variability than that predicted by the model. In some experiments, we observe the trend dA/dtt−1/2 from the first instants, whereas in others, the projected area shows an initially constant growth rate. We attribute this anomalous behavior to the fact that, at such short times, some of the clouds still retain some residual velocity as a result of the implosion of the cavitation bubble that birthed them. As this initial velocity dissipates, the projected area transitions to the regime where it decays as the inverse of the square root of the time. This residual velocity can be observed in the high-speed movies (see Movies in the ESI). Finally, despite the experimental noise, it is possible to infer that the growth rate stabilizes and decays slower than t−1/2, consistent with the model predictions.

5 Conclusions

We report here experiments where dense bubble clouds grow by diffusion in a supersaturated liquid under microgravity conditions. The bubble clouds are created upon the implosion of a cavitation bubble generated by a spark. Although this way of producing the clouds introduces some variability in the experiments, the evolution of the six bubble clouds studied is qualitatively identical, which suggests that the diffusive growth mechanism that we study is fairly independent of the particular size or shape of the cloud. To achieve microgravity conditions, the experiments have been carried out in the drop tower facility of the German Center of Applied Space Technology and Microgravity (ZARM). The absence of gravity allows us to observe the purely diffusive growth dynamics of the bubble for several seconds, more than an order of magnitude longer than what could be observed under normal gravity conditions,7 where bubble buoyancy would quickly become dominant. Moreover, we avoid other effects such as natural convection9,10 inside the liquid phase.

Inspired by the experimental observations, we have developed a mathematical model where each bubble grows competing for the available CO2 with the others in the cloud, which are treated as point mass sinks. Although, strictly speaking, the model is only valid in the limit of a dilute bubble cloud, it predicts reasonably well the qualitative evolution observed in experiments. The reasonably good agreement of the model with the observations has allowed us to explain, using simple heuristic arguments, the experimental observations. Moreover, this model has been used to investigate aspects that could not be determined experimentally, namely the final distribution of bubble sizes inside the cloud and the effect of the saturation level.

Conflicts of interest

There are no conflicts to declare.

Appendix A

Experimental data from all drops

We present the post-processing analysis carried out for the high-speed movies obtained for the six drops. Regarding the image analysis, the projected area of the cloud has been tracked using custom-made image processing software implemented in Matlab. We started tracking this area from 2 or 3 frames after the spark generates the bubble cloud untill the end of the experiment. Thus, the time span during which the area is measured runs from 0.05 to 3.2 s.

As we said in Section 4, it is easier to observe the physics behind these experiments if we represent the time derivative of the projected area. To calculate the time derivative of the cloud size, we use a central finite difference with a 4th order of accuracy with a uniform grid spacing, Δt = 0.001 s, despite our time experimental resolution being 5 × 10−4 s, to avoid the numerical noise. Then, we screen the data using a ten-point centered moving mean filter in order to see the growth trends. The figure below shows the evolution of the projected area (first column) and the filtered time derivative of the two perpendicular views for the six drops (Fig. 12).

image file: d0sm00015a-f12.tif
Fig. 12 Projected area and the time derivatives obtained from the two camera views in the six experiments carried out during the drop tower campaign.


This work was supported by the Netherlands Center for Multiscale Catalytic Energy Conversion (MCEC), an NWO Gravitation programme funded by the Ministry of Education, Culture and Science of the government of the Netherlands and the Spanish FEDER/Ministry of Science, Innovation and Universities-Agencia Estatal de Investigación though grants DPI2017-88201-C3-3-R and DPI2018-102829-REDT, partly funded through European Funds. We also acknowledge the support of the European Space Agency (ESA) for providing access to the drop tower through grants HSO/US/2015-29/AO and HRE/RS-PS/2018-7/AO. We thank the team from the ZARM Drop Tower Operation and Service Company (ZarM FAB mbH) for their valuable technical support during the experimental campaigns, and the technical staff of the Department of Thermal and Fluids Engineering of Universidad Carlos III de Madrid for their assistance in developing the experimental setup. We are indebted to Dr L. Champougny for her valuable comments about the manuscript.


  1. S. Peng, V. Spandan, R. Verzicco, D. Lohse and X. Zhang, J. Colloid Interface Sci., 2018, 532, 103–111 CrossRef CAS PubMed.
  2. X. Zhu, R. Verzicco, X. Zhang and D. Lohse, Soft Matter, 2018, 14, 2006 RSC.
  3. D. Langevin, C. R. Mec., 2017, 345, 47–55 CrossRef.
  4. S. Michelin, E. Guérin and E. Lauga, Phys. Rev. Fluids, 2018, 3, 043601 CrossRef.
  5. F. Garcia-Moreno, E. Solórzano and J. Banhart, Soft Matter, 2011, 7, 9216–9223 RSC.
  6. P. Yazhgur, D. Langevin, H. Caps, V. Klein, E. Rio and A. Salonen, npj Microgravity, 2015, 1, 15004 CrossRef PubMed.
  7. J. Rodríguez-Rodríguez, A. Casado-Chacón and D. Fuster, Phys. Rev. Lett., 2014, 113, 214501 CrossRef PubMed.
  8. P. S. Epstein and M. S. Plesset, J. Chem. Phys., 1950, 18, 1505–1509 CrossRef CAS.
  9. O. R. Enríquez, C. Sun, D. Lohse, A. Prosperetti and D. van der Meer, J. Fluid Mech., 2014, 741, R1 CrossRef.
  10. A. M. Soto, O. Enríquez, A. Prosperetti, D. Lohse and D. van der Meer, J. Fluid Mech., 2019, 871, 332–349 CrossRef CAS.
  11. F. Stuart, P. Harrop, S. Knott and G. Turner, Geochim. Cosmochim. Acta, 1999, 63, 2653–2665 CrossRef CAS.
  12. B. L. Jolliff and M. S. Robinson, Phys. Today, 2019, 72, 44–50 CrossRef.
  13. N. Koursari, O. Arjmandi-Tash, A. Trybala and V. M. Starov, Microgravity Sci. Technol., 2019, 31, 589–601 CrossRef.
  14. F. García-Moreno, S. Tobin, M. Mukherjee, C. Jiménez, E. Solórzano, G. V. Kumar, S. Hutzler and J. Banhart, Soft Matter, 2014, 10, 6955–6962 RSC.
  15. F. García-Moreno, P. H. Kamm, T. R. Neu, F. Bülk, R. Mokso, C. M. Schlepütz, M. Stampanoni and J. Banhart, Nat. Commun., 2019, 10, 1–9 CrossRef PubMed.
  16. P. Vega-Martínez, J. Rodríguez-Rodríguez, D. van der Meer and M. Sperl, Microgravity Sci. Technol., 2017, 29, 297–304 CrossRef.
  17. C. Willert, B. Stasicki, J. Klinner and S. Moessner, Meas. Sci. Technol., 2010, 21, 075402 CrossRef.
  18. B. Goh, Y. Oh, E. Klaseboer, S. Ohl and B. Khoo, Rev. Sci. Instrum., 2013, 84, 014705 CrossRef CAS PubMed.
  19. O. R. Enríquez, C. Hummelink, G.-W. Bruggert, D. Lohse, A. Prosperetti, D. van der Meer and C. Sun, Rev. Sci. Instrum., 2013, 84, 065111 CrossRef PubMed.
  20. S. F. Wu, The Helmholtz Equation Least Squares Method, Springer, 2015, pp. 27–62 Search PubMed.
  21. P. Peñas-López, A. M. Soto, M. A. Parrales, D. van der Meer, D. Lohse and J. Rodríguez-Rodríguez, J. Fluid Mech., 2017, 820, 479–510 CrossRef.
  22. J. Schmelzer and F. Schweitzer, J. Non-Equilib. Thermodyn., 1987, 12, 255–270 Search PubMed.


Electronic supplementary information (ESI) available. See DOI: 10.1039/d0sm00015a

This journal is © The Royal Society of Chemistry 2020