Hierarchical bubble size distributions in coarsening wet liquid foams

We have investigated foam coarsening for liquid fractions between 15% and 50% using the instrument described in ref. 33. Details can be found in the Materials and Methods. The evolution of the bubble size has been studied elsewhere (34, 35). From the sample surface observations (a typical image is shown in Fig. 1A), we measure the bubble sizes using image analysis and determine the bubble size distributions of the radius normalized by its average $\rho =R/⟨R⟩$. The initial size distributions produced by our experimental setup are asymmetric (positive skew) with a maximum at $\rho \approx 0.6$ (see Fig. 1B for foam with $15\%$ liquid fraction as an example). The normalized size distributions broaden with time, and a sharp peak builds up progressively for small bubble sizes, i.e., $\rho \approx 0.3$, until a stationary form is reached, indicating a Scaling State. This is shown in Fig. 1B for times $t>$ 2,000 s. This evolution is typical of the measurements we have made for foams with liquid fractions within the range $15\%\le \varphi <{\varphi }^{\ast }$, with ${\varphi }^{\ast }\approx 39\%$. The small bubbles corresponding to the peak in the distribution are highlighted in Fig. 1A. After an increase in the transient regime, they finally represent about $35\%$ of the total bubble population in the scaling state (inset of Fig. 1B). We also measured the number of those small bubbles per foam vertex to reach a maximum average value of 1.5, due to space limitation in the vertices. As a consequence, the size distribution becomes invariant in time (statistically self-similar) as observed.

Up to now, such an excess of small bubbles has not been reported in the literature (21, 25, 27, 29). In order to check whether distributions with an excess of small bubbles are also found in drier foams, we have performed coarsening experiments using the same surfactant and a liquid fraction of $8\%$, low enough for gravity effects to be compensated in a ground-based experiment by rotating the cell around a horizontal axis (clinostat). As shown in Fig. 1C, we observed a similar excess of small bubbles. The small bubbles were thus seemingly present but not detected in previous studies. This is probably because high spatial resolution together with a careful image analysis is needed (36). The only experimental work we have found that indirectly relates to this is that of Feitosa and Durian (25), which reports the development of transient bidispersity for initially monodisperse bubbles in a steady state column, where drainage and coarsening occur simultaneously. In their simulations of 2D foam coarsening, Khakalo et al. (37) have observed an excess of small bubbles but the gas transfer through interstitial bulk liquid was not taken into account. Other peculiar size distributions such as “Apollonian” distributions were observed during the decay of beer foam (38) and with emulsions (39). In contrast to what happens in our systems, they arose from coalescence events.

For $\varphi >{\varphi }^{\ast }$, we have observed a different scenario: the initial bubble size distribution shrinks until a steady state is reached where the size distribution is notably narrow (SI Appendix, Fig. S1). The latter distribution is reminiscent of the theoretical distribution predicted for the Ostwald regime (13). Around ${\varphi }^{}\ast$, a change in the growth laws for the average bubble size is also shown for the same foam samples (34, 35):

The Sauter mean radius $R_{32}=⟨R^3⟩/⟨R^2⟩$ is defined as the ratio of third to second moments of the bubble radius distribution.

To clarify the origin of the hierarchical bubble population, we have identified bubbles that eventually disappear and tracked the evolution of their area. Fig. 2A shows examples of such measurements in a foam with $\varphi =15\%$. Similar data are shown for other liquid fractions between $\varphi =20\%$ and $\varphi =33\%$ in SI Appendix (cf. SI Appendix, Fig. S2). Over time, the individual bubble area can either increase or decrease, depending on the bubble’s gas exchanges with its neighbors, but most of the observed bubbles eventually shrink (Fig. 2A). The magnitude of the shrinking rate appears to be initially similar to that characterizing the initial growing rate. Then, a transition occurs and the area decreases much more slowly. Actually, the shrinking after this transition can be extremely slow, and we think this is the underlying mechanism explaining why a peak at smaller than average bubbles builds up in the size distribution. Remarkably, the bubble radius at the transition, $R_t$, is such that its area $A_t=\pi R_t^2$ increases linearly with time, which is similar to the evolution of the squared mean radius in the Scaling State (Eq. 1). Moreover, the transition to the very small shrinking rate appears to occur when the bubble has become so small that it fits inside the interstice between three larger bubbles at the surface and possibly loses contacts with them as sketched in Fig. 2B. (Movies S1–S3). They are free to move throughout the interstice without being pressed against multiple neighbors. Such small bubbles can have different configurations in the interstice, i.e., near the center of the interstice or in contact with one bubble or two bubbles, but these configurations do not last for the entire life of the bubbles because their positions are jostled as the jammed bubbles intermittently rearrange due to the coarsening-induced dynamics (9, 40). We call the small bubbles roaming bubbles. Note that they are reminiscent of rattlers (grains carrying no force) in granular media (41). We conjecture that the bubble size at the transition, $R_t$, should scale as the maximum radius of a sphere that can be trapped in such an interstice at the wall surface.

In a coarsening foam that has reached the Scaling State, there is only one independent length scale of the bubble packing structure. Since the bubbles that form the interstices are bigger than the encaged roaming bubbles, we chose to characterize their average size by the Sauter mean radius. With respect to $⟨R⟩$, $R_{32}$ indeed represents mainly the average radius of the larger bubbles of the distribution and minimizes the contribution of the small bubbles. At a time $t$, the maximum radius of a sphere trapped in such a vertex can be written, on average:

where $x_n(\varphi )$ is a dimensionless geometrical coefficient. We show in SI Appendix, Fig. S3 the plots of $R_t$ versus $R_{32}$ for each liquid fraction. The plots are reasonably described by Eq. 3, allowing the determination of the average coefficient $x_n$ for each liquid fraction (Fig. 2C). $x_n(\varphi )$ varies from $0.25$ to $0.55$ as $\varphi$ varies from $15\%$ to $38\%$, respectively. Using those $x_n$ values, the transition radii $R_t$ collapse on a linear master curve when plotted versus $x_n(\varphi )R_{32}$ (cf. SI Appendix, Fig. S3).

We have performed a geometrical calculation of the size of the interstice between a plane and three perfect spheres of equal radius $R_{32}$ in contact together and with the plane (Fig. 2C). This leads to $x_n=1/3$. This value is smaller than what is measured for liquid fractions corresponding to the bubble random close packing fraction, i.e., ${\varphi }_{\text{rcp}}\approx 31\%$ (see Section D for more details), beyond which the bubbles are spherical. As the liquid fraction gets close to ${\varphi }_{\text{rcp}}$, the foam osmotic pressure, which pushes neighboring bubbles against each other at contacts, becomes very low, and it can be inferred that the cage formed by the triplets of bubbles of radius $R_{32}$ loosens. Note that such a geometrical loosening effect is general and independent of friction (42). Therefore, as a correction to the previous calculation, a distance $ϵR_{32}$ is added around each sphere (Fig. 2C). The coefficient now reads: $x_n=\frac{ϵ(2+ϵ)+\frac{4}{3}{(1+ϵ)}^2}{2(2+ϵ)}\approx \frac{1}{3}+ϵ$, and it increases significantly due to the loosening effect: assuming a moderate loosening $ϵ\approx 0.2$ gives $x_n\approx 0.5$ which is in better agreement with our measurements (Fig. 2C and Movie S2). It is reasonable to assume that polydispersity may also impact the size of the interstice. This effect can be estimated by considering two bubbles of size $R_{32}$ and a third one with size $\beta R_{32}$. It can be shown that in such a case, the coefficient reads $x_n\approx \frac{1}{3}+0.11(\beta -1)$. Therefore, the magnitude of the polydispersity effect is much weaker than the previous one, in addition to the fact that it can work in both directions, depending on the value of $\beta$, which we observed to vary in the range $0.3<\beta <1.5$ (SI Appendix, Fig. S4). However, it is worth noting that a significant fraction of bubbles have a radius larger than $R_{32}$, i.e., $1\le \beta \le 1.5$, and that almost half of the nodes are bounded by one such large bubble (see SI Appendix, Fig. S4 as an example for foam with $15\%$ liquid). These findings suggest that the effect of polydispersity is only slightly positive and should only slightly increase $x_n$, i.e., the size of the wall interstice. We conclude that the loosening of the bubble packing is the main effect accounting for the measured $x_n$ values.

To extend our prediction to any liquid fraction $\varphi \le {\varphi }_{\text{rcp}}$, we turn to ref. 43, where the radius of passage of a hard sphere through the liquid channels so-called Plateau borders was determined as a function of $\varphi$ and bubble radius $R$ in a monodisperse foam. Due to the uniformity of the capillary pressure through the foam, which sets the radius of curvature of the channels, and thus their cross-section, the bubble radius at the transition $R_t$ should be proportional to this radius of passage. Following the approach proposed in ref. 43 we refer to the effective pore radius introduced by Johnson et al. (44): $\mathrm{\Lambda }\approx {(8\tilde{k}/\tilde{\sigma })}^{1/2}R$, where $\tilde{k}$ is the dimensionless liquid Darcy’s permeability through the foam structure, i.e., $k/R^2$, and $\tilde{\sigma }$ is the ratio of the electrical conductivity of the foam to that of the foaming liquid. Therefore, the expression sought for $x_n$ is

where $\xi$ is a geometrical coefficient to be determined. Note that the latter is expected to account for the loosening and polydispersity effects discussed previously. To continue, we now need expressions for $\tilde{k}$ and $\tilde{\sigma }$. Since $\mathrm{\Lambda }$ was initially proposed for solid porous media, the permeability should correspond to foam having rigid interfaces to mimic solid-like boundary conditions. As studied by Rouyer et al. (45), its expression is given by $\tilde{k}={\varphi }^2/(312{(1-2.15\varphi +1.37{\varphi }^2)}^2)$ within the range of liquid fractions $1\%\le \varphi \le 40\%$. For foams and bubbly liquids, Feitosa et al. (46) proposed an approximate analytical expression for $\tilde{\sigma }$, i.e., $\tilde{\sigma }=2\varphi (1+12\varphi )/(6+29\varphi -9{\varphi }^2)$. Using these expressions, we set $\xi =2.2$ in Eq. 4 in order to get a predicted value of $x_n$ close to the measured value 0.53 for $\varphi \approx {\varphi }_{\text{rcp}}$ (Fig. 2C). Remarkably, the agreement with our experimental data is very good over the whole range of liquid fractions, which reinforces the physical picture that $R_t$ actually corresponds to the size of the interstices formed by the jammed bubbles around the roaming bubbles. Note that in all of the above, nothing is really specific to the fact that we are looking at the wall. In bulk, typical interstices are formed by four bubbles in a tetrahedral assembly. The geometrical calculation for four bubbles in contact gives $\sqrt{\frac{3}{2}}-1\approx 0.225$, compared to $1/3$ at the wall. Therefore, we can estimate $x_n$ for bulk by using Eq. 4 with coefficient $\xi =2.2\times (0.225/0.333)\approx 1.5$. Provided this value is used, the behavior observed at the wall should be similar to the behavior observed in the bulk of the foam.

In this section, we focus on the dissolution rate of the roaming bubbles in the range $\varphi <{\varphi }^{\ast }$. We first consider the data for times longer than those that mark the intersection of the dissolution curve with $A(t)=\pi R_t^2$ (Fig. 2A and SI Appendix, S2). We follow the evolution of the radius of roaming bubbles $R(t)$ for $R(0)\lesssim R_t$. For comparison, we similarly analyze individual bubbles roaming in the bubbly liquids ($\varphi >{\varphi }^{\ast }$), from the instant they start to continuously shrink. Several examples of the curves are presented in Fig. 3A. We observe that the following function fits well all the curves (48, 49):

where the only fitted parameter ${\mathrm{\Omega }}_r$ represents the dissolution rate of the roaming bubble. Such fits were performed for all the liquid fractions and the average values of ${\mathrm{\Omega }}_r$ are presented in Fig. 3B. ${\mathrm{\Omega }}_r$ is found to depend only weakly on liquid fraction: ${\mathrm{\Omega }}_r\approx$ 1 to 2 $\mathrm{\mu }$m${}^2/s$. We also plot on Fig. 3B the growth rate ${\mathrm{\Omega }}_p$ that characterizes the coarsening of the foam in the Scaling State (Eq. 1). It appears that ${\mathrm{\Omega }}_p\gg {\mathrm{\Omega }}_r$ for $\varphi \lesssim {\varphi }_{\text{rcp}}\approx 31\%$, and ${\mathrm{\Omega }}_p\approx {\mathrm{\Omega }}_r$ for ${\varphi }_{\text{rcp}}<\varphi <{\varphi }^{\ast }$. This comparison reinforces our discussion in Section A: The size of the roaming bubbles, represented on the Left side of the distribution, varies more slowly than the average bubble size. As a result, the roaming bubbles accumulate in the interstices formed by the larger bubbles.

As the dissolution rate ${\mathrm{\Omega }}_r$ plays a crucial role in the accumulation mechanism of the roaming bubbles, we seek here to understand this value. The starting point is the comparison of our data with theory for the dissolution of isolated bubbles (48, 49), which gives the steady dissolution rate far enough from the final instant of bubble disappearance as ${\mathrm{\Omega }}_r=-d{R}^2/dt=2D_m{V}_m\left(c(R)-c_{\infty }\right)=2D_m{V}_m\text{He}{P}_0(1-\zeta )$, where the saturation parameter $\zeta =c_{\infty }/\text{He}{P}_0$ characterizes the gas saturation of the liquid environment, $c(R)$ and $c_{\infty }$ are, respectively, the gas concentrations in the liquid at the bubble surface and at infinity, $P_0$ is the gas pressure at infinity, and $\text{He}$ and $D_m$ are, respectively, the Henry solubility and the diffusion coefficient of the air molecules in the foaming solution. $V_m$ is the molar volume of the gas at the pressure $P_0$. From the measured ${\mathrm{\Omega }}_r$, we deduce an effective value for the saturation parameter: $\zeta =0.973-0.987$, which suggests that the bubbles dissolve faster than if they were isolated, and despite the presence of the large neighboring bubbles which impose at their interface a gas concentration larger than $\text{He}{P}_0$. To explain this apparent contradiction, it is important to understand that the gas transfer is controlled by the concentration gradient and not only by the concentration difference. Due to the short distances involved between the roaming bubble interface and the interfaces of the large neighboring bubbles, the concentration gradient around the roaming bubble reaches relatively high values compared to the case of the isolated bubble. Therefore, to mimic this situation, we consider the configuration illustrated in the inset of Fig. 3B, where a roaming bubble of radius $R$ is centered in a cavity of radius $R_t$ and is surrounded by a liquid shell of thickness $R_t-R$. The local concentration at the outside boundary of the shell is estimated as that at the surface of a bubble of average size $R_{32}$. From Fick’s first law, we then predict the bubble dissolution rate ${\mathrm{\Omega }}_r$ in that shell environment (see more details in SI Appendix). For the range of values of $R_{32}$ in the scaling state in our experiments and typical ratio $R/R_t$, we expect ${\mathrm{\Omega }}_r\approx$ 0.75 to 4 $\mathrm{\mu }$m${}^2$/s which provides boundaries consistent with the measured values of ${\mathrm{\Omega }}_r$ (cf. Fig. 3B).

A drawback of this shell-like model is that the roaming bubble is assumed to remain at the center of the interstice, which is not always the case. Indeed, we often noticed transient apparent contacts between the roaming bubble and either one of the bubbles delimiting the interstice or two larger bubbles forming a corner. These transient contacts can result from adhesive forces. We have indeed observed that under microgravity conditions, persistent aggregates form spontaneously in dilute bubble dispersions. In complementary ground-based experiments, we have observed a contact angle close to $3-4^o$ (34). The underlying configuration may be an adhesive contact with the formation of a liquid film that slightly flattens the bubbles or it can be a near-contact with a small separation distance so that the roaming bubble is spherical. Since it was not possible to distinguish between these two types of contact, we estimated the dissolution rate for both cases (see details of the calculation in SI Appendix). In the range of average bubble sizes $R_{32}$ of our experiments, assuming a film thickness effective for the transport of gas of the order of 40 to 60 nm (34), we found that the expected rates fall within the range of values measured for ${\mathrm{\Omega }}_r$. This remains broadly true if the bubble is in a corner, where the corresponding dissolution rate is twice larger. Therefore, whatever the configuration considered for the roaming bubble in the interstice, we find values for its dissolution rate that are compatible with our measurements, which gives robustness to the proposed mechanism based on the accumulation of long-lasting roaming bubbles in the foam interstices.

Let us analyze now the role of liquid fraction on the distribution shape. Details on the analysis are given in SI Appendix. Fig. 4 shows the normalized bubble size distributions observed in the Scaling State for each sample liquid fraction.

The PDF for $\varphi =15\%$ is the same as that of Fig. 1 in the Scaling State. It exhibits a prominent narrow peak, that we identified to the roaming bubble population in Section A, followed by a broad peak for the foam jammed bubble population. These features qualitatively persist up to $\varphi <38\%$ but the narrow peak progressively shifts toward larger $\rho$ while its height decreases. For $\varphi \ge 40\%$, PDFs exhibit a single peak, which is consistent with the fact that all bubbles should be roaming bubbles. PDFs become narrower as $\varphi$ increases and their peak height increases. This qualitative change is also captured by the abrupt variation of statistical quantities like polydispersity (cf. SI Appendix, Fig. S5). None of the existing theories predict such distributions (15). These findings indicate a cross-over between qualitatively different PDFs occurring for a liquid fraction ${\varphi }^{\ast }\approx 39\%$. This transition coincides with the observed change of growth laws Eqs. 1 and 2 and it is attributed to the onset of the formation of a foam gel due to weak attraction between bubbles as evidenced by finite contact angle at films junctions (34).

The expected jamming liquid fraction for randomly close-packed monodisperse hard spheres is ${\varphi }_{\text{rcp}}=36\%$. However, polydispersity will reduce this value since smaller bubbles can fit into the interstices between larger ones. This effect has been predicted by numerical simulations of polydisperse close packings of spherical particles with lognormal PDF, as a function of the shape parameter $\sigma$ (19, 47). In our foams, the close packing concerns the population of jammed bubbles, which are connected to each other via films. Therefore, we compare the measured shape parameter of the foam bubble distribution ${\sigma }_2$ to the predicted ones (cf. Fig. 3C). We find them to coincide within the range $\varphi =30\%$ and $\varphi =32\%$: We expect the close packing fraction ${\varphi }_{\text{rcp}}$ of our foams to lay inside the range between these 2 values.

To provide an independent result of the close packing fraction of frictionless spheres with the polydispersity observed in our samples in the Scaling State, we have performed molecular dynamics simulations. Since here we are only interested in the geometrical sphere packing problem at the jamming point where the confinement pressure and interaction forces drop to zero with increasing $\varphi$, we expect the nature of the interaction law used in the simulations to have only a minor impact. Using Hertzian interactions, in the framework of the molecular dynamics code LAMMPS (Materials and Methods), we obtained ${\varphi }_{\text{rcp}}=$31.0 $\pm 0.5\%$, in remarkable agreement with our analysis based on the work of Farr and Groot (19). Note that strictly speaking, our simulations only provide an upper bound for the optimal random close packing fraction of such polydisperse spheres, which may be obtained by more sophisticated simulation procedures described in the literature (52). However, in the context of our experiments, the truly relevant packing fraction is the one of a coarsening foam. In this case, we expect a local packing which is not exactly the most compact possible one. A jammed foam regularly undergoes rearrangements, helping it to settle into new minimal energy configurations. This implies that in between rearrangements, the packing is not always optimally close-packed. Simulations of this where we also replace the Hertzian interaction by the more realistic Morse Witten law (53, 54) are the subject of ongoing work.

The roaming bubbles represent a significant proportion of the total number of bubbles. As a result, exclusion of the roaming bubbles from the determination of the average radius leads to an up to 3-fold overestimation of coarsening rates, depending on the liquid fraction (34).

In terms of volume fraction with respect to the liquid volume, the roaming bubbles represent up to ten percent depending on $\varphi$. It can therefore be expected that their impact is important for certain properties. Roaming bubbles can also modify foam drainage, where they slow down the flow of the liquid. A study with solid spheres, located in the nodes of the liquid network of the foam, showed that such an amount of particles in the liquid could reduce the permeability of the foam by $40\%$ (55). This shows the bias of systematically ignoring their presence. Let us mention that to date, this effect has never been taken into account in permeability modeling.

The presence of the roaming bubbles can be highly detrimental for applications of foams where the microstructure is an important parameter, all the more so if they migrate and accumulate in large proportions in certain places. Moreover, for a number of these applications, the interstitial liquid is a complex fluid, possibly with yield stress properties that will prevent gravity from evacuating the roaming bubbles: One expects to find relatively high volume fractions of roaming bubbles in such systems.

By examining some recent papers on foamed construction materials, we recognized traces of the presence of such roaming bubbles. For example in ref. 56, very small pores can be seen in the bulk nodes of cement foam solidified after coarsening, as revealed by the microtomography image. There is every reason to believe that literature is full of such examples.

On the other hand, one can make use of these roaming bubbles. Note that the stakes are high in terms of producing solid foam structures with hierarchical porosity, including both macro- and microscaled pores. Such structures have recently been produced by 3D printing (57) and they were found to present enhanced energy absorption properties and enhanced mechanical resistance to cyclic loading.

S.C.-A., O.P., R.H., E.R., A.S., and D.L. designed research; N.G., M.P., A.M., A.R., S.C.-A., O.P., R.H., E.R., A.S., and D.L. performed research; R.H. contributed analytic tools; N.G., M.P., A.M., S.C.-A., O.P., and R.H. analyzed data; and N.G., M.P., S.C.-A., O.P., R.H., E.R., A.S., D.J.D., and D.L. wrote the paper.

The authors declare no competing interest.