Coarsening and mechanics in the bubble model for wet foams

TL;DR

This paper investigates coarsening dynamics and mechanical response in wet foams near the jamming transition using numerical simulations of a bubble model, revealing a diverging coarsening time and insights into elastic behavior.

Contribution

It introduces a numerical study of coarsening and mechanics in wet foams near jamming, linking coarsening rates to wetness and mechanical relaxation.

Findings

Coarsening reaches a steady state with a characteristic size distribution.

Coarsening time diverges as the foam approaches jamming.

Mechanical response is influenced by the interplay of coarsening and relaxation times.

Abstract

Aqueous foams are an important model system that displays coarsening dynamics. Coarsening in dispersions and foams is well understood in the dilute and dry limits, where the gas fraction tends to zero and one, respectively. However, foams are known to undergo a jamming transition from a fluid-like to a solid-like state at an intermediate gas fraction,$\phi_c$. Much less is known about coarsening dynamics in wet foams near jamming, and the link to mechanical response, if any, remains poorly understood. Here, we probe coarsening and mechanical response using numerical simulations of a variant of the Durian bubble model for wet foams. As in other coarsening systems we find a steady state scaling regime with an associated particle size distribution. We relate the time-rate of evolution of the coarsening process to the wetness of the foam and identify a characteristic coarsening time that diverges approaching jamming. We further probe mechanical response of the system to strain while undergoing coarsening. There are two competing time scales, namely the coarsening time and the mechanical relaxation time. We relate these to the evolution of the elastic response and the mechanical structure.

Equations

$$\frac{\langle R\rangle}{\langle R_{\rm in}\rangle}\simeq\left\{\begin{array}[]{cc}1&t\ll\tau_{c}\\ \left({t}/{\tau_{c}}\right)^{\alpha}&t\gg\tau_{c}\end{array}\right.$$

$${\text{F}}_{ij}=F_0\left(\frac{R_i+R_j-\mid {\text{r}}_i-{\text{r}}_j\mid }{R_i+R_j}\right)\frac{{\text{r}}_i-{\text{r}}_j}{\mid {\text{r}}_i-{\text{r}}_j\mid },$$

$${\text{v}}_i=\frac{1}{{\mu }_0}\sum _{i\ne j}^N{\text{F}}_{ij},$$

$$\frac{d{\text{x}}_i}{dt}=\sum _{i\ne j}^N⟨R_{\mathrm{in}}⟩\left(\frac{R_i+R_j-\mid {\text{r}}_i-{\text{r}}_j\mid }{R_i+R_j}\right)\frac{{\text{r}}_i-{\text{r}}_j}{\mid {\text{r}}_i-{\text{r}}_j\mid },$$

$$\frac{d{V}_i}{dt}=K{A}_{ij}\left(\frac{1}{R_j}-\frac{1}{R_i}\right),$$

$$\text{K}\text{u}(t^{\ast })+\text{B}\dot{\text{u}}(t^{\ast })=\sigma (t^{\ast })V{\text{e}}_{\gamma },$$

$${\sigma }_{\mathrm{el}}$$

$${\sigma }_{\mathrm{visc}}$$

$$\frac{{\tau }_r}{{\tau }_c}\sim \mathcal{T}\left(\frac{t^{\ast }}{{\tau }_c}\right),$$

Full text

Coarsening and mechanics in the bubble model for wet foams

Kseniia Khakalo,∗a Karsten Baumgarten,b Brian P. Tigheb and Antti Puistoa

Aqueous foams are an important model system that displays coarsening dynamics. Coarsening in dispersions and foams is well understood in the dilute and dry limits, where the gas fraction tends to zero and one, respectively. However, foams are known to undergo a jamming transition from a fluid-like to a solid-like state at an intermediate gas fraction,$\phi_{c}$. Much less is known about coarsening dynamics in wet foams near jamming, and the link to mechanical response, if any, remains poorly understood. Here, we probe coarsening and mechanical response using numerical simulations of a variant of the Durian bubble model for wet foams 1. As in other coarsening systems we find a steady state scaling regime with an associated particle size distribution. We relate the time-rate of evolution of the coarsening process to the wetness of the foam and identify a characteristic coarsening time that diverges approaching jamming. We further probe mechanical response of the system to strain while undergoing coarsening. There are two competing time scales, namely the coarsening time and the mechanical relaxation time. We relate these to the evolution of the elastic response and the mechanical structure.

††footnotetext: a* Aalto University, School Science, Laboratory of Applied Physics B.O.B 11100, FI-00076 AALTO, Finland; E-mail: Kseniia.Khakalo@aalto.fi*††footnotetext: b* Delft University of Technology, Process & Energy Laboratory, Leeghwaterstraat 39, 2628 CB Delft, The Netherlands*

1 Introduction

Foams are composed of repulsively interacting gas bubbles dispersed in a liquid phase. Based on gas fraction they are categorized as “wet” or “dry” 2. In wet foams, the liquid concentration is higher and the bubbles mostly retain their spherical shape. In dry foams, only thin films of liquid separate the bubbles due to the limited amount of liquid available. Therefore, the bubbles in dry foams appear in polyhedral shapes. The mechanical properties and rheology of a foam also depend on its gas fraction. At gas fractions well below the critical value $\phi_{c}$ ($\approx 0.84$ in 2D and 0.64 in 3D), their mechanical response is mainly determined by the background fluid 2. At gas fractions above $\phi_{c}$, they form amorphous, jammed assemblies. Beyond this limit, their mechanical response can be expected to depend on the bubble properties, such as their interaction potentials and size distribution.

Foams are not thermodynamically stable, and the bubble size distribution evolves due to destabilization of the foam. The three principal mechanisms governing foam destabilization are drainage, coarsening, and coalescence 2. Drainage removes liquid between the bubbles via gravity or evaporation increasing the gas concentration. Coarsening occurs when gas diffuses from smaller bubbles to larger ones, thanks to the difference in their internal pressure 3. Coalescence, where bubbles join when thin films rupture between them, primarily takes place in dry foams.

Self-similar scaling arguments provide valuable insight into the prolonged coarsening observed in foams 2, 4, 5, 6. In the scaling state, the growth of bubbles is expected to reach an asymptotic limit with the average radius following the scaling law

[TABLE]

where $\langle R_{\rm in}\rangle$ is the average bubble radius at time $t=0$ and $\tau_{c}$ is a time scale characteristic of the coarsening dynamics. Experiments have indeed confirmed the asymptotic limit 4, 5, 6, 3. The scaling exponent is $\alpha=1/2$ in the dry limit ($\phi\rightarrow 1$) and $\alpha=1/3$ in the limit of bubbly fluids ($\phi\rightarrow 0$). There is some numerical evidence that $\alpha$ interpolates between these values for intermediate $\phi$ 7, but to date there have been no studies that systematically probe coarsening in the vicinity of the jamming point.

Here, we model coarsening using a variant of the Durian bubble model 8. This is a soft sphere model, describing the evolution of the spatial configuration of overlapping spheres resembling bubbles, droplets or soft particles in a dispersed system. Eventhough the model is rather simple, it has been proven very useful in the research of jamming in foams and emulsions 9, 10. To further advance its capabilities to study coarsening, we have implemented inter-bubble gas diffusion as an additional degree of freedom to our set of dynamical equations. This is done in the spirit of Gardiner et al. 1, who studied coarsening in the bubble model at gas fractions far above the jamming limit. At this range the model successfully reproduced the asymptotic limit, Eq. 1, with the appropriate scaling of the average radius.

The present work represents the first numerical study of coarsening in the bubble model near jamming. We present several main results. First, the value of the critical volume $\phi_{c}$ is altered by the coarsening dynamics. We find that aged samples are particularly efficient at filling voids, and $\phi_{c}$ increases to the unusually large value of $0.87$ in 2D. Next, near jamming the scaling state of Eq. (1) holds, with an exponent $\alpha\approx 0.45$ and a coarsening time $\tau_{c}$ that diverges on approach to $\phi_{c}$. Finally, coarsening dramatically influences mechanical response near jamming. We characterize the storage modulus and the viscous relaxation time, both of which scale with sample age and the distance to $\phi_{c}$.

2 Bubble model with coarsening

We model foams using Durian’s bubble model 8, 11. The bubble model describes foams at the bubble level as packings of randomly distributed soft spheres (or disks in 2D). These spheres interact via a harmonic pair potential proportional to their overlap

[TABLE]

as illustrated in Fig. 1a. The model assumes the fully overdamped limit: the bubbles are massless, and the force due to bubble overlap ${\textbf{F}}_{i}$ must at all times be compensated by the drag force ${\textbf{F}}^{d}_{i}=-\mu_{0}({\textbf{v}}_{i}-\langle{\textbf{v}}_{j}\rangle)$. Here $\langle{\textbf{v}}_{j}\rangle$ is the velocity of the background fluid often computed as the average velocity of the neighboring bubbles, and $\mu_{0}$ is the viscosity of the fluid. Since we impose no external deformation to our system, we set $\langle{\textbf{v}}_{j}\rangle$ = 0. Then, the bubbles follow a quasistatic equation of motion with

[TABLE]

which can be integrated in a molecular dynamics fashion to obtain the evolution of the foam. From considerations of dimensionality assume that $F_{0}=\sigma_{0}\langle R_{\rm in}\rangle$. In this case, the dynamics of the system is determined by ratio $\mu_{0}/\sigma_{0}$. Therefore, the time scale in the simulation can be chosen so that Eqs. 2 and 3 become

[TABLE]

where $t$ is dimensionless time, scaled with $\mu_{0}/\sigma_{0}$. The main merits of the model are that it is sufficiently simple, while it still allows to easily vary the foam properties such as polydispersity, volume fraction, and dimensionality.

In coarsening, gas diffuses from smaller bubbles to larger ones driven by the difference of their respective Laplace pressures. To take this into account in the bubble dynamics model, Gardiner, Dlugogorski, and Jameson (henceforth GDJ) proposed a scheme where, in addition to their elastic and viscous interactions, the bubbles are allowed to exchange gas 1. The gas exchange rate of two contacting bubbles is proportional to their interaction area $A_{ij}$ and the difference in their Laplace pressures. Each bubble’s volume $V_{i}$ changes according to

[TABLE]

where $K$ is the diffusion parameter, encapsulating the properties of the liquid film and the bubbles, such as the effective permeability, surface tension, and temperature. We integrate this set of differential equations using a second order adaptive step size predictor-corrector scheme with error tolerances set to $1\cdot 10^{-6}$.

The simulation procedure is as follows. The simulations begin by first randomly distributing 3000 bubbles in a periodic rectangle (2D) at the initial volume fraction of $\phi=0.45$. For each bubble an initial radius is assigned according to a Gaussian distribution with the mean $\langle R_{\rm in}\rangle=0.006$ and variance $21\%$. To reach the target volume fraction $\phi_{0}$ we then compress the structure by rescaling the dimensions of the simulation cell at a constant velocity. After the compression, we equilibrate the system by allowing the bubble positions to relax until the energy changes less than 0.0001% for 1000 iterations. This gives us an initial structure, such as the one shown in Fig. 1b. Finally, we turn on the coarsening and run the simulations until the number of bubbles is smaller than a cut-off, which we take to be 300, and the bubble size distribution has reached the scaling state (Fig. 1c). In all the plots involving time, we have used the time scale $t^{*}=\frac{K}{\langle R_{\rm in}\rangle^{2}}t$.

3 Scaling state

During an induction period at the beginning of the coarsening, the probability density function (PDF) of bubble sizes $R/\langle R(t^{*})\rangle$ broadens (Fig. 2), while the average bubble size $\langle R(t^{*})\rangle$ remains approximately constant – see Fig. 3a. After the induction time the system reaches a scale-independent regime where ${\rm PDF}(R/\langle R(t^{*})\rangle)$ ceases to evolve with time (see Fig. 2). At each volume fraction, the PDF in the scaling state exhibits a significantly broader shape, along with a high number of small bubbles of size around $R\approx 0.1R_{max}$. This most likely relates to the way the system optimizes the occupied space by filling the voids between the big bubbles with fitting small ones. The small bubbles have no overlap with any of their neighbors (are so-called rattlers). As a result, having no neighbors, they will not experience any gas exchange until their surrounding bubbles change their configuration. Furthermore, the rattler radius geometrically allowed in a foam increases with the average bubble size.

The average bubble radius in the scaling state follows a power law $\langle R\rangle\sim{t^{*}}^{\alpha}$, as observed in Fig. 3a for a range of volume fractions. While not all volume fractions reach the scaling state, among those that do we find an exponent $\alpha\approx 0.45$. This is strikingly close to the value of $1/2$ required by von Neumann’s law in perfectly dry foams ($\phi=1$) – despite the fact that von Neumann’s law does not hold exactly in the bubble model. Due to dimensional considerations and the conservation of total volume, one also anticipates related scaling relations for the mean contact “area” (length in $D=2$ dimensions) between bubbles, $\langle A\rangle\sim t^{\alpha}$, as well as the total number of bubbles, $N\sim{t^{*}}^{-D\alpha}$. These are verified in Fig. 3b and c.

It is apparent from Fig. 3a that the onset of the scaling state is determined by a volume fraction-dependent time scale $\tau_{c}(\phi)$. The volume fraction $\phi\approx 0.87$ appears to be a marginal case; lower values of $\phi$ show no evidence of scaling within the simulation runtime $T=100$, while larger values show clear power law scaling at long times. In order to test for connections to the jamming transition, in Fig. 3d we plot the mean coordination number for varying $\phi$ as a function of time. In determining $z$, we assume that the system’s evolution is quasistatic, so that “rattlers” can be meaningfully identified and removed. Recall that jamming occurs at the critical coordination number $z_{c}=2D=4$ in two dimensions, in accord with a constraint counting argument that dates to Maxwell. One clearly sees that packings for $\phi\geq 0.88$ satisfy $z>z_{c}$ for the entire runtime, while $\phi=0.87$ dips below 4 approximately half way through the run – the critical value $\phi_{c}$ has increased due to coarsening. As all initial conditions are jammed, the increase in the jamming volume fraction, and the corresponding drop in the average coordination number, relate to the change in the bubble size distribution under coarsening. When the bubble size distribution reaches the scaling state the bubbles fill space more optimally, by reducing the bubble overlap via the gas diffusion. In order to estimate $\phi_{c}$ more accurately, we have probed the interval $0.87\leq\phi\leq 0.88$ in steps of $0.001$ for longer runs to $T=1000$. At these long times the system reaches sizes $N\sim{\cal O}(100)$ (recall that $N\sim 1/{t^{*}}^{D\alpha}$) and the coordination number experiences large fluctuations, which restricts our ability to determine a precise jamming volume fraction $\phi_{c}$. We find that for volume fractions $\phi\leq 0.872$ the mean coordination unambiguously drops below $z_{c}$ within the runtime, while for $\phi\geq 0.877$ it clearly remains above $z_{c}$. For $0.873\leq\phi\leq 0.876$, the coordination number fluctuates on either side of $z_{c}$, suggesting that the system’s evolution passes through both jammed and unjammed configurations.

The time scale $\tau_{c}(\phi)$ where the scaling state is reached can be estimated from the evolution of the average bubble radius. In Fig. 4 we estimate $\tau_{c}$ as the time where $\langle R\rangle/R_{\rm in}$ is equal to 2, i.e. the average bubble has increased its radius by $100\%$. The same can be done by computing time where initial number of bubbles decreased by a factor of 2. The data in both cases are reasonably fit by a power law $\tau_{c}\propto 1/(\phi-\phi_{c})^{\alpha}$, with fitting parameters $\phi_{c}=0.873$ and $\alpha\approx 0.45$. Deviations at small values of $\Delta\phi\equiv\phi-\phi_{c}$ are clearly present; we expect they are associated with the fluctuations between jammed and unjammed states identified above.

4 Mechanics

In the previous Section we identified signatures of the jamming transition in the scaling dynamics of coarsening foams. We now seek to correlate these effects with features of the mechanical response.

Bubbles store energy when distorted and dissipate energy when sliding past each other, which gives rise to viscoelastic response.12, 13 Foam viscoelasticity can be probed both experimentally and numerically by measuring the complex shear modulus $G^{*}(\omega)=G^{\prime}(\omega)+\imath G^{\prime\prime}(\omega)$, defined as the complex ratio of shear stress and strain amplitude under harmonic forcing at angular frequency $\omega$.14 The real and imaginary parts of $G^{*}$ are known as the storage and loss moduli, respectively. For reference, we recall that simple viscoelastic solids (frequently referred to as Kelvin-Voigt solids) have a constant storage modulus $G_{0}$ and a linear loss modulus $\eta_{0}\omega$; their ratio $G_{0}/\eta_{0}$ selects a characteristic relaxation time. Deviations from the Kelvin-Voigt form are associated with a spectrum of relaxation times (rather than a single time scale), but the time scale $\tau_{r}$ where the loss and storage moduli cross, i.e. $G^{\prime}(1/\tau_{r})=G^{\prime\prime}(1/\tau_{r})$ is still an important reference point. It indicates a crossover from predominantly solid-like response at low frequencies to predominantly liquid-like response at high frequencies.

While the complex shear modulus is meant to describe steady state dynamics, the structure of coarsening foams in the scaling state continuously evolves (“ages”).15, 16 Aging violates time translational invariance, which is assumed in the usual definition of the complex shear modulus.17 Nevertheless, oscillatory rheology is often used to characterize soft materials, with a focus on frequencies that are fast compared to the evolution of the structure.15, 16 Here we determine the storage and loss moduli via a numerical “experiment” that disentangles viscous relaxation from coarsening-induced rearrangements. States are sampled at varying times $t^{*}$ from a coarsening simulation at fixed $\phi$. Coarsening dynamics are then turned off (all particle sizes are held fixed) and the complex shear modulus is measured according to the method described below. In this way it is possible to obtain moduli over the full range of frequencies $0\leq\omega\leq\infty$; however, we focus on time scales that are short compared to the system’s age.

In order to measure the complex shear modulus, we employ a linearization scheme introduced in Ref. 18.

Given a particular configuration of bubbles, its collective response to shear is described by the $DN+1$-component vector ${\textbf{u}}=(\vec{u}_{1},\vec{u}_{2},\ldots\vec{u}_{D},\gamma)^{T}$, where $\vec{u}_{i}$ is the displacement of bubble $i$ from its reference position, and $\gamma$ is the shear strain experienced by the unit cell. The response to a shear stress $\sigma$ is given by the solution to the first order differential equation

[TABLE]

where ${\rm\textbf{e}}_{\gamma}$ is a unit vector along the strain coordinate. The stiffness matrix K and damping matrix B describe the elastic and viscous forces on the particles, respectively. K consists of second derivatives of the elastic potential energy with respect to the particle and strain degrees of freedom; B is similarly defined in terms of the Rayleigh dissipation function. Details are available in Ref. 18. Linearization is strictly valid only when deformation amplitudes are infinitesimal; nevertheless, numerical studies indicate that moduli calculated in this way remain accurate over a finite strain interval 19, 20, 21. By Fourier-transforming Eq. 6 and solving for the complex shear strain in response to a sinusoidally oscillating shear stress with frequency $\omega$, one can determine complex shear modulus.

In Fig. 5, we plot the storage and loss moduli (solid and dashed curves, respectively) for a system at $\phi=0.878$, close to but above the jamming transition. The storage modulus displays a low frequency plateau, indicating that the sampled configurations are jammed solids. (We stress again that the linearization scheme employed here “turns off” coarsening, hence any softening at asymptotically low frequencies due to coarsening-induced rearrangements will not be captured.) At high frequencies there is a second plateau in $G^{\prime}$, associated with affine deformations 18, 21. There is a gradual crossover between these two plateaus. In previous work it was shown that this crossover occupies a widening window in frequency as $\phi\rightarrow\phi_{c}^{+}$, such that $G^{\prime}\sim G^{\prime\prime}\sim\omega^{1/2}$.18 The loss modulus $G^{\prime\prime}$ is comparatively simple; it is nearly linear over the entire frequency range, consistent with a Kelvin-Voigt solid. Experimental measurements of the loss modulus in foams often show a plateau at low frequency.15, 16 This feature is absent from our data, and from prior studies of Durian’s bubble model without coarsening.18 We suggest that the plateau results from physics that is not incorporated in the bubble model, e.g. due to thin film flow.

The storage modulus shows a clear dependence on the sample age, with an overall downward shift with increasing age. Their intersection point defines the mechanical relaxation time $\tau_{r}$, which we measure at varying volume fraction and age – see Fig. 6a. The relaxation time clearly depends on both the system age $t^{*}$ and the distance to jamming; it grows larger with increasing age and decreasing distance to jamming $\Delta\phi$. The dependence of $\tau_{r}$ on $t^{*}$ and $\Delta\phi$ can be rationalized in two steps, beginning with its growth with age.

The age dependence of $\tau_{r}$ is controlled by the scaling of $G^{\prime}$ and $G^{\prime\prime}$ with $t^{*}$. These can be anticipated on dimensional grounds, by which one expects there to be characteristic elastic and viscous stress scales

[TABLE]

Here $N$, $R$, $F_{\rm el}$, and $F_{\rm visc}$ are typical values of the particle number, particle radii, and elastic and viscous forces, respectively; $L^{D}$ is the volume of the unit cell, which is constant. The elastic force law scales with the dimensionless overlap (c.f. Eq. 2) and should therefore be independent of time in the scaling state; hence the time-dependence of the typical elastic stress scales as $\sigma_{\rm el}\sim N(t^{*})R(t^{*})\sim{t^{*}}^{-\alpha}$. By contrast, the typical viscous force is set by the bubble velocity $V\propto R(t^{*})\,\omega$. Hence $\sigma_{\rm visc}\sim{t^{*}}^{0}$, consistent with observations. Turning back to $\tau_{r}$, we note that the relaxation time in systems without coarsening is insensitive to the distance to jamming; hence the dependence here is likely to be inherited from the coarsening dynamics. And indeed, a simple balancing of the viscous and elastic stress scales would suggest a frequency $\omega_{r}\equiv 1/\tau_{r}$ that depends on age as $\omega_{r}\sim R\sim{t^{*}}^{\alpha}$, consistent with the data in Fig. 6a.

In order to understand the dependence of $\tau_{r}$ on $\Delta\phi$, we postulate that the coarsening time $\tau_{c}$ sets the natural units for both relaxation time and the age of the system. This is expressed most naturally in the form of a scaling ansatz,

[TABLE]

for some function ${\cal T}(x)$. Indeed, in Fig. 6b we obtain good data collapse when plotting $\tau_{r}\,\Delta\phi^{\alpha}$ versus $t^{*}\,\Delta\phi^{\alpha}$. Treating $\alpha$ as a free parameter, the best collapse is found for $\alpha=0.55$, close to the value $0.45$ determined independently above. (Good data collapse for $\phi<0.95$ can also be obtained using $0.45$.) As expected, ${\cal T}\sim x^{\alpha}$ for large values of $x$, when the system’s age is large compared to the coarsening time $\tau_{c}$. It follows that, in the scaling state, the mechanical relaxation time obeys $\tau_{r}\sim{t^{*}}^{\alpha}/\Delta\phi^{\alpha(1-\alpha)}$. Deviations from a slope of $\alpha$ occur when the age is smaller than the coarsening time.

5 Conclusions

We have presented the results of numerical simulations of the coarsening of foams close to, and slightly above, the jamming volume fraction. For this purpose, we implemented the Durian bubble model, with extensions to incorporate gas diffusion between the bubbles.

Our main observation are: (i) The model captures the expected ${t^{*}}^{\alpha}$ scaling, with $\alpha\approx 0.5$, of the average radius well above the jamming limit. (ii) The critical volume fraction where jamming occurs is shifted to $\phi_{c}\approx 0.873$, significantly higher than that found, e.g., in typical bidisperse mixtures 22. (iii) Approaching $\phi_{c}$, the characteristic coarsening time $\tau_{c}$ diverges as $1/\Delta\phi^{\alpha}$. (iv) The foam’s mechanical response is radically influenced by the coarsening; the mechanical relaxation time $\tau_{r}$, where the oscillatory rheology shows a crossover from liquid- to solid-like behavior, shows scaling both with the system’s age and $\Delta\phi$. Points (ii-iv) are all related to the change of the bubble distribution due to coarsening, which increases the jamming volume fraction $\phi_{c}$, unjamming the system partially or completely at times. The final point establishes a clear connection between mechanical relaxation and the coarsening dynamics.

Finally, our results suggest directions for future work. In order to more deeply understand the enhanced packing efficiency of the scaling state, it would be necessary to model the form of the bubble size distribution directly. There may be fruitful connections to systems undergoing rupture and/or Apollonian packings.23, 24 An additional question concerns the dominant mechanism of gas exchange between bubbles. In the present simulations, fluxes are between pairs of particles in contact. Obviously this slows down as the gas fraction decreases. In bubbly liquids, by contrast, gas exchange is mediated by the fluid. A third possibility is that, close to the jamming volume fraction, there is a non-negligible flux through the Plateau borders (i.e. through the packing’s fluid-filled voids, rather than through the increasingly smaller thin film interfaces durianAPS, 25. This would yield additional terms in Eq. 5. Additionally, in technological applications foams often undergo shear flow during coarsening; how are these two forms of driving coupled? In addition, industrial foams are often formed of thixotropic complex fluids, such as (nano)particulate suspensions. Coarsening dynamics changes due to the non-linear dynamics of the suspending liquid, ultimately stopping completely26, 25 – how can coarsening and mechanics be characterized and modeled in such cases? Many of these issues and questions can potentially be addressed with straightforward extensions of the present model.

6 Acknowledgements

A.P. and K.K. are grateful to the Academy of Finland for financial support through project No. 278367. K.B. and B.P.T. acknowledge funding from the Netherlands Organization for Scientific Research (NWO).