Effects of liquid fraction and contact angle on structure and coarsening in two-dimensional foams

Abstract

Aqueous foams coarsen with time due to gas diffusion through the liquid between the bubbles. The mean bubble size grows, and small bubbles vanish. However, coarsening is little understood for foams with an intermediate liquid content, particularly in the presence of surfactant-induced attractive forces between the bubbles, measured by the interface contact angle where thin films meet the bulk liquid. Rigorous bubble growth laws have yet to be developed, and the evolution of bulk foam properties is unclear. We present a quasistatic numerical model for coarsening in two-dimensional wet foams, focusing on growth laws and related bubble properties. The deformation of bubble interfaces is modelled using a finite-element approach, and the gas flow through both films and Plateau borders is approximated. We give results for disordered two-dimensional wet foams with $256$ to $1024$ bubbles, at liquid fractions from $2\,\%$ to $25\,\%$, beyond the zero-contact-angle unjamming transition, and with contact angles up to $10^\circ$. Simple analytical models for the bubble pressures, film lengths and coarsening growth rates are developed to aid interpretation. If the contact angle is non-zero, we find that a prediction of the coarsening rate approaches a non-zero value as the liquid fraction is increased. We also find that an individual bubble's effective number of neighbours determines whether it grows or shrinks to a good approximation.

1. Introduction

Aqueous foams are packings of gas bubbles in liquid, as illustrated in figure 1. They have elasticity from the surface tension of the interfaces, and plasticity due to bubble rearrangements (Weaire & Hutzler 1999). A foam's properties hence differ substantially from those of its components (Cantat et al. 2013), and it is a model rheological material whose microstructure is accessible in experiments (Denkov et al. 2012; Stewart & Hilgenfeldt 2023). Foams also have a multitude of applications, from foods and drinks (Weaire & Hutzler 1999) to soil treatment (Géraud et al. 2016) and fire suppression (Martin 2012b).

Figure 1. Simulated foam structure for (a) liquid fraction $\phi = 10\,\%$ and no bubble attraction, and (b) $\phi = 3\,\%$ and contact angle $\theta = 10^\circ$, taken from our simulations. The components of the foam are labelled, and the interface discretisation (schematically), local interface separation $h$ at a mesh vertex, and $\theta$ are shown. The surface tension in the films and Plateau borders is $\gamma _{f}$ and $\gamma _\infty$, respectively, with $\gamma _{f} = \gamma _\infty \cos \theta$ (Langevin 2020, p. 88).

However, foams are unstable, ageing due to coarsening in addition to film breakage (Chae & Tabor 1997) – we consider the latter no further, noting that it can be suppressed in experiments (Roth, Jones & Durian 2013; Pasquet et al. 2023b). Coarsening arises from the diffusion of dissolved gas through the liquid, primarily through the thin films between bubbles; the gas is transported from small to large bubbles due to the higher pressures of the former (Schimming & Durian 2017). Thus, the mean bubble volume increases as the small bubbles vanish (Lambert et al. 2010). This process may be detrimental in fire-suppression applications, for example, due to the increased mixing of air and fuel vapour within larger bubbles (Martin 2012b), and can hasten the perishing of foods (Martin 2012a). While the dynamics of a confined foam under coarsening differ from bulk foams, the pressure difference required to initiate flow through a porous medium, during soil treatment for example, is also affected by coarsening (Jones, Getrouw & Vincent-Bonnieu 2018). Additionally, foam coarsening is an accessible model for grain growth in crystalline solids (Smith 1952).

In the dry limit of small liquid fraction $\phi$, which is the ratio of liquid area to total foam area (Cantat et al. 2013), the coarsening growth rate of a bubble in a two-dimensional foam is determined only by its number of neighbours, through von Neumann's law (von Neumann 1952) which we discuss in § 3.3.1. Experiments and simulations (Glazier & Weaire 1992; Stavans 1993) both show that the dry foam approaches a scaling state, in which the bubble size and neighbour-number distributions scale uniformly with time $t$ (Mullins 1986; Lambert et al. 2010). For a given bubble, let $R$ be the radius of a circle with the same area (Princen 1988), i.e. its effective radius. It may be shown that the mean effective radius increases as $\langle R \rangle \sim t^{1/2}$ in such a state, since diffusion through thin films between bubbles is the dominant mode of gas transfer in a dry foam (Mullins 1986). Similarly, the scaling state of a three-dimensional dry foam is also well established, with the same behaviour of $\langle R \rangle$ (Thomas, de Almeida & Graner 2006; Lambert et al. 2010), where $R$ is now the radius of a sphere with the same volume. However, the individual bubble growth rates depend upon their precise geometries in this case (MacPherson & Srolovitz 2007; Cantat et al. 2013).

But real foams always have non-zero $\phi$, and large values may be encountered in applications such as fire suppression (Laundess et al. 2011) and the fabrication of solid foams (Cantat et al. 2013; Galvani et al. 2023). For coarsening in the wet limit $\phi \to 1$, i.e. Ostwald ripening (Stavans 1993), a scaling state is obtained with $\langle R \rangle \sim t^{1/3}$ in two and three dimensions (Cantat et al. 2013, p. 77). The exponent differs from the dry case because diffusion through the bulk liquid is now the dominant gas transfer mechanism (Mullins 1986). However, less is known about coarsening at moderate $\phi$. Experiments are difficult to control because drainage of the foam's liquid under gravity occurs on shorter time scales than coarsening (Weaire & Hutzler 1999; Born et al. 2021). Nevertheless, experiments have been performed using diamagnetic levitation of foams (Isert, Maret & Aegerter 2013), as well as in microgravity on the International Space Station (ISS) (Born et al. 2021; Galvani et al. 2023; Pasquet et al. 2023a,b). A narrow transition was found between the limiting growth exponents $1/2$ and $1/3$ of $\langle R \rangle$, over an interval of around $15\,\%$ or less in $\phi$ (i.e. $\Delta \phi = 0.15$) near the unjamming transition $\phi = \phi _{c}$ at which the bubbles (in a foam without bubble attraction) lose contact. This occurs at $\phi _{c} \approx 16\,\%$ in two dimensions, and $\phi _{c} \approx 36\,\%$ in three (Cantat et al. 2013, p. 195).

The theory of coarsening at moderate $\phi$ also remains limited. In two dimensions, bubble growth laws have been developed for foams to which the decoration theorem applies (Bolton & Weaire 1991; Roth et al. 2013; Schimming & Durian 2017). This theorem states that a two-dimensional equilibrium wet foam containing only Plateau borders which meet exactly three films, typically holding when $\phi \lesssim 3\,\%$ (Jing et al. 2021), is also an equilibrium dry foam when the Plateau borders are omitted (Bolton & Weaire 1991). These growth laws have been compared with experiments in Hele-Shaw cells (Roth et al. 2013; Chieco & Durian 2021), accounting for the Plateau borders along the bounding plates. An interpolation between the known growth laws for zero and large $\phi$ was proposed, and found, when averaged, to agree with simulations at intermediate $\phi$ (Fortuna et al. 2012). The gas flow rates between adjacent circular or spherical bubbles, present near the unjamming transition $\phi _{c}$ (in the absence of bubble attraction), have also been derived (Schimming & Durian 2017). The latter work has been extended to three-dimensional bubbles with films, although the growth rates of individual bubbles are not yet predicted (Durian 2023). We are not aware of any general and rigorous growth laws for $0 < \phi < \phi _{c}$, in either two dimensions (except when the decoration theorem applies) or three. Nor are we aware of a fully developed theory that predicts quantitatively the transition in growth exponents observed experimentally, although progress has been made by Durian (2023).

Several simulations of coarsening in wet foams have been performed (Bolton & Weaire 1991; Gardiner, Dlugogorski & Jameson 2000; Fortuna et al. 2012; Thomas et al. 2015; Khakalo et al. 2018). Evidence for scaling states has been found over a range of $\phi$, in two and three dimensions. However, the Potts model simulations (Fortuna et al. 2012; Thomas et al. 2015) predict a different form for the transition of the growth exponent of $\langle R \rangle$ with $\phi$, compared with experiments (Isert et al. 2013; Pasquet et al. 2023b) and simulations (Khakalo et al. 2018) using the bubble model (Durian 1995). As $\phi$ increases from zero, an immediate decrease of the exponent from $1/2$ is observed in the Potts model simulations, resulting in a discrepancy for $\phi < \phi _{c}$. While the bubble-model simulations qualitatively reproduce the experimental transition (albeit in two dimensions), the relative rate of diffusion across films and Plateau borders is not predicted.

Most numerical studies have used models which are suited to simulating large numbers of bubbles – $3000$ to $10\,000$ for the bubble model (Khakalo et al. 2018), and approximately $10^5$ for the Potts model (Fortuna et al. 2012; Thomas et al. 2015) – but which do not include accurate bubble deformation. A bubble's growth rate depends on the portion of its surface in contact with other bubbles, along with its pressure and those of its neighbours (Roth et al. 2013), each of which is determined by the bubble geometry. The pressure is related to interface curvature via the Young–Laplace law (Weaire & Hutzler 1999). To our knowledge, only Bolton & Weaire (1991), Benzi et al. (2015) and Pelusi, Sbragaglia & Benzi (2019) have performed numerical coarsening studies that accurately model the bubble shapes in wet foams. The first study is limited to small $\phi$, while the latter two, which used a lattice Boltzmann approach, were primarily concerned with bubble rearrangements during coarsening.

Furthermore, the effect of the foam's contact angle $\theta$ on coarsening has not been widely investigated. As illustrated in figure 1(b), this is the angle between the tangents of the film and Plateau-border interfaces where they meet (Ivanov & Toshev 1975; Denkov, Petsev & Danov 1995), which arises from an imbalance in their surface tensions (Cantat et al. 2013, p. 49). We neglect for the moment the transition region between these interfaces (Kralchevsky & Ivanov 1985b). The contact angle is determined by the surfactant, via the disjoining pressure (Langevin 2020, p. 88), and increases with the degree of attraction between bubbles (Princen 1983). A characteristic property of foams with $\theta > 0$ is flocculation, whereby bubbles cluster due to their attraction (Princen 1983; Cox et al. 2018). While $\theta$ may be negligible in typical foams (Höhler, Seknagi & Kraynik 2021), a contact angle of $\theta \approx 4^\circ$ is thought to have affected the results of the ISS coarsening experiments, by delaying the transition in growth exponents to $\phi > \phi _{c}$ (Pasquet et al. 2023b). Furthermore, experiments have produced films with $\theta > 10^\circ$ (Princen 1968; Seknagi 2022). Contact angles also occur in emulsions (Bibette et al. 1993), to which we expect our results also apply, due to their similar structure to foams (Weaire & Hutzler 1999). Prior work has been done to characterise foam structure at $\theta > 0$ (Cox et al. 2018; Feng et al. 2021; Jing et al. 2021; Jing & Feng 2023), along with foam rheology (Menon, Govindarajan & Tewari 2016). We are not aware of numerical studies investigating the effects of $\theta > 0$ on coarsening, although non-zero $\theta$ has been used for technical reasons (Fortuna et al. 2012; Thomas et al. 2015).

In this article, we present a quasistatic numerical model of coarsening in two-dimensional wet foams, which accurately models the bubble geometries, and allows $\theta > 0$. Two-dimensional foams are widely studied as more tractable models of the real three-dimensional systems (Kähärä, Tallinen & Timonen 2014; Cox et al. 2018; Khakalo et al. 2018), and can be approximately realised in Hele-Shaw cells (Smith 1952; Roth et al. 2013). Our approach to the foam structure is adapted from the models of Kähärä et al. (2014) and Boromand et al. (2018), the latter having been widely applied, including to foams and emulsions (Boromand et al. 2019; Golovkova et al. 2021). These methods are suited to relatively small systems, with approximately $1700$ bubbles or fewer. We implement our simulations in Kenneth Brakke's Surface Evolver software (Brakke 1992,  2013).

We vary $\theta$ by altering the attractive component of the disjoining pressure. Droplet attraction in emulsions has previously (Golovkova et al. 2021) been implemented in the model of Boromand et al. (2018) using a similar approach.

Our coarsening model is inspired by that applied analytically by Marchalot et al. (2008) and Schimming & Durian (2017), and approximates the gas flow through both the liquid films and the Plateau borders.

Our guiding assumption is that results from accurately modelled small systems, for which the scaling state is likely inaccessible (Thomas et al. 2006), may give insight into coarsening in macroscopic foams, including by refining the necessary approximations used in simulations of large systems.

We describe our numerical methods in § 2. Section 3 gives our results from simulating two-dimensional disordered foams at various $\phi$, for several values of $\theta$ comparable to those observed experimentally (Princen 1968), and we conclude in § 4. Instead of simulating the time evolution of foams under coarsening, we consider the instantaneous bubble growth rates, and related properties, in foams at structural equilibrium, due to our relatively small system sizes of $256$ and $1024$ bubbles. Analytical approximations are developed to assist with interpretation, and we give a model for the bubble film lengths in Appendix A. Methods we use for measuring simulated foam and bubble properties are described in Appendix B, and a test of foam equilibration is discussed in Appendix C.

2. Numerical methods

We begin by summarising our methods for equilibrating foams, before defining our coarsening model, and our approach for generating the initial foam structure.

2.1. Structural model

2.1.1. Discretisation

The foam structure is modelled using the approach of Kähärä et al. (2014) and Boromand et al. (2018). The bubbles are bounded by a closed interface, with arbitrary shape, of mesh vertices connected by straight edges. This is illustrated schematically in figure 1(b). Hence, all liquid–gas interfaces are explicitly included – the bubbles are disconnected, and their rearrangements can occur without adjustments in the discretisation. We emphasise that vertices and edges refer here (and throughout) to elements of the discretisation, rather than to the infinitesimal Plateau borders and films, respectively, in a dry foam.

The bubble areas are fixed, unlike in the models of Kähärä et al. (2014) and Boromand et al. (2018), taking the foam's gas to be effectively incompressible (Cantat et al. 2013, p. 27). The liquid is identified with the region outside the bubbles. Since coarsening occurs on substantially longer time scales than structural equilibration (Thomas et al. 2015), we use a quasistatic approach like Boromand et al. (2018). This contrasts with Kähärä et al. (2014), who applied their model to flowing foams. We also neglect gravity, and therefore choose to measure all pressures relative to the uniform pressure of the liquid (thus taken to be zero).

Another approach (Bolton & Weaire 1992; Jing et al. 2015; Cox et al. 2018) is instead to assume that the films separating bubbles have zero thickness, and to treat them separately from the Plateau-border interfaces. While this likely has advantages in numerical efficiency, since the interfaces are all circular arcs (Bolton & Weaire 1992), handling bubble rearrangements in this model is complex for wet foams (Cox et al. 2018), particularly in three dimensions (Weaire & Hutzler 1999, p. 80). Furthermore, a small amount of bubble adhesion is often needed for numerical stability (Jing et al. 2015; Cox et al. 2018), and the criteria for adding or removing a film between slightly contacting bubbles could induce artefacts near the unjamming transition.

However, a disadvantage of the model we use is that, for numerical convergence, the film thickness $h_0$ must not be too small, though the minimum thickness can be decreased by refining the mesh. For the degree of refinement we use (stated in § 3), the minimum usable film thickness is approximately $10^{-2} \langle R \rangle$, where $\langle R \rangle$ is the mean bubble radius. This is substantially larger than in real foams, for which $h_0 \lesssim 10^{-4} \langle R \rangle$ (Cantat et al. 2013, pp. 18–20). The film thickness is set by the interactions between different interfaces, which we now describe.

2.1.2. Disjoining pressure

Similarly to Kähärä et al. (2014), we implement bubble interactions through a disjoining pressure acting between nearby liquid–gas interfaces. This ensures that bubbles do not overlap at equilibrium. In order to allow $\theta$ to be varied, we select a disjoining pressure inspired by the Derjaguin–Landau–Vervey–Overbeck (DLVO) theory for parallel flat interfaces with separation $h$ (Cantat et al. 2013, p. 94), which we also use between curved interfaces. The interface separation in the latter case (see figure 1) is discussed further in § 2.1.4.

The DLVO model includes an electrostatic repulsion term dominant at smaller $h$, which has the form ${\rm e}^{-\kappa h}$ for constant $\kappa > 0$. A van der Waals attraction proportional to $1/h^3$ is present which dominates at larger $h$ (Langevin 2020, pp. 89–93). Since we expect transient bubble overlaps ($h \leq 0$) in the simulations, during structural relaxation only, we alter the form of the second interaction to $h {\rm e}^{-\kappa h}$ to improve stability. We also select $\kappa = 1/h_0$, where $h_0$ is the equilibrium film thickness – a choice which is again made for stability, as discussed at the end of this subsection. Hence, we use the disjoining pressure

(2.1)\begin{equation} \varPi_{D}(h) = A (1 - \alpha h) \, {\rm e}^{{-}h / h_0} , \end{equation}

where $A$, $\alpha$ and $h_0$ are positive constants. Positive disjoining pressure corresponds to interface repulsion. Our aim in selecting this form is a qualitative model for $\varPi _{D}$, which allows arbitrary $\theta$ to be set via the constant $\alpha$ (i.e. by varying the relative strengths of the attractive and repulsive components). While (2.1) agrees qualitatively with DLVO theory for $h$ of moderate value and above, it does not reproduce the dominance of van der Waals attraction over electrostatic repulsion in the latter model at very small $h$, nor the consequent local maximum of $\varPi _{D}(h)$. This omission corresponds to the assumption that only common black films, for which $\varPi _{D}$ is below this maximum, are present in the foams we simulate (Cantat et al. 2013, pp. 97–98). We note that the real functional form of $\varPi _{D}$ depends upon the surfactant (Langevin 2020).

We are not aware of prior studies using this form for $\varPi _{D}$. Kähärä et al. (2014) and Boromand et al. (2018) use a repulsive harmonic interaction. Golovkova et al. (2021) incorporate attractive interactions, but use a piecewise-linear form for $\varPi _{D}$ without relating its attractive components to a contact angle. We performed test simulations using a similar piecewise-linear $\varPi _{D}$, varying the range of attractive interactions, to confirm that our results are not sensitive to the disjoining pressure model.

The constants $A$ and $\alpha$ in (2.1) are set as follows. We enforce that the equilibrium film thickness $h_0$, which is a parameter of the simulations, is attained by a flat film separating two bubbles whose pressures equal the foam's capillary pressure $\varPi _{C}$. This is the area-weighted mean bubble pressure – if the $i$th bubble has area $A_i$ and pressure $p_i$, then (Höhler et al. 2021)

(2.2)\begin{equation} \varPi_{C} = \frac{\sum_i A_i p_i}{\sum_i A_i} , \end{equation}

recalling that our liquid pressure is zero. Hence, $\varPi _{D}(h_0) = \varPi _{C}$ by balancing the pressures on the flat film interfaces (Toshev & Ivanov 1975). We note that $\varPi _{C}$ is not known before the simulations are run, so this condition is applied iteratively, as described in § 2.3. Films separating bubbles with different pressures will have different equilibrium thicknesses (Princen 1988), as discussed later.

The remaining degree of freedom is used to obtain the desired $\theta$. Let $\gamma (h)$ be the surface tension (in dimensions of force, since our model is two-dimensional) of a liquid–gas interface in a flat film of thickness $h$, and let $\gamma _\infty$ be the tension of an isolated interface outside a film. Then (Langevin 2020, p. 88)

(2.3)$$\begin{gather} \gamma(h) = \gamma_\infty + \frac{1}{2} \int_h^\infty \varPi_{D}(s) \, {\rm d} s , \end{gather}$$

(2.4)$$\begin{gather}\cos[\theta(h)] = \gamma(h) / \gamma_\infty . \end{gather}$$

The latter result arises from balancing the surface tension forces at one end of a flat film (Kralchevsky & Ivanov 1985a), with $\theta$ defined as in figure 2. There is strictly no discontinuity in the tangents of simulated liquid–gas interfaces which corresponds to $\theta$, due to the transition regions between films and Plateau borders (Kralchevsky & Ivanov 1985b). An interesting consequence of (2.3) and (2.4) is that $\varPi _{D}$ must have an attractive component even for $\theta (h) = 0$ (Ivanov & Toshev 1975), so that $\gamma (h) = \gamma _\infty$. This is because $\varPi _{D}(h) > 0$ in an equilibrium film (Cantat et al. 2013, p. 96). The contact angle is undefined if $\gamma (h) > \gamma _\infty$.

Figure 2. Definition of the contact angle $\theta$ accounting for the transition region between films and Plateau borders. A circular arc describing the Plateau-border interface outside the transition region is depicted, along with a line similarly describing the film interface. Their angle of intersection gives $\theta$ (Kralchevsky & Ivanov 1985b; Denkov et al. 1995).

The contact angle $\theta _{m}$ measured by approximating the simulated bubble interfaces by a collection of circular arcs (Kralchevsky & Ivanov 1985b; Denkov et al. 1995), as described in Appendix B.3, differs from the value of $\theta$ expected from (2.4) with $h = h_0$. When $\theta (h_0) = 0$, we find a typical $\theta _{m} \approx 3^\circ$ for $\phi \approx \phi _{c}$, which manifests the attractive component of $\varPi _{D}$ in this case, while the discrepancy decreases as $\theta (h_0)$ increases. As discussed further in § 3.3.2, we find that $h > h_0$ for films in foams near the unjamming transition. Hence, from (2.3), and the short-range repulsion in $\varPi _{D}$, (2.4) then predicts a larger contact angle. We parametrise our foams with respect to $\theta (h_0)$, denoted $\theta$ henceforth, which we recall determines the interface tension in a flat film between bubbles with the capillary pressure.

By substituting (2.1) into (2.3) and (2.4) for $h = h_0$, and applying the above condition $\varPi _{D}(h_0) = \varPi _{C}$, our disjoining pressure is

(2.5)\begin{align} \varPi_{D}(h) &= \frac{2 \gamma_\infty}{h_0} \left[\left(1 - \cos\theta + \frac{\varPi_{C} h_0}{\gamma_\infty}\right)\right. \nonumber\\ &\quad \left. - \left(1 - \cos\theta + \frac{\varPi_{C} h_0}{2 \gamma_\infty}\right) \frac{h}{h_0}\right] {\rm e}^{1 - h / h_0} . \end{align}

In order to model foams with no bubble attraction, we also implement a repulsive disjoining pressure

(2.6)\begin{equation} \varPi_{D}(h) = \dfrac{\varPi_{C}}{e - 1} \begin{cases} {\rm e}^{2 - h / h_0} - 1 , & \text{for }h \leq 2 h_0 ;\\ 0, & \text{otherwise} . \end{cases} \end{equation}

This is obtained by omitting the term in $\alpha$ from (2.1), and including a cutoff such that $\varPi _{D}(h) = 0$ for $h > 2 h_0$ (so bubble neighbours can be reliably calculated). A constant is subtracted for $h \leq 2 h_0$ so $\varPi _{D}(h)$ is continuous, and the condition $\varPi _{D}(h_0) = \varPi _{C}$ is then applied. We note that $\gamma (h) > \gamma _\infty$ and $\theta$ is undefined (inevitable for a repulsive $\varPi _{D}$), by (2.3) and (2.4). Equations (2.5) and (2.6) are plotted in figure 3.

Figure 3. The form of disjoining pressure $\varPi _{D}$ we use, given by (2.5) and (2.6), versus film thickness $h$ (relative to its equilibrium value $h_0$). For $\theta > 0$, the capillary pressure $\varPi _{C}$ does not scale out. The larger value is representative of a simulated foam with $\phi = 2\,\%$, and the smaller for a flocculated foam with larger $\phi$ (see § 3.1). Each curve satisfies $\varPi _{D}(h_0) = \varPi _{C}$, as stated in the text.

Equation (2.3) is not exact for curved interfaces, which are ubiquitous in foams, and a realistic disjoining pressure would depend upon the interface curvature (Denkov et al. 1995). However, we neglect these effects as a simplifying assumption. This is the Derjaguin approximation, justified when the radii of curvature of the interfaces are large compared with their separation $h$ (Denkov et al. 1995). The latter assumption does not hold in general for our simulations, although it is appropriate for the weakly curved interfaces of thin films. The mean ratio of film-interface radius of curvature to film thickness is measured to be approximately $1500$ in our simulations (at $\phi = 2\,\%$ and $10\,\%$ for repulsive $\varPi _{D}$, and $\phi = 2\,\%$ and $25\,\%$ at $\theta = 10^\circ$), but the ratio varies from around $9$ up to approximately $10^5$ for individual films. This variation is due to that in curvature rather than in film thickness. We measure both film properties using the methods described in Appendix C.

We now explain further the requirement for larger $h_0$ in our simulations than in real foams, along with our selection of $\kappa = 1/h_0$ for (2.1). For a given mesh refinement, it appears necessary for convergence that $\varPi _{D}$ not vary too rapidly with $h$ near $h_0$. This may be related to a discretisation-induced difference between $h$ on the two sides of a film, discussed later in Appendix C. In the context of the simulations, $h_0$ could in principle be set arbitrarily small, but, to avoid unphysical interface overlaps at equilibrium, $|\varPi _{D}'(h_0)|$ would need to be correspondingly large (recalling that variations in bubble pressure correspond to variations in the value of $\varPi _{D}$ in equilibrium films). Furthermore, even allowing for our larger $h_0$, convergence requires that $|\varPi _{D}'(h_0)|$ is much smaller than found in real foams (Bergeron & Radke 1992), giving rise to larger variations in film thickness between bubbles of different pressure. However, we note that some variations are expected in real foams (Princen 1988). The derivative $|\varPi _{D}'(h_0)|$ is set by the choice of $\kappa$, defined above – we take $\kappa = 1/h_0$ in (2.1), which is large enough that bubbles do not overlap at equilibrium, but sufficiently small for the simulations to converge at the mesh refinement we use (stated in § 3). We discuss the resulting variations in the film thickness of bubbles in § 3.3.2.

The disjoining pressure $\varPi _{D}(h)$ is applied to each mesh vertex on the liquid–gas interfaces, unlike Kähärä et al. (2014) and Boromand et al. (2018) who apply it to the edges in their foam models. We set $h$ equal to the shortest distance to another interface (the local interface separation), as shown in figure 1(b).

2.1.3. Vertex neighbours

To determine the local interface separation $h$ at a mesh vertex $v_i$, neighbour searches are first performed to find the vertex closest to $v_i$ which lies on a different interface, illustrated in figure 4(a).

Figure 4. Schematic of (a) the local neighbour search, and (b) the piecewise-linear interface extrapolation. The outward unit normals to the edges adjoining the neighbouring vertex in (b) are $\boldsymbol {n}_1$ and $\boldsymbol {n}_2$, the displacement of the vertex from its neighbour is $\boldsymbol {d}$, and $h$ is the extrapolated shortest distance from the vertex to an opposing interface.

In order to improve the efficiency of these searches, compared with a brute-force approach, we use the fact that all vertices lie on particular bubble interfaces, and their adjoining vertices thereon do not change. Our nearest-neighbour algorithm is as follows.

1. (i) Cover each bubble with a circle whose centre is the centroid of its vertices, and whose radius is the minimum required to cover each of these vertices.

2. (ii) For each bubble, determine the bubbles whose circles overlap its own, illustrated in figure 5. These are the neighbouring bubbles, and are found by calculating the distance between each circle centre. This search is sufficiently fast since our systems contain relatively few bubbles (no more than $1024$). For these overlap checks, the circle radii are scaled by $125\,\%$ (an arbitrary value, but sufficiently large). If a scaling were not performed, then interface overlaps would occur during bubble rearrangements, as some bubbles would intersect one another before registering as neighbours.

3. (iii) For each vertex $v_i$ on each bubble, calculate the closest vertex on each neighbouring bubble. If the latter bubble was a neighbour during the last pass of the algorithm, then the previously closest vertex provides an initial guess. Otherwise, we take the vertex closest to the centroid of the bubble on which $v_i$ lies (since this vertex need only be found once for all vertices on the bubble of $v_i$). The closest vertex is then obtained through a local search along the neighbouring bubble's interface, starting at the guessed vertex, and, at each step, moving to an adjoining vertex if it is closer to $v_i$. This is illustrated in figure 4(a). The process is only guaranteed to give a vertex at a local minimum of distance, but works well in practice due to the smooth bubble geometries (see figure 1) and the choice of initial guess.

4. (iv) For each vertex $v_i$ on each bubble, compare the distance to the closest vertex on each neighbouring bubble. The closest of these is then the nearest neighbour of vertex $v_i$.

Figure 5. Schematic of the covering circles used to determine neighbouring bubbles for the purpose of finding vertex neighbours. The circle for a particular bubble is shown, including its fringe used to improve numerical stability, along with the circles of its neighbouring bubbles (i.e. those that overlap the first circle).

We are not aware of descriptions of the nearest-neighbour approaches used in prior simulations of this type, although we note that Okuda & Hiraiwa (2023a) use octrees in comparable three-dimensional simulations of biological cells (Okuda & Hiraiwa 2023b). The efficiency of these algorithms appears to be of importance in scaling such foam simulations to larger systems and higher mesh refinement.

2.1.4. Interface extrapolation

Having determined the nearest neighbour of each vertex, we then calculate the local interface separation $h$. We could set $h$ equal to the distance from the vertex to its neighbour, but this would result in an unrealistic roughness in the interfaces – this approach has previously been used to model static friction in other materials (Boromand et al. 2018).

Issues of interface extrapolation or interpolation also arise in other fluid dynamics simulations. For example, Bazhlekov, Anderson & Meijer (2004) use a spherical interpolation between vertices to obtain the local interface separation in three-dimensional boundary integral simulations.

Instead, we use a piecewise-linear interface extrapolation, such that the local interface separation is the minimum distance to the half-infinite extensions of the two edges adjoining the neighbour vertex, as illustrated in figure 4(b). Let $\boldsymbol {n}_1$ and $\boldsymbol {n}_2$ be outward unit normals to the edges, and let $\boldsymbol {d}$ be the displacement of the considered vertex from its neighbour. Then we define

(2.7)\begin{equation} h = \begin{cases} \max\left(\boldsymbol{d} \boldsymbol{\cdot} \boldsymbol{n}_1, \boldsymbol{d} \boldsymbol{\cdot} \boldsymbol{n}_2 \right), & \text{if the neighbouring interface is convex} ,\\ \min\left(\boldsymbol{d} \boldsymbol{\cdot} \boldsymbol{n}_1, \boldsymbol{d} \boldsymbol{\cdot} \boldsymbol{n}_2\right), & \text{if the neighbouring interface is concave} . \end{cases} \end{equation}

The two cases are shown in figure 6. Equation (2.7) ensures that $h < 0$ if the vertex overlaps the neighbour interface during structural relaxation, so the disjoining pressure acts to oppose the overlap. However, (2.7) differs from the exact piecewise-linear interface extrapolation when the considered vertex is closest to its neighbouring vertex, rather than another point on the latter's adjoining edges. In the absence of overlap, this can occur only for convex neighbouring interfaces, when the considered vertex lies in the shaded sector of figure 6(a), where the effect on a distance contour is shown: $h$ is underestimated by (2.7) here.

Figure 6. Illustrations of (2.7) for (a) convex, and (b) concave, neighbouring interfaces, with notation from figure 4(b). The edges adjoining the neighbour vertex are shown, as are the infinite extensions thereof (dotted lines). The considered vertex may in principle lie anywhere relative to the neighbour, including in the shaded sector of (a) where $h$ would be the distance to the neighbour vertex itself in an exact implementation of the piecewise-linear extrapolation. The dotted and dashed lines in this sector show a given $h$ contour using the exact extrapolation and (2.7), respectively.

The difference from the exact extrapolation is expected to cause only small errors in $h$, due to the smooth interfaces of equilibrium foams and our degree of mesh refinement. Hence, the angle of the sector in figure 6(a) is anticipated to be small in practice. For the $1024$-bubble foams discussed in § 3 (for repulsive $\varPi _{D}$ and $\theta = 10^\circ$, and at liquid fraction $2\,\%$ and $25\,\%$), approximately $4\,\%$ of vertices are affected by the difference at liquid fraction $\phi = 2\,\%$, rising to approximately $22\,\%$ at $\phi = 25\,\%$. For these affected vertices, the resulting root-mean-square relative error in $h$ is below around $0.2\,\%$, and the maximum magnitude of the relative error is approximately $2.3\,\%$. These errors are judged to be small enough to justify using the simplified (2.7).

We use the above extrapolation due to its closeness to the geometry of the discretisation. Other approaches might allow improved convergence at lower mesh refinement, and may avoid the imbalance in $h$ between the interfaces of a curved film noted in Appendix C. Kähärä et al. (2014) and Boromand et al. (2018) bypass the issue of interface extrapolation by explicitly calculating the shortest distance between edges. However, we determine neighbouring vertices, rather than neighbouring edges, under the assumption of greater efficiency, due to the simpler distance calculations involved.

2.1.5. Implementation

We implement our simulations using the Surface Evolver, developed by Kenneth Brakke (Brakke 1992,  2013). This software is frequently applied in the study of foams (Kraynik, Reinelt & van Swol 2003; Jing et al. 2015; Höhler et al. 2021), usually under the assumption that the liquid films have zero thickness.

The Surface Evolver may be extended using its scripting language. We have implemented finite film thickness in this manner, via a disjoining pressure and neighbour searches. Our scripts are not compiled, so it is likely that our simulations could be made faster by implementing the routines in the software's public source code. We also do not take advantage of any significant parallel processing.

We note that pairwise repulsion between vertices already exists in the Surface Evolver as ‘knot’ energies (Brakke 2013). However, it is not straightforward to adapt these to our purposes, due to our desire for interface extrapolation and interactions only between nearest neighbours.

As is usual in the Surface Evolver, local energy minimisation is used to obtain the equilibrium foam structures. Let $\varGamma _{F}$ be the union of all liquid–gas interfaces in the simulated foam. Then the foam's total energy (recalling that its liquid and gas are treated as incompressible) is given by

(2.8)\begin{equation} E = \int_{\varGamma_{F}} \gamma(h) \, {\rm d} l , \end{equation}

where $\gamma (h)$ is obtained from (2.3), and $h$ is the local interface separation from (2.7) (which varies around $\varGamma _{F}$). This quantity is minimised using conjugate gradient iterations, as implemented in the Surface Evolver (Brakke 1992,  2013).

Let $\boldsymbol {x}_i$ be the position of the $i$th vertex, $\boldsymbol {x}_{i_1}$ and $\boldsymbol {x}_{i_2}$ those of its adjoining vertices on the same interface, $\boldsymbol {x}_n$ that of its nearest neighbour, and $\boldsymbol {x}_{n_1}$ and $\boldsymbol {x}_{n_2}$ those of its neighbour's adjoining vertices. We will consider foams with periodic boundary conditions, so these positions are those of the nearest copies of the vertices to $\boldsymbol {x}_i$. Define $l_i = (|\boldsymbol {x}_{i_1} - \boldsymbol {x}_i| + |\boldsymbol {x}_{i_2} - \boldsymbol {x}_i|) / 2$ as the length of interface associated with the $i\text {th}$ vertex (half the length of its adjoining edges), and $h_i(\boldsymbol {x}_i, \boldsymbol {x}_n, \boldsymbol {x}_{n_1}, \boldsymbol {x}_{n_2})$ as the corresponding local interface separation (calculated as in § 2.1.4). Also, let $L_i = |\boldsymbol {x}_{i_1} - \boldsymbol {x}_i|$ be the length of the $i\text {th}$ edge, between the $i\text {th}$ vertex and its adjoining vertex at $\boldsymbol {x}_{i_1}$. If there are $N$ vertices, the energy $E$ is hence expressed in the simulations as

(2.9)\begin{equation} E(\{\boldsymbol{x}_i\}) = \sum_{i=1}^N \left(\gamma_\infty L_i(\boldsymbol{x}_i, \boldsymbol{x}_{i_1}) + \frac{l_i(\boldsymbol{x}_i, \boldsymbol{x}_{i_1}, \boldsymbol{x}_{i_2})}{2} \int_{h_i(\boldsymbol{x}_i, \boldsymbol{x}_n, \boldsymbol{x}_{n_1}, \boldsymbol{x}_{n_2})}^\infty \varPi_{D}(s) \, {\rm d} s \right) , \end{equation}

where the integrals are determined explicitly from (2.5) or (2.6). If $x_k$ is the $k\text {th}$ coordinate in $\{\boldsymbol {x}_i\}$, then we approximate the gradient $\boldsymbol {\nabla } E(\{x_k\})$ by (in component form)

(2.10)\begin{equation} \frac{\partial E}{\partial x_k} \approx \sum_{i=1}^N \left(\gamma_\infty \frac{\partial L_i}{\partial x_k} - l_i \varPi_{D}(\bar{h}_i) \frac{\partial \bar{h}_i}{\partial x_k} \right) , \end{equation}

where $\bar {h}_i(\boldsymbol {x}_i) \equiv h_i(\boldsymbol {x}_i, \boldsymbol {x}_n, \boldsymbol {x}_{n_1}, \boldsymbol {x}_{n_2})$ for fixed neighbour positions (i.e. these are taken as parameters rather than variables). This replacement is made because our implementation does not allow derivatives of $h_i$ with respect to the neighbour coordinates. We compensate by doubling the second summed term in (2.9) when calculating $\boldsymbol {\nabla } E$, as done in (2.10). Using Newton's third law, this is equivalent to including such derivatives under the approximation that each vertex is the nearest neighbour of its own nearest neighbour, with each vertex having the same $l_i$ as its neighbour (these lengths are equal to a tolerance of $20\,\%$ for around half of vertices), and that the corresponding distances $h_i$ are those between the vertices. We find that, of the vertices close enough to their neighbour that $\varPi _{D}$ is not negligible (i.e. $h_i < 8 h_0$, by figure 3), approximately $3/4$ are the neighbour of their neighbour in our simulations. As an approximation, we have also neglected the contributions to $\boldsymbol {\nabla } E$ from derivatives of $l_i$ in (2.10), interpreted as corrections from $\varPi _{D}$ to the surface tension in the first summed term of (2.10). These corrections are expected to be small for our $\theta \lesssim 10^\circ$ (see § 3) by (2.4), and to increase with $\theta$.

The bubble area constraints mentioned in § 2.1.1 are set using Lagrange multipliers via the Surface Evolver's usual procedures (Brakke 1992), which gives the gas pressure of each bubble as the value of the corresponding multiplier.

Equation (2.10), including the mentioned constraint terms and an added restriction that vertices can move only in the direction normal to their interface, is used to evolve the vertex coordinates $x_k$ at each conjugate gradient iteration, using the Surface Evolver's internal routines. Because the bubbles are expected to rearrange slowly, we perform the time-consuming neighbour searches every $20$ iterations. However, the local interface separations $\bar {h}_i$ are recalculated before each iteration, using the current coordinates of the vertex neighbours. After the same interval of $20$ iterations, the edges on each individual bubble interface are set to a uniform length, with a tolerance of $5\,\%$, by using the Surface Evolver's vertex averaging routines (edge elasticity was instead implemented by Kähärä et al. (2014)). The energy $E$, obtained from (2.9), is checked every $20$ iterations to see whether it has changed by less than a threshold (specified in § 3) since the last check. If so, then the foam is considered equilibrated.

Transiently, and usually during the early stages of structural relaxation, the assumptions of the interface extrapolation may fail. For example, spikes in the discretisation of a bubble's interface (i.e. vertices with a large displacement from those adjoining them on the interface) may result in a large interface overlap ($h \ll 0$) being falsely registered for a vertex. This would cause the vertex to experience an extremely large repulsive force, resulting in further spikes which may render the mesh unusable. To avoid this, we detect instances of large overlap (according to the interface extrapolation), and set the coordinates of the overlapping vertex and its neighbour to the mean coordinates of the vertices adjoining them on their interfaces. We also detect spikes using the angle turned by successive edges on an interface, suppressing them by the same method. Small interface overlaps, which do not result in the overlapping vertex swapping neighbours, can be resolved without such interventions (see §§ 2.1.2 and 2.1.4).

Our use of the Surface Evolver's relaxation routines to implement vertex interactions (we term this method 1 for structural relaxation) have resulted in the following caveat: the vertex neighbour properties (and hence the $\bar {h}_i$) are not recalculated during a given conjugate gradient iteration, interfering with the automatic selection of an ideal step size (Brakke 1992). This slows relaxation, and means the energy is no longer guaranteed to decrease. In practice, we find that this method is still considerably faster than gradient descent with a fixed step size (method 2, to which the caveat does not apply), and that systematic energy increases are rare for our system parameters (with their occurrence also being reduced at higher mesh refinement). We tried implementing our vertex interactions in the Surface Evolver source code to allow an exact use of the conjugate gradient algorithm (method 3), but we found the result unreliable due to the step size decreasing to zero prior to convergence being clearly established. This may be due to an incompatibility between the system energy and its calculated gradient, noting that our vertex interactions are asymmetric, but the issue remains unresolved. Therefore, we use method 1 in the results reported here.

We have compared foams equilibrated by each of methods 1 to 3 (for $256$ bubbles, $\phi = 2\,\%$, and $\theta = 10^\circ$), and find them consistent, with method 1 resulting in the lowest final energy (this energy differs by less than $0.002\,\%$ among the methods). For method 3, we included the contributions to $\boldsymbol {\nabla } E$ from derivatives of $l_i$ in (2.9) (recalling that we expect these contributions to increase with $\theta$), and so this comparison also tests our mentioned neglect of these contributions in method 1, which we use hereafter. We provide a further check of the equilibration of our simulated foams in Appendix C.

2.2. Coarsening model

Having given our methods for structural relaxation, we now describe our coarsening model.

Coarsening occurs due to the transport of dissolved gas through the foam's liquid (Cantat et al. 2013). Let $c(\boldsymbol {r}, t)$ be the gas concentration at position $\boldsymbol {r}$ in the liquid, and at time $t$. Following Schimming & Durian (2017), we assume that almost all the gas transfer between bubbles occurs when the gas concentration is close to equilibrium. Hence, $c$ satisfies Laplace's equation

(2.11)\begin{equation} \nabla^2 c = 0 . \end{equation}

Boundary conditions are given by Henry's law $c = H p$ on the interface of each bubble, where $p$ is the bubble's pressure and $H$ is Henry's constant, related to the solubility (Cantat et al. 2013, p. 109). Only differences in $c$ are relevant, so pressures relative to that of the liquid may be used.

Let $\delta F$ be the gas flow rate (in dimensions of area per time) across an element $\delta l$ of a bubble's interface, due to a nearby bubble. Using the approach of Marchalot et al. (2008) and Schimming & Durian (2017), we approximate this rate as that between two infinite, parallel, straight interfaces, separated by the local interface separation $h$ (the same as used to determine $\varPi _{D}$). Let the pressure difference between the two bubbles be $\Delta p$. The solution of Laplace's equation is linear between such interfaces, so (Cantat et al. 2013, p. 109)

(2.12)\begin{equation} \delta F \approx D H \, \delta l \Delta p / h , \end{equation}

where $D$ is a diffusion coefficient. The gas flows towards the bubble with lower pressure. We apply this approximation to each vertex, taking $\delta l$ to be half the sum of the lengths of the two adjoining edges (i.e. $l_i$ for the $i\text {th}$ vertex, from § 2.1.5). The relevant quantities are illustrated in figure 7.

Figure 7. Illustration of the quantities used in (2.12). With these definitions, $\delta F$ is a contribution to the loss of gas from the bubble with pressure $p_1$.

By summing $\delta F$ over each vertex on a bubble's interface, we obtain an approximation to its area growth rate $\dot {A}$, which accounts for gas transfer through both the thin films and Plateau borders.

Both contributions were previously included in lattice Boltzmann simulations (Benzi et al. 2015; Pelusi et al. 2019), although these studies were focused on the motion of bubbles during coarsening, rather than their growth rates. The contributions were also present in the phase-field simulations of Fan et al. (2002). However, their focus was not on foam-like systems. Otherwise, to our knowledge, prior simulations of wet foams with accurate bubble geometries have only included gas flow through the films (Bolton & Weaire 1991) – this is the border blocking assumption (Roth et al. 2013). Simulations using simplified bubble geometries have implemented approximate contributions of both types (Fortuna et al. 2012; Thomas et al. 2015; Khakalo et al. 2018), although more work is needed to justify their forms for the Plateau border contributions, which were originally derived for $\phi \gg \phi _{c}$.

Schimming & Durian (2017) compared their approximate gas flow rates through Plateau borders, obtained by a similar approach to the above, with those given by solving Laplace's equation numerically. They obtained close agreement for thin films. However, our simulations use relatively thick films to ensure convergence, and $\delta l$ is finite (around one hundredth of the bubble perimeter), whereas $\delta l \to 0$ in their analytical results. Therefore, the error in our bubble growth rates is not fully characterised, though we expect the approximation to be effective in thin films, since it is based on the solution to Laplace's equation between parallel straight interfaces. We also note that each element of interface can exchange gas with only one neighbouring bubble, which may induce further errors in the gas flow through Plateau borders, where neighbouring bubbles meet (R. Höhler, personal communication). However, comparing the predictions of our coarsening model with existing growth laws in § 3.3, we find fair agreement for the gas flow through Plateau borders at small and large liquid fractions, thus supporting our approximate approach.

We note that our coarsening model does not conserve the total area of gas in the foams. Let $\sigma (\dot {A})$ denote the standard deviation of the bubble growth rate distribution. Since the mean growth rate satisfies $\langle \dot {A} \rangle = 0$ when the gas is conserved (Cantat et al. 2013, p. 105), a dimensionless measure of non-conservation is $|\langle \dot {A} \rangle | / \sigma (\dot {A})$. For the systems described later in § 3, this relative error takes a value below $1\,\%$ for $\theta = 10^\circ$. For repulsive $\varPi _{D}$, the error increases from less than $1\,\%$ at $\phi = 2\,\%$ to approximately $2\,\%$ when $\phi \approx \phi _{c}$. At $\phi = 24.5\,\%$, for which the bubbles are out of contact, the error is around $8\,\%$. For small $\phi$, we believe that this error comes partly from a difference in $h$ between the interfaces of a thin film, caused by the discretisation, which is discussed in Appendix C.

2.3. Generation of disordered foams

We now describe our process for generating the initial disordered foam structure.

As is usual (Cox et al. 2018; Jing et al. 2021), we begin with a Voronoi tessellation of the plane, with a domain containing the desired number of bubbles subject to periodic boundary conditions. These tessellations are obtained using the vor2fe software by Kenneth Brakke (Brakke 1986). The resulting dry foam is then equilibrated using standard techniques (Brakke 1992), with one straight edge per liquid film for efficiency. Relaxing the dry foam before setting the liquid fraction $\phi$ (Jing et al. 2015) is useful for computational efficiency.

To obtain a more representative bubble area distribution for coarsening foams than that of a Voronoi tiling, we sample areas $A$ from the distribution fitted by Roth et al. (2013) to experimental data for quasi-two-dimensional foams (i.e. foams in Hele-Shaw cells). This is a compressed exponential distribution, defined (up to normalisation) by (Roth et al. 2013)

(2.13)\begin{equation} \rho(A / \langle A \rangle) = \left(\frac{A}{\langle A \rangle}\right)^{\beta_1 - 1} \exp\left[-\beta_2 \left(\frac{A}{\langle A \rangle}\right)^{\beta_1}\right] , \end{equation}

where $\langle A \rangle$ is the mean area, $\beta _1 = 1.21$, and $\beta _2 = 0.926$. The sampled areas are randomly assigned to the bubbles in the Voronoi foam (Jing et al. 2015), without regard to their initial areas. The bubble areas are changed gradually, in tenths of the difference between the old and new values, with equilibration at each step. We note that the data of Roth et al. (2013) are for foams with small $\phi$, and may differ from the truly two-dimensional case we model.

The sampled areas are more polydisperse than those obtained from a Voronoi tessellation. Let $R_{2 1} \equiv \langle R^2 \rangle / \langle R \rangle$ (Cantat et al. 2013, p. 251). Then the bubble distributions have a polydispersity $\mathcal {P} = R_{2 1} / \sqrt {\langle R^2 \rangle } - 1$ (Cox et al. 2018) of $0.093$ for the sampled areas and $0.036$ for the Voronoi tessellations. The former value was evaluated numerically from (2.13), and the latter calculated for a single 16 384-cell tessellation (though $\mathcal {P}$ fluctuates for the relatively small foams we study in § 3). We have not observed qualitative differences in the foam or bubble properties considered in § 3 when the Voronoi-tessellation areas are used. The same relationships between growth rate and radius, for example, are observed at the lower polydispersity, but only over the narrower interval of realised bubble radii. Hence, even if the sampled area distribution is unrealistic for larger $\phi$, its high polydispersity remains useful for characterising the relationships for a wider range of bubble radii.

Next, the system is annealed following the approach of Kraynik et al. (2003), adapted to two dimensions, so that the distribution of bubble topology better approximates that in a real foam (Kraynik et al. 2003). By applying linear transformations to the foam's periodic domain, along with the mesh vertices therein, a sequence of extensional step strains with size $\{6/5, 5/6, 5/6, 6/5\}$ (Kraynik et al. 2003) is applied along both axes in turn, with relaxation after every strain. Further transformations are then applied to the periodic domain and the vertices, by applying cycles of simple and extensional shear (Weaire & Hutzler 1999), to relax the deviatoric stress, measured according to Appendix B.1, to zero (within a tolerance $10^{-5}\,\gamma _\infty / \langle R \rangle$). We do this by applying small shears of the respective type to estimate the modulus, before using a shear which would zero the strain given this modulus and the present stress, were the system linear. Stress relaxation of the initial structure was performed in the dry foam coarsening simulations of Herdtle & Aref (1992). We have not investigated the effects of varying this preparation process.

The liquid fraction $\phi$ is then set, by duplicating the mesh edges and vertices associated with each bubble (since these are shared by neighbours in the dry model), and uniformly expanding the simulation domain to $1/(1 - \phi )$ of its original area. The same transformation is applied to the bubble centroids, while the bubble shapes and sizes are maintained. Hence, $\phi$ is set by effectively adding liquid to the system, keeping the gas area fixed (Cox et al. 2018).

We next subdivide all edges on the bubble interfaces $m$ times in order to obtain the desired mesh refinement. Some of the initial edges are very short, so these are skipped to improve convergence. The longer edges are further subdivided to compensate. The mean number of edges (and vertices) per bubble is then $n_{e} = 6 \times 2^m$, since $6$ is the mean number of bubble neighbours (and hence mesh edges) in the initial dry foam (Weaire & Hutzler 1999, p. 29). After initial structural relaxation, vertices are added or deleted to ensure small and large bubbles have acceptable numbers of vertices. The criteria are that we refine once all edges longer than $1.05/n_{e}$ of the perimeter of a circle with the mean bubble area (i.e. not much longer than edges in a typical bubble with its initial refinement), and we delete half of each bubble's edges which are shorter than $4.76/n_{e}$ (approximately $1/20$ for $m = 4$, as used in § 3) of the perimeter of a circle with the smallest bubble's area.

The foam's structure is relaxed using the approach of § 2.1. We recall, from § 2.1.2, that the strength of $\varPi _{D}$ (determining the equilibrium film thickness $h_0$) is set via the foam's capillary pressure $\varPi _{C}$, which is not known beforehand. Hence, we initially relax the foam using a rough estimate of $\varPi _{C}$ adapted from Princen (1979). This is the result for hexagonal foams except with $\phi _{c}$ shifted to $16\,\%$, suitable for disordered foams (Cantat et al. 2013, p. 195), and $R$ replaced by $R_{2 1}$:

(2.14)\begin{equation} \varPi_{C} \approx \frac{\gamma_\infty}{R_{2 1}} \sqrt{\frac{1 - \phi}{1 - \phi_{c}} \frac{\phi_{c}}{\phi}} . \end{equation}

A measurement of $\varPi _{C}$ in the relaxed foam, obtained from (2.2) directly, is used to reset the strength of $\varPi _{D}$, and the foam is relaxed again. We find that around five iterations of this process are needed for $\varPi _{C}$ to converge to precision $10^{-2} \gamma _\infty / \sqrt {A_{F}}$ (where $A_{F}$ is the area of the foam's periodic domain). Once this is done, the desired equilibrium film thickness $h_0$ is approximately obtained. However, as mentioned in § 2.1.2, the film thickness will vary between bubbles, due to their variations in pressure (Princen 1988).

Finally, the deviatoric stress is relaxed again (now with a larger tolerance $10^{-3}\gamma _\infty / \langle R \rangle$, due to the extra computational cost for a wet foam), giving an initial wet foam structure that is intended to be representative of a bulk foam. In the next section, we analyse coarsening-related bubble properties within such foam samples. We use five distinct samples containing $256$ bubbles, along with a single sample of $1024$ bubbles, which are analysed at different liquid fractions and contact angles, as described below. As noted in § 1, we do not coarsen these samples with time (i.e. evolve the bubble areas according to their growth rates) due to the small system sizes (Thomas et al. 2006).

3. Results

We now present our main results, on coarsening-related properties of two-dimensional disordered wet foams, for various liquid fractions $\phi$ and contact angles $\theta$.

Our simulations follow a similar method to that of Bolton & Weaire (1990) in gradually increasing $\phi$, and relaxing the foam's structure at each step (we also relax the film thickness and applied stress at every $\phi$). Their approach has been applied by Cox et al. (2018), Jing et al. (2021), Feng et al. (2021) and Jing & Feng (2023). We thereby capture the flocculation (Cox et al. 2018) discussed in § 3.1, whereas qualitatively different foam structures result if $\phi$ is set to a large value in a single step.

We perform these liquid fraction sweeps from $\phi = 2\,\%$ to $25\,\%$. The upper bound is arbitrary, and chosen so that data is obtained for foams with liquid fractions above the unjamming transition without bubble attraction, i.e. $\phi _{c} \approx 16\,\%$ (Cantat et al. 2013, p. 195). The lower bound is determined by convergence. Previous simulations (with $\theta > 0$) using a different structural model have used the same upper bound (Cox et al. 2018; Jing et al. 2021).

Simulations are run for $\theta \in \{0, 2.5^\circ, 5^\circ, 7.5^\circ, 10^\circ \}$, and for repulsive $\varPi _{D}$. Our range for $\theta$ is comparable to prior numerical studies (Feng et al. 2021; Jing et al. 2021), and does not exceed contact angles observed in foam experiments (Princen 1968; Seknagi 2022). Larger values of $\theta$ could be studied, but result in slower convergence and require higher mesh refinement. We also note $\theta = 10^\circ$ is small enough that $\sin \theta \approx \theta$ and $\cos \theta \approx 1$. We use these approximations while obtaining theoretical results below.

We set the equilibrium film thickness $h_0$ to $\langle R \rangle / 100$. Our simulations cannot converge for values much below $\langle R \rangle / 200$ at the selected degree of refinement, stated below, recalling our comments in § 2.1.

The simulated foams contain $256$ or $1024$ bubbles. Our energy-change threshold for halting structural relaxation is $10^{-6} \gamma _\infty \sqrt {A_{F}}$ (over $20$ conjugate gradient iterations), where $A_{F}$ is the area of the foam's periodic domain, and we use $m = 4$ refinements of the initial dry-foam mesh. Hence, the mean number of vertices per bubble is initially $96$, before the discretisation of large or small bubbles is adjusted. The liquid fraction is incremented in steps of $\Delta \phi = 0.5\,\%$. These parameters have been varied to check convergence. The growth rates at $\phi = 2\,\%$ for $m = 4$ differ from their values at $m = 5$ by approximately $10\,\%$. The difference corresponds to a uniform scaling for all bubbles to a good approximation, and may be related to convergence of the Young–Laplace law as applied to the films. This is discussed in Appendix C where we describe a convergence test in detail.

The execution time for one liquid fraction sweep with $256$ bubbles is approximately 100 h on a personal computer with a recent $16$-core Intel i7 processor, 16 GB of random-access memory, and a solid-state drive, recalling that no parallel processing is used. A sweep with $1024$ bubbles takes around 700 h. The execution times are usually lower for smaller $\theta$ and larger $\phi$, and when bubble areas are taken directly from the initial Voronoi tessellation.

In § 3.1 below, we describe the global properties of the simulated foams. Next, in § 3.2, we discuss coarsening-related properties of individual bubbles, and develop some analytical approximations. Section 3.3 concerns the variation of the simulated bubble growth rates with $\phi$, and their comparison with prior growth laws. Finally, in § 3.4, we consider the aggregated properties of the foam's bubbles, in order to better quantify their variation with $\phi$ and $\theta$.

3.1. Foam properties

Let $\bar {\boldsymbol {\tau }}_{F}$ be the stress tensor spatially averaged, after Batchelor (1970), over a two-dimensional foam. Let $A_{F}$ be the foam's total area, $A_k$ the area and $p_k$ the pressure of the $k\text {th}$ bubble, and let $\boldsymbol {t}$ be the unit tangent and $\gamma$ the surface tension at a point on the union $\varGamma _{F}$ of all liquid–gas interfaces in the foam. Also let $\delta _{i j}$ be the Kronecker delta symbol. Then the components of $\bar {\boldsymbol {\tau }}_{F}$ are (Batchelor 1970; Cantat et al. 2013, pp. 175–177)

(3.1)\begin{equation} (\bar{\boldsymbol{\tau}}_{F})_{i j} = \frac{1}{A_{F}} \left(\int_{\varGamma_{F}} \gamma t_i t_j \,{\rm d} l - \sum_k p_k A_k \delta_{i j}\right) , \end{equation}

recalling that we set the liquid pressure to zero. Our means of measuring $\bar {\boldsymbol {\tau }}_{F}$ in the simulations is given in Appendix B.1. The osmotic pressure (Princen 1979) is then given by (Hutzler & Weaire 1995; Höhler et al. 2021)

(3.2)\begin{equation} \varPi_{O} ={-}\tfrac{1}{2} \operatorname{Tr} \bar{\boldsymbol{\tau}}_{F} . \end{equation}

This is the average pressure in the foam, above that of the liquid (Höhler et al. 2021). It is also approximately the pressure that a piston would need to exert on the foam to maintain its present $\phi$, where this piston is permeable to the liquid but not the bubbles (Princen 1979; Höhler et al. 2021).

In figure 8, we plot the simulated variation of $\varPi _{O}$ with effective liquid fraction (defined below) for several $\theta$. The osmotic pressure is positive for $\phi$ up to a threshold value, at all considered $\theta$. This threshold is approximately the unjamming transition $\phi _{c} \approx 16\,\%$ (Cantat et al. 2013, p. 195) for small $\theta$, as expected (Hutzler & Weaire 1995), and decreases for larger $\theta$, as for hexagonal foams (Princen 1979). For $\phi$ above the threshold, $\varPi _{O} \approx 0$, so the foam is no longer under compression. For repulsive $\varPi _{D}$ and $\theta = 0$, the variation of $\varPi _{O}$ with $\phi$ is in qualitative agreement with prior simulation results (Hutzler & Weaire 1995), including, via figure 8(b), the apparent scaling $\varPi _{O} \sim (\phi - \phi _{c})^2$ near $\phi _{c}$. At larger $\theta$, the greater differences between individual runs at small $\varPi _{O}$, reflected in the error bars of figure 8(b), are interpreted to arise from an observed dependence on initial conditions, mentioned below, in the manner in which flocculation progresses.

Figure 8. Osmotic pressure $\varPi _{O}$ versus effective liquid fraction, defined in (3.3), in foams with different contact angles $\theta$. Data is given for $1024$-bubble runs in (a), while (b) gives the variation of $\sqrt {\varPi _{O}}$ for $256$-bubble runs to clarify the zero of $\varPi _{O}$. The data in (b) is averaged over five different initial foams, with the square root then being taken, and the error bars give the propagated sample standard deviation of the linear data. For comparison, the solid (dashed) curve is for a $1500$-bubble foam with $\theta = 10.8^\circ$ ($\theta = 5.1^\circ$) and $h_0 = 0$, simulated with the model of Cox et al. (2018). The curves in (b) end when $\varPi _{O} < 0$.

Alongside our data, we plot that obtained using the same simulation model as Cox et al. (2018), i.e. the Surface Evolver with $h_0 = 0$. Good agreement is observed. To make the comparison, we correct for our finite film thickness by plotting against the effective liquid fraction $\phi _{eff}$, instead of $\phi$, obtained by subtracting from the foam's liquid area a strip of width $h_0/2$ around the perimeter of each bubble (Princen 1979; Khan & Armstrong 1989). If $\langle P \rangle$ is the mean bubble perimeter (the perimeter $P$ is to be distinguished from the foam polydispersity $\mathcal {P}$), and $\langle A \rangle$ the mean bubble area, then

(3.3)\begin{equation} \phi_{eff} = \phi - (1 - \phi) \frac{h_0 \langle P \rangle} {2 \langle A \rangle} . \end{equation}

For a $1024$-bubble foam with $\phi = 2\,\%$, we find $\phi _{eff} \approx 1.14\,\%$, with the deviation decreasing as $\phi$ increases. The contact angle $\theta$ does not have a strong effect. A three-dimensional version of (3.3) was used by Mason et al. (1997).

In figure 9, we show the foam structure at $\phi = 25\,\%$ for several $\theta$. These foams all have $\varPi _{O} \approx 0$, recalling figure 8, and so it is clear that a variety of structures can have zero osmotic pressure. For repulsive $\varPi _{D}$, the osmotic pressure is zero because the bubbles are not touching. However, when there is bubble attraction the foam flocculates, as found by Cox et al. (2018) in simulated foams, and as is well known in emulsions (Bibette et al. 1993). A clustered structure with many contacting bubbles is energetically favourable, due to the decrease of interface tension within films, recalling (2.4) (Princen 1983; Cox et al. 2018). The presence of bubble attraction, and hence flocculation, for $\theta = 0$ is a consequence of finite film thickness, as discussed in § 2.1.2. At $\theta = 10^\circ$, different initial foam structures result in different manners of flocculation (e.g. some will form a single fissure, while others form multiple large Plateau borders).

Figure 9. Simulated structure of $256$-bubble periodic foams at $\phi = 25\,\%$, for (a) repulsive $\varPi _{D}$, (b) $\theta = 0$ and (c) $\theta = 10^\circ$. A square subset of each foam, with area equal to that of the periodic domain, is shown. Since the latter domain is deformed due to stress relaxation (§ 2.3), some simulated bubbles may be omitted or may appear twice. Each foam was generated from the same initial dry structure.

We note, from figure 8, that $\varPi _{O} < 0$ for a range of $\phi$ near $16\,\%$ at contact angle $\theta = 10^\circ$ (including for the data with $h_0 = 0$). This would correspond to the foam exerting suction on the confining piston mentioned above. As $\phi$ is further increased, $\varPi _{O}$ increases towards zero. We interpret this behaviour as arising from the energy cost for a bubble to lose a neighbour at high $\theta$. Bubbles therefore experience some extension before they detach, and this regime, in which the foam is not fully flocculated, has $\varPi _{O} < 0$. Our interpretation is supported by the observation that, when $\phi$ is decreased again from $25\,\%$, $\varPi _{O}$ remains non-negative (not shown). This hysteresis is in contrast to the behaviour for repulsive $\varPi _{D}$ (Hutzler & Weaire 1995).

Next, in figure 10, we plot the mean number of bubble neighbours $\langle n \rangle$. We measure $n$ for a bubble $B_i$ by counting the distinct bubbles from whose interfaces the vertices of $B_i$ experience a positive (i.e. repulsive) $\varPi _{D}$. The behaviour of $\langle n \rangle$ has been studied in detail for $\theta > 0$ (Feng et al. 2021; Jing et al. 2021; Jing & Feng 2023), and by Winkelmann et al. (2017) and Boromand et al. (2019) without bubble attraction (the latter using a comparable model to ours). Unlike the prior studies that we are aware of, we have used the same model for both cases.

Figure 10. Mean number of neighbours $\langle n \rangle$ versus effective liquid fraction, defined in (3.3), for foams with various degrees of bubble attraction. Data from single runs are plotted for $1024$-bubble foams, while the data for $256$-bubble foams is averaged over five runs (the error bars give the sample standard deviation). The solid curve is the PLAT data at $\theta = 0$ and $h_0 = 0$ given by Jing et al. (2021), and we also compare with $h_0 = 0$ simulations at $\theta = 10.8^\circ$ (dashed curve) and $\theta = 5.1^\circ$ (dotted curve) using the approach of Jing et al. (2021) and Cox et al. (2018) – the same runs as in figure 8.

The simulation data given by Jing et al. (2021), obtained using the PLAT software, is in good agreement for smaller $\phi _{eff}$. The deviation at larger values may be due to our finite film thickness. In figure 10, we also compare with simulations using the same model as Jing et al. (2021) for $\theta > 0$. The disagreement is much larger (surprisingly, given the consistency in figure 8), although the qualitative behaviour (Feng et al. 2021) is the same. The approximate plateaus to which $\langle n \rangle$ decreases are considerably lower in our simulations, and the failure of the decoration theorem, resulting in $\langle n \rangle < 6$ (Bolton & Weaire 1990), occurs at smaller $\phi _{eff}$. Our deviations from Jing et al. (2021) are slightly smaller, but still large, and we find less of a difference than them between zero and slight bubble attraction.

Although the simulations with $\theta > 0$ that we compare with in figure 10 have a much lower polydispersity ($\mathcal {P} = 0.008$ for $\theta = 10.8^\circ$ and $\mathcal {P} = 0.041$ for $\theta = 5.1^\circ$, for $\mathcal {P}$ as defined in § 2.3) than we use, we do not believe that this can account for the whole discrepancy. For identical initial dry foam structures, we found in tests that the variation of $\langle n \rangle$ still showed a large difference between the predictions of the simulations described here, and of those following the techniques of Jing et al. (2021) (which have zero film thickness $h_0$). Instead, we suggest that the cause may be the finite size of the transition region (Kralchevsky & Ivanov 1985b) between films and Plateau borders when $h_0 > 0$, which should limit the minimum film length (which is zero for $h_0 = 0$) to finite values, and thus decrease the maximum attractive force that a bubble may exert on a neighbour. Attractive contributions to the force come from the ends of the film, where the interface separation passes through the minimum of $\varPi _{D}$ (see figure 3), while the contributions from the rest of the film are repulsive since $\varPi _{D} > 0$ there (Denkov et al. 1995). The contributions from the film ends are constant for fixed $\theta$ (see Appendix A), whereas the repulsive contribution increases (in magnitude) with film length, as it is given by integrating $\varPi _{D}$ along the film. Hence, the maximum attractive force decreases with increasing minimum film length. A reduction in the maximum force would result in earlier detachment of adhered bubbles for $h_0 > 0$ as $\phi$ is increased, thus lowering $\langle n \rangle$ compared with the case for $h_0 = 0$.

Finally, we note that care is needed when increasing $\phi$ in our simulations with bubble attraction, due to the fragility of the bubble clusters, which can easily be broken apart if the step $\Delta \phi$ is too large. This would decrease $\langle n \rangle$, although our convergence tests imply that this is not the cause of our differences from prior data.

3.2. Bubble properties

We now turn to the properties of individual bubbles. As noted in § 1, the growth rate of a bubble is determined partly by its pressure $p$ and the length $L$ of its perimeter along which it is in contact with its neighbours, i.e. its film length (Lemlich 1978; Roth et al. 2013). Therefore, we focus on these properties, before considering the growth rates themselves in § 3.3.

3.2.1. Bubble pressures

In figure 11, we plot the scaled pressures of the bubbles against their effective radii $R$. Let $r_{PB}$ be the radius of curvature of a bubble's Plateau borders, given via the Young–Laplace law $p = \gamma _\infty / r_{PB}$ (Weaire & Hutzler 1999). Then the plotted property is $p R / \gamma _\infty = R / r_{PB}$, which is a measure of bubble deformation. From the figure, this ratio approaches unity as $R \to 0$, i.e. the smallest bubbles are circular and undeformed. These are analogous to ‘rattlers’ which exist in the Plateau borders of a polydisperse foam (Khakalo et al. 2018; Galvani et al. 2023). At smaller liquid fractions, $R / r_{PB}$ increases with $R$, so larger bubbles experience relatively more deformation (Bolton & Weaire 1991). This is illustrated in figure 12. At higher liquid fractions, $R / r_{PB} \approx 1$ for all $R$, as expected for repulsive $\varPi _{D}$, because all bubbles are approximately circular at $\phi$ near and above that of the unjamming transition (Durian 1995). However, we see from figure 11(b) that $R / r_{PB} \approx 1$ also holds in foams with bubble attraction at large $\phi$ (i.e. flocculated foams), consistent with theory for ordered foams (Princen 1979), although there is a greater spread in $R / r_{PB}$ for larger $\theta$ at all $\phi$.

Figure 11. Scaled pressure $p R / \gamma _\infty = R / r_{PB}$ versus effective radius $R$, where $r_{PB}$ is the Plateau border curvature radius, for individual bubbles in $1024$-bubble foams at various liquid fractions. Data for (a) repulsive $\varPi _{D}$ and (b) $\theta = 10^\circ$ is shown, and compared with (3.8) using the osmotic pressure measured in the simulations. We also compare with a version of (3.8) for which the factor $1/(1 - \phi )$ mentioned in the text is retained.

Figure 12. Illustration of the observation that larger bubbles are more deformed in a wet foam (Bolton & Weaire 1991), using a simulation with $\phi = 3.5\,\%$ and repulsive $\varPi _{D}$. The smallest bubbles are almost circular.

The relationship between $R$ and $p R / \gamma _\infty$ appears linear to a good approximation. To interpret this, we use the approach of Höhler et al. (2021). They generalised to arbitrary polydispersity an established relation (Princen 1988) between the osmotic and capillary pressures in three-dimensional monodisperse foams. The same proof in two dimensions, extended from $\theta = 0$ to our considered contact angles $\theta \lesssim 10^\circ$, for which $\cos \theta \approx 1$, gives

(3.4)\begin{equation} \varPi_{C} \approx \frac{\varPi_{O}}{1 - \phi} + \frac{\gamma_\infty}{R_{2 1}} . \end{equation}

This incorporates an approximation of each bubble's perimeter by that of a circle with equal area (Kraynik et al. 2012). A verification that (3.4) holds in our simulations is given in figure 13. We now adapt the argument of Höhler et al. (2021), which they apply over the entire foam, to a single bubble of effective radius $R$ and pressure $p$, and make additional simplifying assumptions. Let $\bar {\boldsymbol {\tau }}$ be the stress tensor spatially averaged over the considered bubble, including its interface $\varGamma$. Let $h$ be the local interface separation at a point on $\varGamma$, with $\gamma (h)$ the corresponding surface tension and $\boldsymbol {t}$ the unit tangent. Let $A = {\rm \pi}R^2$ be the bubble's area. Hence, in component form, (Batchelor 1970; Cantat et al. 2013, pp. 175–177)

(3.5)\begin{equation} \bar{\tau}_{i j} = \frac{1}{A} \oint_\varGamma \gamma(h) t_i t_j \, {\rm d} l - p \delta_{i j} , \end{equation}

where $\delta _{i j}$ is the Kronecker delta. We note that averaged stresses over individual bubbles in a dry foam were considered by Kraynik et al. (2012). Also let $\bar {p} = -\tfrac {1}{2} \operatorname {Tr} \bar {\boldsymbol {\tau }}$ be the bubble's reduced pressure, defined analogously to $\varPi _{O}$ (Höhler et al. 2021). This reduced pressure differs from the bubble's gas pressure $p$. Using (3.5), we obtain

(3.6)\begin{equation} \bar{p} = p - \frac{1}{2 {\rm \pi}R^2} \oint_\varGamma \,{\rm d} l \,\gamma (h) . \end{equation}

Next, we evaluate the integral, using that $\gamma (h) \approx \gamma _\infty$ in the films, from (2.4) and recalling that we have $\cos \theta \approx 1$ since $\theta \lesssim 10^\circ$. We also approximate the bubble perimeter as that of a circle with the same area. The error in doing so is around $6\,\%$ for a simulated dry foam (Kraynik et al. 2012), and expected to be smaller for larger $\phi$. Hence, we obtain

(3.7)\begin{equation} \bar{p} \approx p - \frac{\gamma_\infty}{R}. \end{equation}

As a first approximation, we take the reduced pressure $\bar {p}$ to be the same for all bubbles. Since $\bar {p}$ is related to the forces applied externally to the bubble (Landau et al. 1986, p. 7), as shown in Appendix A, this corresponds to a mean-field assumption (Lemlich 1978; Feitosa et al. 2006) that the environment of each bubble is comparable. Substituting (3.7) into (3.4), and using (2.2) for $\varPi _{C}$, gives $\bar {p} = \varPi _{O} / (1 - \phi )$. Hence, for $\theta \lesssim 10^\circ$,

(3.8)\begin{equation} p \approx \varPi_{O} + \frac{\gamma_\infty}{R} , \end{equation}

omitting the factor of $1-\phi$ for simplicity, which we find is a good approximation for all considered $\phi$ in two dimensions (see figure 11). This is the single-bubble analogue of the more rigorous (3.4). Its agreement with the data, shown in figure 11, indicates that bubbles approximately obey an effective Young–Laplace law, with a mean-field external pressure $\varPi _{O}$ (whose value is currently taken from the simulations, but not fitted). Agreement also exists when the bubble areas are taken directly from a Voronoi tessellation (data not shown). We therefore suggest that the pressure environment of the bubbles, determined by their neighbours, is comparable throughout the foam (Jorjadze, Pontani & Brujic 2013), and does not have a strong systematic variation with the bubble size.

Figure 13. Capillary pressure versus liquid fraction for $1024$-bubble foams with repulsive $\varPi _{D}$ and $\theta = 10^\circ$. Comparison is made to (3.4) for both cases, using $\varPi _{O}$ as measured in the simulations at each liquid fraction step. The prediction of $\varPi _{C}$ is linearly interpolated between these steps for clarity.

We are not aware of prior discussions of (3.8) for polydisperse wet foams in the literature, although a comparable exact result for dry foams (when approximated as having straight films) is known (Li et al. 2023). The osmotic pressure diverges in this limit, so the relation is between the radius of a bubble and its pressure relative to some zero. It would be interesting to see, through similar simulations, whether this type of simple approximation also holds in three-dimensional polydisperse wet foams. As shown in the next section, (3.8) may be applied to develop models for other bubble properties.

3.2.2. Film lengths

We now consider the total length $L$ of the thin films adjoining a given bubble in the simulated foams, which we recall is another determinant of the bubble's growth rate. Noting figure 1(b), $L$ is the length of the bubble's liquid–gas interface in contact with other bubbles, and we describe our method for measuring $L$ in Appendix B.3. Comparably to Denkov et al. (1995), we measure $L$ in our simulations by approximating each film interface with a circular arc, and determining the latter's intersection with (or closest approach to) circular arcs approximating the two adjoining Plateau-border interfaces (cf. figure 2). These interfaces are expected to approach circular arcs for high levels of mesh refinement, by the Young–Laplace law (Weaire & Hutzler 1999). The bubble film lengths are plotted in figure 14 for small and large $\phi$. We plot the ratio of $L$ to the bubble perimeter $P$ in figure 14(a), i.e. the proportion of bubble interface adjoining a thin film. The length per film, relative to $P$, is given in figure 14(b).

Figure 14. Ratio of film length to perimeter versus effective radius for individual bubbles (with $n$ neighbours) in a $1024$-bubble foam, (a) at $\phi = 2\,\%$ with repulsive $\varPi _{D}$, and (b) at $\phi = 25\,\%$ with $\theta = 10^\circ$. Comparison is made with (3.9), with $\varPi _{O}$ taken from the simulations in (a), and assumed to be zero in (b) due to flocculation. The relative length per film $L / (n P)$ is plotted in (b), since this is predicted directly by (3.9) for $\varPi _{O} = 0$.

At $\phi = 2\,\%$, figure 14(a) shows that $L / P$ decreases with decreasing bubble radius $R$, in agreement with the observation that smaller bubbles experience less deformation (Bolton & Weaire 1991). As $R$ increases, $L / P$ appears to saturate. No qualitative variation in this behaviour with $\theta$ was observed. At large $\phi$, where $L$ is smaller due to decreased bubble deformation, the film proportion depends strongly on $\theta$ due to flocculation, as shown in § 3.4.1.

Certain features of $L / P$ may be captured by an analytic approximation, as for the bubble pressures. Let $\bar {R} \equiv R / R_{2 1}$ and $\bar {\varPi } \equiv \varPi _{O} R_{2 1} / \gamma _\infty$ (Princen 1988) be dimensionless measures of bubble size and osmotic pressure, respectively. In Appendix A, we argue that

(3.9)\begin{equation} \frac{L}{P} \approx 2 \frac{\bar{\varPi} \bar{R} + n \theta / {\rm \pi}}{1 + \left(2 \bar{\varPi} + 1\right) \bar{R}} . \end{equation}

This is derived by relating the reduced bubble pressure $\bar {p}$ to the forces applied to the bubble's interface by its neighbours, and approximating this pressure by $\varPi _{O}$ as previously. We assume $\theta \lesssim 10^\circ$, as considered in our simulations, so that $\sin \theta \approx \theta$ and $\cos \theta \approx 1$. The predictions of this equation are compared with the simulated values of $L$ in figure 14.

The general variation of $L$ with $R$ at low liquid fractions and for repulsive $\varPi _{D}$ is captured, although the spread at a given $R$ is not. The approximations used to derive (3.9) effectively average the explicit dependence upon $n$. In particular, we neglect the detailed geometry of a bubble's films, which we expect to vary with $n$ due to constraints from Plateau's laws (Weaire & Hutzler 1999; Roth et al. 2013), at small $\phi$ where the decoration theorem holds (Bolton & Weaire 1991). Hence, the bands in figure 14(a) are not reproduced.

An isotropic model (Roth et al. 2013) might predict the spread in $L$ at small $\phi$, including the bands for different neighbour numbers $n$. However, its reliance on the decoration theorem precludes applicability to moderate $\phi$, which is our regime of interest.

For foams with $\theta > 0$, an approximate relation between $n$ and $R$ in wet foams is needed for (3.9) to predict $L$ as a function of $R$ only. This is beyond our scope, and so we gauge (3.9) via $L / (n P)$ (the relative length per film) in flocculated foams for which $\varPi _{O} \approx 0$. Again, a broad agreement with the data is seen in figure 14(b). We emphasise that, in general, the explicit $n$ term does not account for the whole dependence of $L$ on $n$, such as the bands in figure 14(a), due to the mentioned averaging. These bands are present at small $\phi$ for all $\theta$ in the considered range.

Equation (3.9) predicts that $L/P$ increases with $n$ (all else being fixed) when $\theta > 0$, whereas figure 14(a) shows a decrease with $n$ for repulsive $\varPi _{D}$. The latter pattern is still present for $\theta > 0$, at small $\phi$ (data not shown). We interpret it as follows. Consider an idealised Plateau border adjoining $n$ bubbles, each of which has equal pressure so that the border has $n$-fold rotational symmetry. The area of this Plateau border grows with $n$ for fixed pressures, and so a bubble with given radius $R$ will experience less deformation (i.e. will have smaller $L / P$) when inserted into the border if $n$ is larger. However, this level of geometric detail is omitted from the derivation of (3.9), as noted above. Hence, the remaining dependence on $n$ in (3.9) should be interpreted as an approximate correction for $\theta > 0$, increasing the length of the films due to attraction between the bubbles (Feng et al. 2021), as discussed further in § 3.4.1.

The behaviour of $L$ appears more complex than that of $p$, described previously. Nevertheless, we show in § 3.4.1 that an approximate averaging of (3.9) over all bubbles describes the mean film lengths very well.

3.3. Bubble growth rates

We now summarise the simulated bubble growth rates, recalling that these include gas transfer across films and through Plateau borders, as described in § 2.2. We describe in § 3.3.1 how the growth rates change as $\phi$ is increased. We then compare our simulations with prior models: in § 3.3.2, we consider the relationship between a bubble's growth rate and its effective number of neighbours $n_{eff}$ (Fortuna et al. 2012), before studying the gas flow through Plateau borders in § 3.3.3.

3.3.1. Effects of liquid fraction

At $\phi = 0$, which is inaccessible to our simulations, the growth rate of a bubble with $n$ neighbours obeys von Neumann's law (von Neumann 1952),

(3.10)\begin{equation} \dot{A} = \frac{2 {\rm \pi}D H \, \gamma(h_0)}{3 h_0} \left(n - 6\right) , \end{equation}

where $D$ and $H$ are, respectively, a diffusion coefficient and Henry's constant (see § 2.2), and $\gamma (h_0)$ is the film-interface tension from (2.4). The prefactor multiplying $n - 6$, which is constant for a particular foam, is derived, for example, by Cantat et al. (2013, p. 109). Hence, a bubble's growth rate is determined only by its neighbour number $n$ in a dry foam, so that the growth rates fall into discrete bands.

When $\phi$ is increased to small values, such that almost all Plateau borders still connect exactly three films (i.e. the decoration theorem is approximately satisfied (Bolton & Weaire 1991)), the growth rates remain predominantly determined by $n$, although the reduction of gas flow by the Plateau borders induces a dependence upon bubble size and shape (Bolton & Weaire 1991; Roth et al. 2013). In figure 15(a), we plot the bubble growth rates at $\phi = 2\,\%$, noting that a sequence of bands, as predicted by (3.10) for $\phi = 0$, is still apparent at this small liquid fraction. These plotted growth rates are in qualitative agreement with experiments and theory for quasi-two-dimensional foams – in particular, for fixed $n$, the absolute growth rates tend to increase slightly with increasing bubble size (Roth et al. 2013). This implies that larger bubbles experience less blocking of gas transfer by the Plateau borders (Bolton & Weaire 1991), again consistent with their greater observed deformation. As discussed in § 3.3.2, the absolute growth rates of the smallest bubbles are enhanced by their thinner films in our simulations.

Figure 15. Individual bubbles growth rates versus effective radius in a particular $1024$-bubble foam at different liquid fractions $\phi$, with repulsive $\varPi _{D}$. The number of neighbours $n$ of each bubble is shown. In (a) $\phi = 2\,\%$; in (b) $\phi = 6.5\,\%$; and in (c) $\phi = 15.5\,\%$. The growth rates are scaled so that the rates from von Neumann's law, (3.10), are the integers. A bubble with $n = 14$ lies outside the plotted domain in (a), to better show the banding of the growth rates.

When $\phi$ is increased further, the Plateau borders grow, and the degree of border blocking increases (Roth et al. 2013). The primary effect of this is to reduce the absolute growth rates, as shown in figure 15(b). However, the decoration theorem also fails for more bubbles as they detach from their neighbours with increasing $\phi$. The geometric constraints from Plateau's laws (Weaire & Hutzler 1999) are relaxed for such bubbles, resulting in growth rates that lie between the bands hitherto occupied.

Eventually, $\phi$ is high enough that most of the bubbles adjoin at least one Plateau border which connects more than three films. As shown in figure 15(c), banding in the growth rates is no longer discernible.

The absolute growth rates continue to decrease with increasing $\phi$, except in the presence of adhesion, where they eventually decrease to a plateau due to flocculation. This is considered further in § 3.4.2. Due to the inclusion of gas flow through the Plateau borders in our model, the growth rates do not shrink to zero at the unjamming transition $\phi _{c} \approx 16\,\%$ for repulsive $\varPi _{D}$. Instead they continue to decrease in magnitude as the bubbles move apart. We expect that our approximation for the growth rates, (2.12) from § 2.2, declines in accuracy as $\phi$ increases beyond this transition, as the gas dissolved in the bulk liquid becomes the main source and sink for bubbles, rather than pairwise transfer between adjacent bubbles (Fortuna et al. 2012). However, as justified in § 3.3.3, we believe that the coarsening approximation remains valid for flocculated foams, since bubbles remain in contact.

3.3.2. Effective neighbour number

We now compare our simulated growth rates with results from the literature. Fortuna et al. (2012) introduced the effective number of neighbours $n_{eff}$ of a bubble in a two-dimensional wet foam, defined as follows. Let $\varTheta$ be the angle turned by the bubble interface's tangent as its film interfaces are traversed, illustrated in figure 16. Contributions to $\varTheta$ from a convex interface, as in this figure, are positive. Then (Fortuna et al. 2012)

(3.11)\begin{equation} n_{eff} = 6 - 3 \varTheta / {\rm \pi}. \end{equation}

This is equivalent to the bubble's actual number of neighbours $n$ for $\phi = 0$, and is equal to $6$ in a wet foam when $n = 0$ (Fortuna et al. 2012). If $n_{eff} < 6$, then (3.11) implies that the film interfaces of the bubble tend to be convex, and so its pressure tends to be larger than those of its neighbours, i.e. it is smaller than average, in a loose sense, from (3.8). The converse is true for $n_{eff} > 6$. In fact, by repeating the derivation (von Neumann 1952) of von Neumann's law in a wet foam, assuming small uniform film thickness $h_0$ and $\theta = 0$, it can be shown that the border-blocking growth rate of a bubble (i.e. including only gas flow through the films) is exactly

(3.12)\begin{equation} \dot{A}_{BB} = \frac{2 {\rm \pi}\gamma_\infty D H}{3 h_0} (n_{eff} - 6) . \end{equation}

This is von Neumann's law with $n$ replaced by $n_{eff}$, and was derived in a slightly different form by Bolton & Weaire (1991). Fortuna et al. (2012) validate this result, in conjunction with an approximation for the Plateau border gas flow, using Potts model simulations. For $\theta > 0$, the only change to (3.12) is that $\gamma _\infty$ becomes $\gamma (h_0)$ from (2.4). Since we consider $\theta \lesssim 10^\circ$, for which $\cos \theta \approx 1$, we have $\gamma (h_0) \approx \gamma _\infty$.

Figure 16. Schematic of the definition of the angle $\varTheta$ (Fortuna et al. 2012), described in the text.

We measure $\varTheta$, and hence $n_{eff}$, by noting that the contribution to $\varTheta$ from a bubble's $j\text {th}$ film is $l_j / r_j$ (Fortuna et al. 2012), where $l_j$ is the film's length (Appendix B.3), and $r_j$ is its radius of curvature. The latter is positive if the film is convex with respect to the bubble, and is defined using the augmented Young–Laplace equation (Kralchevsky & Ivanov 1985b; Denkov et al. 1995) according to Appendix B.2.

Figure 17 shows our bubble growth rates $\dot {A}$ plotted against $n_{eff}$ (including contributions from gas flow through the Plateau borders). The two quantities are closely related at all $\phi$ and $\theta$ we have investigated, provided that the bubbles have neighbours (including in flocculated foams). Otherwise $n_{eff} = 6$ identically, while the growth rates still vary.

Figure 17. Bubble growth rate versus $n_{eff}$ in a $1024$-bubble foam (a) at $\phi = 2\,\%$ for repulsive $\varPi _{D}$, and (b) at $\phi = 25\,\%$ for different $\theta$. A large bubble with $n = 14$ neighbours is omitted from (a), for clarity. The black lines are (3.12) for the set film thickness $h_0$.

This close relationship between $\dot {A}$ and $n_{eff}$ (with the sign of $\dot {A}$ being well determined by the sign of $n_{eff} - 6$) was not apparent to us from the simulations performed by Fortuna et al. (2012), likely due to the fact that the Potts model does not accurately capture bubble geometry. In principle, deviations from (3.12) give the gas transfer through Plateau borders, although the following caveats should be considered.

A separate group of bubbles may be observed in figure 17 at $\phi = 25\,\%$ and $\theta = 0$, with $n_{eff} < 6$ and growth rates considerably smaller than other bubbles at the same $n_{eff}$ (close to the vertical line). These are interpreted to result from difficulties in defining and measuring film lengths and radii of curvature in extremely short films with a film thickness comparable to their length, due to the transition regions between films and Plateau borders (Kralchevsky & Ivanov 1985b). The group is smaller at higher mesh refinement, but is not practical to eliminate.

We previously mentioned film thickness variations due to our form for $\varPi _{D}$, in § 2.1.2. The effect of these on the bubble growth rates is shown in figure 18, where we plot our simulated border-blocking growth rates against (3.12), for which the thickness of each film is the set value $h_0$. We measure the border-blocking growth rate of a bubble with pressure $p$ as follows. Let the bubble's $j\text {th}$ film (of $n$) have length $l_j$ (Appendix B.3), thickness $h_j$ (Appendix B.2), and adjoin a bubble with pressure $p_j$. Then (Bolton & Weaire 1991; Cantat et al. 2013)

(3.13)\begin{equation} \dot{A}_{BB} = D H \sum_{j=1}^n \frac{l_j}{h_j} (p_j - p) . \end{equation}

We observe from figure 18 that large deviations from (3.12) are mainly restricted to small bubbles. We interpret this by noting that smaller bubbles have higher gas pressure by (3.8), and so the disjoining pressure in their films must be larger to maintain equilibrium, recalling § 2.1.2 (Toshev & Ivanov 1975). Hence, by figure 3, these films are thinner than those of larger bubbles, enhancing the absolute growth rates of small bubbles (Princen 1988). Bubble pressure increases nonlinearly with decreasing size, so large bubbles do not show a decrease in their absolute growth rates to the same degree.

Figure 18. Simulated border-blocking bubble growth rates in a $1024$-bubble foam with repulsive $\varPi _{D}$ at $\phi = 2\,\%$, compared with their scalings to approximately account for film thickness variations between bubbles. Here $\langle h_{f} \rangle$ is a bubble's mean film thickness – a representative thickness is measured (according to Appendix B.2) for each of its films, and the mean of these thicknesses is calculated. The black line is (3.12) with $h_0$ replaced by $h$.

We approximately correct for the film thickness variations by scaling the growth rates with the measured mean film thickness $\langle h_{f} \rangle$ (see caption of figure 18), reducing the differences for small bubbles. This simple scaling cannot be done for the total growth rates, since the dependence on film thickness of Plateau border transfer is different, as seen in the next section. While discrepancies are still present, part of which we expect to come from thickness variations between the individual films of each bubble (agreement is not substantially improved at higher mesh refinement), our results are consistent with (3.12).

In figure 17(b) at $\phi = 25\,\%$, the growth rates are also affected by the fact that the foam-wide mean film thickness is larger than the set value $h_0$ (by a factor of approximately $1.59$ for $\theta = 0$). We interpret this by observing that the interfaces of short films between bubbles with equal pressure retain curvature over their whole lengths in our simulations, interpreted as being due to the finite transition regions between films and Plateau borders (Kralchevsky & Ivanov 1985b). The augmented Young–Laplace equation (Kralchevsky & Ivanov 1985b; Denkov et al. 1995), given in Appendix B.2, then requires a weaker disjoining pressure than for zero curvature, and hence a thicker film.

3.3.3. Gas transfer through Plateau borders

Schimming & Durian (2017) derive a bubble growth law in two-dimensional wet foams at small $\phi$ (satisfying the decoration theorem), which accounts for gas flow through the Plateau borders. It is assumed that $\theta = 0$, and the geometry is simplified by taking each side of each border to have equal curvature. Let $r_{PB}$ be the radius of curvature of a bubble's Plateau border interfaces (given by the Young–Laplace law in terms of the bubble pressure $p$), let $r_j$ be the signed radius of curvature of its $j\text {th}$ film (measured as in Appendix B.2), and define the circularity $C = (R / n) \sum _{j = 1}^n 1/r_j$ (Roth et al. 2013) where $n$ is the bubble's number of neighbours. Then the contribution to the bubble's growth rate from Plateau-border transfer, termed border crossing, is (Schimming & Durian 2017)

(3.14)\begin{equation} \dot{A}_{BC} \approx{-}2 {\rm \pi}\gamma_\infty D H \frac{n C}{R} \sqrt{\frac{r_{PB}}{h_0}} . \end{equation}

We believe that a typographical error resulted in an extra factor of $1/\sqrt {3}$ for Schimming & Durian (2017), which we have corrected. We plot a comparison between this approximation and our results for a repulsive disjoining pressure at $\phi = 2\,\%$ in figure 19. Part of the difference is believed to come from the transition regions at the ends of the films (Kralchevsky & Ivanov 1985b) – the approach we use to measure the film lengths (Appendix B.3) tends to include part of these regions, which are omitted in the derivation of (3.14). Therefore, the film lengths are overestimated in the simulations, and thus also the proportion of the growth rates due to gas transfer through the films. We note that the definition of film length is ambiguous for finite film thickness, due to these transition regions. Other contributions to the difference in figure 19 come from the approximations made in our coarsening model (see § 2.2) and in (3.14) from Schimming & Durian (2017). Agreement is not improved at higher mesh refinement (and is worse for the smallest bubbles, in the cluster with smallest $\dot {A}_{BC}$).

Figure 19. Border-crossing growth rate, defined as the difference between the total and border-blocking rates, versus $n_{eff}$ in a $1024$ bubble foam at $\phi = 2\,\%$ with repulsive $\varPi _{D}$. The simulations are compared with (3.14).

We also recall that our selected $h_0$ is large compared with real foams. This is instead expected to exaggerate the border-crossing rates relative to the gas flow through the films in our simulations (Schimming & Durian 2017), which is a separate effect from the underestimation at a given $h_0$, discussed above.

Near to the unjamming transition, at which the film lengths decrease to zero (Cantat et al. 2013, p. 195), the bubble growth rate due to one neighbour was given by Schimming & Durian (2017) for $\theta = 0$. Let $R_i$ be the effective radius of the bubble's $i\text {th}$ neighbour, and $h_i$ the minimum distance between their interfaces. By adding contributions from each neighbour, and neglecting the effect of other neighbours on each contribution, which would tend to reduce the absolute growth rates (Evilevitch et al. 2002),

(3.15)\begin{equation} \dot{A} \approx {\rm \pi}\gamma_\infty D H \sum_{i = 1}^n \sqrt{\frac{2}{h_i}} \frac{R - R_i}{\sqrt{R R_i (R + R_i)}}. \end{equation}

This is compared with our simulations in figure 20, which show good agreement, provided that contributions from all nearby bubbles are included – not just those in contact, but also those nearly in contact. Nearby bubbles are taken here to be those counted as neighbours by the vertex neighbour algorithm (see § 2.1.3 and figure 5), and contacting bubbles are those which exert $\varPi _{D} > 0$ on a vertex of the considered bubble.

Figure 20. Total bubble growth rates plotted against (3.15) (Schimming & Durian 2017) for a $1024$-bubble foam at $\phi = 16\,\%$ (set in one step, rather than increased gradually) with repulsive $\varPi _{D}$. We do not expect the means of setting $\phi$ to have a large effect for this $\varPi _{D}$. The contacting and nearby bubbles are as specified in the main text. The black line has unit gradient.

While the issue of accounting for long-range gas transfer in simulations (R. Höhler, personal communication) remains unresolved (for which our approximation described in § 2.2 is likely poor), and may be particularly important due to our thick films, we observe here that the two different approaches to the bubble gas flow rates are consistent. Equation (3.15) uses an approximation for small $h_i$, although we find the exact prediction (Schimming & Durian 2017), which gives larger absolute growth rates, differs by less than $10\,\%$ for $h_1 \leq 2 R$ when $R_1 \leq 10 \, R$ and $n = 1$. We expect that most of the difference between the two datasets plotted in figure 20 is due to bubbles which are only just out of contact, noting from figure 5 that the vertex neighbour algorithm can count fairly distant bubbles as neighbours.

Our use of the growth rate approximation (§ 2.2) for flocculated foams, in which bubbles retain close contacts as $\phi$ is increased, is also supported by the above results for foams with repulsive $\varPi _{D}$ near unjamming. This is because the geometry is similar except for the retention of short films in the flocculated foams, while we expect the coarsening approximation to be effective for films as noted in § 2.2.

3.4. Statistics of bubble properties

We now give averaged properties of the bubbles in the simulated foams, in order to quantify some of our observations regarding their variation with $\phi$ and $\theta$.

3.4.1. Film lengths

First, we consider the averaged bubble film lengths $L$. Recalling our approximation, (3.9), we may estimate its prediction of the mean ratio of film length to perimeter, for $\theta \lesssim 10^\circ$, as

(3.16)\begin{equation} \left\langle \frac{L}{P} \right\rangle \approx \frac{\bar{\varPi} + \langle n \rangle \theta / {\rm \pi}}{1 + \bar{\varPi}} \end{equation}

by setting $R = R_{2 1}$ and $n = \langle n \rangle$ on the right-hand side. For $\theta = 0$, this has an equivalent form to an existing approximation for the film areas in a monodisperse three-dimensional foam (Höhler et al. 2021; Pasquet et al. 2023b).

Equation (3.16) is compared with our simulations in figure 21(a), and we observe a very close agreement, including in flocculated foams. We note that our simulations with repulsive $\varPi _{D}$, rather than $\theta = 0$, should be compared with the predictions of (3.16) for $\theta = 0$, since the simulations with $\theta = 0$ include some bubble attraction due to our finite film thickness $h_0$. We have assumed $h_0 \to 0$ in the derivation of (3.9) (see Appendix A), and thus (3.16). From figure 21, the film lengths initially decrease with liquid fraction, as the bubbles become less deformed. For repulsive $\varPi _{D}$, the lengths vanish as the bubbles separate at the unjamming transition, while they decrease to an approximate plateau for $\theta \geq 0$ due to flocculation. Our results for $\theta > 0$ are in qualitative agreement with those of Feng et al. (2021), who measured the total length of all films in foams simulated using a different model, and are consistent with the interpretation that $\theta > 0$ caused films to be retained at $\phi > \phi _{c}$ in the ISS coarsening experiments (Pasquet et al. 2023b).

Figure 21. Mean ratio of film length to perimeter versus $\phi$, for repulsive $\varPi _{D}$ alongside various $\theta$. In (a), simulations with $1024$ bubbles are compared with (3.16), and error bars give the standard deviation of the simulated $L / P$ distribution. The values of $\bar {\varPi }$ and $\langle n \rangle$ in (3.16) are taken from the simulations at each liquid fraction step, with linear interpolation of the predicted $\langle L / P \rangle$ between these steps. In (b), single runs with $1024$ bubbles are shown alongside the mean of five runs with $256$ bubbles. Error bars in (b), giving the sample standard deviation of $\langle L / P \rangle$ in the latter simulations, are omitted since they are close to the markers in size.

We note that (3.16) predicts that the plateau satisfies $\langle L / P \rangle \approx \langle n \rangle \, \theta / {\rm \pi}$, since $\varPi _{O} \approx 0$ in a flocculated foam (Princen 1983). Similar predictions of the mean film area in flocculated three-dimensional foams have been made, which instead vary as $\theta ^2$ (Durian 2023; Pasquet et al. 2023b).

The error in (3.16) is larger for smaller $\phi$, likely due to the greater deformation of bubbles, and the geometric simplifications made while deriving (3.9) (see Appendix A). We note that the exact mean $\langle L / P \rangle$ predicted by (3.9) behaves similarly, although is a slightly worse approximation to the data in the plateau.

3.4.2. Growth rates

Next, we consider the averaged bubble growth rates. Noting that $\langle \dot {A} \rangle \equiv \langle {\rm d} A / {\rm d} t \rangle = 0$ because the system is closed and the gas is taken as incompressible (neglecting that our coarsening approximation does not conserve the gas exactly, as discussed in § 2.2), we consider the root-mean-square growth rate $\sqrt {\langle \dot {A}^2\rangle } \equiv \sqrt {\langle ({\rm d}A / {\rm d}t)^2 \rangle }$. Since this measures typical instantaneous bubble growth rates (in absolute value), we use it as a proxy for the rate at which the foam would coarsen at early times.

In figure 22, we plot the variation of $\sqrt {\langle \dot {A}^2\rangle }$ with $\phi$, for several values of $\theta$. As expected, the coarsening rate decreases with increasing $\phi$ for small liquid fractions, since the Plateau borders grow larger and further frustrate the flow of gas (Roth et al. 2013). The contact angle does not have a strong effect in this regime.

Figure 22. Root-mean-square growth rate versus $\phi$, for repulsive $\varPi _{D}$ and various $\theta$. Results for single $1024$-bubble runs are shown, along with the mean of five runs with $256$ bubbles (the error bars give the sample standard deviation).

For repulsive $\varPi _{D}$, the coarsening rate decreases for all $\phi$, since the gas flow between separated bubbles (see figure 9) decreases as their separation increases (Schimming & Durian 2017). However, as observed for the film lengths, the presence of bubble attraction causes the coarsening rate to decrease to a plateau. As noted in § 3.1, the bubbles flocculate (Cox et al. 2018), thereby retaining films that allow for efficient gas transfer. The typical film length increases with $\theta$, and so the plateau values likewise increase. We recall that there is bubble attraction for $\theta = 0$ in our model, due to the finite film thickness. Hence, we are consistent with reports that flocculation in emulsions increases the coarsening rate (Djerdjev & Beattie 2008), along with the interpretation that $\theta > 0$ caused the $\langle R \rangle$ growth exponent (see § 1) to remain equal to $1/2$ (i.e. its value for dry foams) at higher $\phi$ than expected in experiments on the ISS (Pasquet et al. 2023b). However, obtaining growth exponents from our simulations would require following coarsening through time, which is not feasible for large enough systems at present.

For similar reasons, we have not investigated the relationship between $\sqrt {\langle \dot {A}^2\rangle }$ and the coarsening rate at early times $t$. For a system with many bubbles, and a sufficiently coarse time resolution, the discontinuous increases of $\langle A \rangle$ when bubbles vanish, with ${\rm d} \langle A \rangle / {\rm d} t = 0$ between these events, are averaged to give a differentiable function $\langle A \rangle$ of time (Cantat et al. 2013, p. 105). In this limit, the coarsening rate could be measured by ${\rm d} \langle A \rangle / {\rm d} t$ for example, whereas $\langle \dot {A} \rangle = 0$ as noted above.

The early-time coarsening rate might differ qualitatively from $\sqrt {\langle \dot {A}^2\rangle }$ in its dependence on $\phi$, since, for example, it has been suggested that the structure of flocculated foams may change substantially under coarsening (Pasquet et al. 2023b). Furthermore, based on a comparison with simulations using bubble areas directly from a Voronoi tessellation, we believe that the variation shown in figure 22 depends quantitatively on polydispersity, but not qualitatively. Since we expect the foam's polydispersity to initially change with time during coarsening (Thomas et al. 2006), we do not believe that $\sqrt {\langle \dot {A}^2\rangle }$ will be constant in time.

Finally, we recall (3.13) for a bubble's border-blocking growth rate, which states that the contribution from each film is proportional to the pressure difference between a bubble and its neighbour, multiplied by the ratio of film length to thickness (Bolton & Weaire 1991; Cantat et al. 2013). Equation (3.8) can be used to approximate the bubble's pressure $p$, and we estimate the neighbour pressures as equal and constant throughout the foam (Pasquet et al. 2023b). Again using (3.8), the neighbour pressure is expressed in terms of a critical radius $R_0$ (Lemlich 1978; Chieco & Durian 2023), whose value is obtained as described below. Equation (3.9) approximates the bubble's total film length $L$, and we use $P \approx 2 {\rm \pi}R$ (Kraynik et al. 2012) therein. Hence, we obtain an approximate border-blocking bubble growth law

(3.17)\begin{equation} \dot{A}_{BB} \approx \frac{4 {\rm \pi}\gamma_\infty D H}{h_0} \left(\frac{R}{R_0} - 1\right) \frac{\bar{\varPi} \bar{R} + n \, \theta / {\rm \pi}}{1 + (2 \bar{\varPi} + 1) \bar{R}} , \end{equation}

where the foam's critical bubble radius $R_0$ is determined by the condition $\langle \dot {A}_{BB} \rangle = 0$ mentioned above (Lemlich 1978; Pasquet et al. 2023b). We assume $\theta \lesssim 10^\circ$ in (3.8) and (3.9), and hence also in (3.17). This gives the growth rate solely in terms of the bubble radius and foam properties, notwithstanding the aforementioned need to relate $n$ to $R$, and thus is a (two-dimensional) generalisation of the Lemlich (1978) approach to allow for a variation in relative film length with bubble size. Recalling figure 14, smaller bubbles have shorter films relative to their size, and thus experience a larger degree of border blocking (Bolton & Weaire 1991). While (3.17) does not accurately give the growth rates of individual bubbles, as shown in figure 23, it may be fruitful to analyse the size distributions it predicts in the coarsening scaling state (Lemlich 1978; De Smet, Deriemaeker & Finsy 1997; Yarranton & Masliyah 1997), particularly whether it can reproduce the large population of small bubbles observed recently by Galvani et al. (2023) through the enhanced border blocking of such bubbles.

Figure 23. Border-blocking growth rates versus effective radius in a $1024$-bubble foam with repulsive $\varPi _{D}$. The liquid fraction is (a) $\phi = 2\,\%$ and (b) $\phi = 10\,\%$. Comparison is made with (3.17), where $R_0$ is fitted as stated in the text, and $\varPi _{O}$ is as measured in the simulations. The simulated growth rates are scaled by the average film thickness $\langle h_{f}\rangle$ for each bubble, as described in § 3.3.2, and $h_0$ is replaced by $\langle h_{f} \rangle$ in (3.17), to improve agreement for small bubbles.

In proposing (3.17), we follow the suggestion of Pasquet et al. (2023b) to investigate growth laws which generalise the Lemlich (1978) approach by accounting for correlations between bubble pressures and film lengths. A three-dimensional growth law of a different form, derived by considering gas flow through the bulk liquid (rather than the films), but which also incorporates correlations between pressure and efficiency of diffusion, was studied numerically by Yarranton & Masliyah (1997). Another three-dimensional growth law, again of a form different to (3.17), was studied by De Smet et al. (1997).

4. Conclusions

We have described simulations for studying coarsening in two-dimensional wet aqueous foams, with a structural model based upon the work of Kähärä et al. (2014) and Boromand et al. (2018), and a coarsening approximation inspired by the analytical work of Marchalot et al. (2008) and Schimming & Durian (2017). Our methods allow for arbitrary liquid fractions $\phi$, and support contact angles $\theta$ up to approximately $10^\circ$. The interfaces interact through a model disjoining pressure, although our selected form could be swapped for another.

We have applied this model to polydisperse foams, containing $256$ or $1024$ bubbles. The liquid fraction $\phi$ has been varied by gradually increasing its value within particular foam samples, following the approach of Bolton & Weaire (1990) and Cox et al. (2018), but using a different structural model. We have analysed the coarsening-related properties of these foams. For the considered range of contact angles $\theta \lesssim 10^\circ$, the bubble pressures $p$ were found to obey an effective Young–Laplace law relating them to the osmotic pressure $\varPi _{O}$, and we derived an approximation for the bubble film lengths $L$. When the averaged film length is estimated from the latter result, it takes an equivalent form to a previously proposed approximation in three dimensions for $\theta = 0$ (Höhler et al. 2021; Pasquet et al. 2023b), and also agrees well with our simulations.

Turning to the coarsening-induced bubble growth rates, we showed that their root-mean-square values decrease to a plateau with increasing $\phi$ in foams with attractive bubble interactions, due to the previously simulated (Cox et al. 2018) flocculation of bubbles. The effective neighbour number $n_{eff}$ of Fortuna et al. (2012) was found to be closely related to the growth rate (determining whether bubbles grow or shrink to a good approximation) whenever bubbles had films between them, caused by either the osmotic pressure or flocculation.

We also compared our simulations of gas flow through Plateau borders with existing predictions (Schimming & Durian 2017). Qualitative agreement was found for bubbles at small $\phi$, along with quantitative agreement at the unjamming transition for zero bubble attraction. An underestimate of the border gas flow rates at small $\phi$ is believed to arise partly from ambiguity in the film lengths due to our finite film thickness. Approximations made in our coarsening model and in the predictions of Schimming & Durian (2017) also contribute.

Our simulations remain computationally expensive, and so we have not considered the time-dependence of foams under coarsening here. However, it may be possible to obtain the time scale over which coarsening occurs in small systems, and to observe its relationship to the root-mean-square growth rates we considered. Such time scales have been investigated by Khakalo et al. (2018), without bubble attraction.

The approach we use seems most suited to accurately simulating foams on a small scale (with around $1000$ bubbles), in order to develop models of their properties. These could then be implemented in large-scale coarsening simulations using the bubble model (Durian 1995), for example. The analytical approximations we have presented here may be a starting point for this process, recalling that our expressions for the bubble pressures and film lengths can be combined to give an approximate growth law. As noted above, the properties of the resulting scaling states could be analysed (Lemlich 1978; De Smet et al. 1997; Yarranton & Masliyah 1997), which may contribute to the interpretation of experiments of coarsening in wet foams performed in microgravity (Born et al. 2021; Galvani et al. 2023; Pasquet et al. 2023b).

As noted in § 2.1.5, the implementation of our simulation model could be improved, including the equilibration algorithm, which may result in larger foams becoming feasible to simulate. The coarsening approximation could also be tested against accurate solutions of Laplace's equation within the foam's liquid, partly to gauge the importance of long-range gas transfer (R. Höhler, personal communication) for which the approximation is expected to be poor. The equilibrium film thickness $h_0$ could be varied to determine its effect on our results, along with our small derivative of the disjoining pressure at $h_0$. These changes to the disjoining pressure could be investigated using smaller systems at higher mesh refinement, to improve convergence. An alternative interface extrapolation or interpolation might also resolve the imbalance in film thickness noted in Appendix C. We expect our large $h_0$ to have increased the proportion of gas transfer through Plateau borders (Schimming & Durian 2017).

Our analytical approximations could be extended by investigating the relationship between bubble radius and neighbour number (Durand et al. 2011), and between the osmotic pressure and the mean neighbour number. By using an independent empirical expression for the osmotic pressure, similar to those described by Höhler et al. (2021), these approximations could be made fully predictive. Finally, as noted by Pasquet et al. (2023b), coarsening in flocculated foams and emulsions at large $\phi$ has not been widely studied. We note that the structure of flocculated systems differs greatly as the contact angle is varied (Bibette et al. 1993).

Another extension is to perform similar studies of three-dimensional wet foams (albeit for fewer bubbles). We have adapted the simulation model to this case, noting that the approach of Boromand et al. (2018) has been similarly adapted by Wang et al. (2021). A comparable three-dimensional approach for biological cells was also used by Van Liedekerke et al. (2020). Such models have not yet been applied to the simulation of foams, to our knowledge, although simulations with zero film thickness have been performed by A. M. Kraynik (personal communication).