Implementation of Morse-Witten theory for a polydisperse wet 2D foam simulation

TL;DR

This paper demonstrates how Morse-Witten theory can be practically applied to simulate polydisperse wet 2D foams, producing results consistent with direct calculations and outlining future steps for 3D modeling.

Contribution

It introduces a practical implementation of Morse-Witten theory for polydisperse 2D foam simulations, extending its applicability.

Findings

Equilibrated 2D foam structures match direct calculations

The method accommodates bubble size polydispersity

Outlines development of a 3D model in future work

Abstract

The Morse-Witten theory (D. Morse and T. Witten, EPL 22 (1993) 549-555) provides a formulation for the inter-bubble forces and corresponding deformations in a liquid foam, accurate in the limit of high liquid fraction. Here we show how the theory may be applied in practice, including allowing for polydispersity in the bubble sizes. The resulting equilibrated 2D structures are consistent with direct calculations, within the limitations of the theory. The path to developing a 3D model is outlined for future work.

Equations

$$\rho (\theta )=R_0+\delta R(\theta ),$$

$$\delta R(\theta )=\frac{R_{0}F}{2\gamma \pi }g(\theta )\text{,with }g(\theta )=(\pi -\theta )\sin (\theta )-\frac{\cos (\theta )}{2}-1.$$

$$-\delta R_1(0)=3R_{1}f/4\pi$$

$$F=2l{p}_0$$

$$F\simeq {\theta }_i{R}_i\gamma \left(\frac{R_1+R_2}{R_1{R}_2}\right),$$

$$x_i=R_i-\cos ({\theta }_i)\rho ({\theta }_i)\simeq -\delta R_i({\theta }_i)$$

$$x_i(f)=\frac{3f{R}_i}{4\pi }-\frac{f^2{R}_i}{2(2+R_1/R_2+R_2/R_1)},$$

$${\Delta }_{12}=(R_1-x_1(f))+(R_2-x_2(f)).$$

$$1-\frac{{\Delta }_{12}}{2R_0}=\frac{x_1(f)+x_2(f)}{2R_0}=\frac{3f}{4\pi }-\frac{f^2}{8}\left[1-{\left(\frac{\Delta R}{R_0}\right)}^2\right].$$

$$x_{ij}=-\frac{R_i}{2\pi \gamma }\left(\sum _k^N{F}_{ik}g(\Delta {\theta }_{jk})\right)-\frac{R_i{F}_{ij}^2}{2{\gamma }^2(2+R_i/R_j+R_j/R_i)}.$$

$$R_i-x_{ij}+R_j-x_{ji}=|{\mathbf{c}}_i-{\mathbf{c}}_j|.$$

$$F_{ij}=F_{ji}$$

$$\sum _j^N{F}_{ij}\frac{{\mathbf{c}}_j-{\mathbf{c}}_i}{|{\mathbf{c}}_j-{\mathbf{c}}_i|}=0.$$

$$F^{(n+1)}=a{F}^{(n+1)}+(1-a)F^{(n)}$$

$$x^{(n+1)}_{ij}=x^{(n)}_{ij}+\frac{R_{j}}{R_{i}+R_{j}}\Big{[}(R_{i}-x^{(n)}_{ij}+R_{j}-x^{(n)}_{ji})-|\mathbf{c}_{j}-\mathbf{c}_{i}|\Big{]}$$

$${\mathbf{c}}_i^{(n+1)}={\mathbf{c}}_i^{(n)}+b\sum _j^N{F}_{ij}^{(n+1)}\frac{{\mathbf{c}}_i-{\mathbf{c}}_j}{|{\mathbf{c}}_i-{\mathbf{c}}_j|},$$

$$x_{ij}=\frac{R_{j}}{R_{i}+R_{j}}\Big{[}\rho(\theta_{ij})+\rho(\theta_{ji})-|\mathbf{c}_{j}-\mathbf{c}_{i}|\Big{]}$$

$$\epsilon =\frac{1}{4\pi R_0\gamma }\sum _{i=0}^N\sum _j{x}_{ij}{F}_{ij},$$

$$\frac{F}{\gamma }=\frac{6Z}{\pi }\frac{\delta R}{R_0}$$

$$\epsilon (\Delta \varphi )={\left(\frac{3}{\sqrt{2\pi }}\frac{\Delta \varphi }{\Delta \varphi +{\varphi }_h}\right)}^2,$$

$$\delta R(\theta )=\frac{-F}{\gamma R_0}G(\theta )$$

$$G(\theta )=-\frac{1}{4\pi }\left\{\frac{1}{2}+\frac{4}{3}\cos \theta +\cos \theta \ln [{\sin }^2(\theta /2)]\right\},$$

$$x_1(F)=\frac{F}{4\pi R_1\gamma }\left[\frac{11}{6}-\frac{2R_2}{R_1+R_2}+\ln \left(\frac{F{R}_2}{4\pi R_1\gamma (R_1+R_2)}\right)\right],$$

$$1-\frac{{\Delta }_{12}}{2R_0}=\frac{F}{4\pi \gamma R_0}\left(\frac{5}{6}+\ln \left(\frac{F}{8\pi \gamma R_0}\right)\right).$$

$$x_{ij}(F)=\frac{F}{4\pi R_i\gamma }\left[\frac{11}{6}-\frac{2R_j}{R_i+R_j}+\ln \left(\frac{F{R}_j}{4\pi R_i\gamma (R_i+R_j)}\right)\right]+\sum _{k\ne j}G(\Delta {\theta }_{jk})\frac{F_{ik}}{R_i\gamma }.$$

Full text

1 Introduction

The theory of Morse and Witten [1, 2] yields formulae relating forces, distortions, and energies of a bubble (or droplet in the analogous case of an emulsion), under the action of forces due to contacts with walls or other bubbles. It proceeds from the case of a single bubble, pressed against a wall by buoyancy. An extension to the case of multiple contacts (and hence a foam), also in static equilibrium, was indicated in the original paper [1]. However, this has never been fully developed; for example, it does not account for polydispersity. Hence, while there have been some limited trials of the theory [3, 4] they have been restricted to monodisperse, or near monodisperse systems. In the present paper we take some steps towards a full implementation of the theory of Morse and Witten allowing for an arbitrary degree of polydispersity.

The theory reduces a foam (or emulsion) to a set of representative points (the centres of mass of the bubbles) with central forces between them, depending on their separation but not as simple local relations. Höhler and Weaire [2] have provided a review of the Morse–Witten theory, to which reference may be made for a more detailed understanding if necessary.

The present paper deals mainly with the case of a 2D foam, for which Weaire et al.[5] have developed theory analogous to the original 3D case; this is the starting point for the present work.

The 2D foam, while not completely realised in practical systems (such as that of bubbles trapped between two plates) is a familiar test ground for the theory of foams [6]. Generally similar to 3D foam, in terms of its properties, it is simpler in many respects, and more readily simulated and visualised. We anticipate analogous methods and results for the 3D foam, albeit with some important differences in detail, and a greater challenge to practical simulation.

Relatively dry (less than 10% liquid) 2D foam has been successfully simulated in the past with the Plat software [7, 8, 9, 10, 11, 12]. It is not based on an energy minimisation routine, but instead directly implements local equilibrium for a wet 2D foam. It models the films and liquid-gas interfaces as circular arcs, constrained to meet smoothly at vertices. This makes it quite an accurate model of 2D foam; however, the software suffers from a failure to converge for liquid fractions close to the wet limit. Therefore, we seek a method for 2D foams that is successful in that limit.

2 Morse–Witten theory in two dimensions

2.1 Basic Theory

In the primitive version of the 2D theory, a 2D bubble is pressed against a fixed line by a buoyancy force [5]. Just as in the 3D case, the distortion of the bubble shape from circular may be found in approximate analytic form, by solving a linearised Young–Laplace equation. This solution can be used to build up a description of the foam of many bubbles, and the forces between them.

It is expressed in terms of the radius $\rho(\theta)$, whose deviation from the unperturbed value $R_{0}$ is $\delta R(\theta)$,

[TABLE]

where $\theta$ is a polar angle relative to the point of contact.

The solution of the linearised Young–Laplace equation then results in [5]

[TABLE]

Here $F$ is the magnitude of the total contact force, $\gamma$ is the line tension, and $g(\theta)$ encapsulated the deformation of a bubble in response to a force as in Weaire et al.[5]. In the following we will often use the dimensionless force $f=F/\gamma$.

Equation (2) represents the deformation of the bubble in such a way that its centroid (or centre of mass), which represents its location, is kept fixed. The profile $\rho(\theta)$ (Equation (1)) is shown in Figure 1; it may be considered to represent a 2D bubble under the action of a point force $F$ at $\theta=0$, but can be used more generally. Real bubbles would not support such a singular deformation. Nevertheless, $g(\theta)$ can be used to predict the shape of a bubble subjected to realistic force distributions.

If this model is used to describe a contact with a straight line, analogous to a flat hydrophobic wall in 3D, then only part of this function is used, the profile being “capped” by a straight line [2]. This is the only case considered (in 3D) by Morse and Witten: hence the previous restriction to monodisperse foams. When describing polydisperse foams we require to deal with contacts with a curved boundary, appropriate to contacts between bubbles of different size (and pressure).

The reader unfamiliar with this subject may wonder why a body force (which we call buoyancy) has been introduced, while it has no place in the problem posed (a foam in the absence of gravity). In fact, the solution for a bubble under the action of several forces in equilibrium may be developed as a combination of the solution given here for the contacts of each bubble, with the effect of body forces cancelling out [2].

2.2 Contact between two bubbles of different sizes

Here we provide a generalisation of the Morse–Witten method to account for contacts between 2D bubbles of different sizes. We require to find the relation between the contact force and the deformation of each bubble, represented by $x_{i}$, i.e. the distance along the centre–centre line from the undeformed bubble to the contact point (see the inset of Figure 2). To lowest order, $x_{i}$ is the distance that the point at the cusp of the contact is displaced, which is

[TABLE]

from Equation (2). This is indicated in red in the inset of Figure 2. However, this overestimates the deformation at the contact (see Figure 3).

A simple derivation of the required relation between $F$ and $x_{i}$ follows. As with many other aspects of the theory, this deals with lowest-order expressions only, and can be developed most expeditiously by using these from the outset (and verifying by a more cautious method if necessary). Thus we can take for the force between two bubbles, to lowest order,

[TABLE]

where $2l$ is the width of the contact (Figure 2) and $p_{0}$ is the mean of the two (lowest order) bubble pressures $p_{i}=\frac{\gamma}{R_{i}}$. Hence

[TABLE]

where $2\theta_{i}$ is the opening angle of the contact, for $i=1,2$. We can also express $x_{i}$ in terms of $\theta_{i}$, as

[TABLE]

This improved expression for the deformation of bubble $i$ (= 1 or 2) is then expanded to $\mathcal{O}(f^{2})$ to give

[TABLE]

which is indicated in the inset of Figure 2. The relative deformation $x_{i}/R_{i}$ is the same for each of the two bubbles. The centre–centre distance $\Delta_{12}$ is then given by

[TABLE]

For two bubbles with radii $R_{1}=R_{0}+\Delta R$ and $R_{2}=R_{0}-\Delta R$, this results in the dimensionless change in separation as

[TABLE]

Thus, terms of order $f^{2}$ or higher need to be considered in the expansion of Equation (6) to account for polydispersity, This is in contrast to the situation in 3D, see Section 6. (Note that, in the monodisperse case, i.e. $\Delta R=0$, the correction term in Equation (9) reduces to $f^{2}/8$, not to $f^{2}/4$, as erroneously stated in the appendix of Weaire et al.[5]. This had no consequences for the results presented in that paper.)

In order to test the accuracy of Equation (7) we proceed as follows. For a given force, $f$, we draw two overlapping bubbles with facing contacts, using Equations (1) and (2). The centres of these are moved apart until their area of overlap is zero, giving the separation for that force.

For the range of normalised force shown ($0<f<0.5$), Equation (9) produces a relative error $<2\%$ (the relative error when considering only its linear part is up to $25\%$, see Figure 3).

Equation (7) may also be used to describe the case of a bubble $i$ in contact with multiple bubbles (of different radii). This results in a set of deformations at its contacts with its neighbours, $j$. The deformation $x_{ij}$ of bubble $i$ due to its contact with bubble $j$, is determined by the sum over all the contacts of bubble $i$,

[TABLE]

Here $\Delta\theta_{jk}$ is the angle the between the centre-centre lines of bubbles $i$ to $j$ and bubbles $i$ to $k$, where $k$ enumerates all the contacts of bubble $i$ (including $j$). $F_{ij}$ is the force experienced by bubble $i$ at its contact with bubble $j$.

The equivalent expression for three dimensions is given in Section 6.

Note however that the linearised theory contains errors of order $f^{2}$ from the outset which we do not claim to eliminate. Given that, the theory is surprisingly successful in improving the lowest order estimate. The situation is similar to that which was encountered in the application of Morse–Witten theory to the pendant drop, although different in detail [13].

3 Formulation of the Morse–Witten model

3.1 Description of a foam

We will proceed to apply Equation (10) to find an equilibrium structure of a polydisperse foam, in a numerical simulation. We consider $N$ bubbles in equilibrium in a square box with periodic boundary conditions. The bubbles are represented by their centroid positions ($\mathbf{c}_{i}$) and radii $R_{i}$. A contact between bubbles $i$ and $j$ has an associated contact force of magnitude $F_{ij}$.

This nonlinear problem is naturally approached by iterative methods. While its defining equations are simple, its implementation is challenging, because of the role of the contact network, which needs to be continually monitored and updated, as explained below.

3.2 Defining Equations

We seek an equilibrium configuration which satisfies the conditions A-D below where the variables to be yielded by iteration are

• •

the centre positions $\mathbf{c}_{i}$,

• •

the contact force magnitudes $F_{ij}$,

• •

and the contact deformations $x_{ij}$.

A) Force-deformation relation.

Forces and deformations must be consistent, that is, satisfy Equation (10).

B) Deformation-displacement relations.

For each contact, the separation of centres of mass, located at positions $\mathbf{c}_{i}$ and $\mathbf{c}_{j}$ must be consistent with the deformations $x_{ij}$ and $x_{ji}$, according to

[TABLE]

C) Action-reaction.

[TABLE]

D) Equilibrium of forces.

The vector net forces on each bubble $i$ must satisfy

[TABLE]

3.3 The contact network

As the system approaches equilibrium the shapes and positions of bubbles change. The contact network is not finally determined until equilibrium is reached, consistent with the above conditions. It requires to be updated as the approach to equilibrium proceeds. Buzza and Cates [3] applied the Morse–Witten theory to the case of an emulsion where the drops are arranged on a simple cubic lattice, for which this difficulty does not arise. Höhler and Cohen-Addad [4], while including a slight polydispersity, also used crystalline systems in which contact changes were excluded. For the disordered foams discussed here a new methodology is thus needed to deal with bubble rearrangements (topological changes).

4 Implementation of the Morse–Witten model

4.1 Iterative scheme

We have developed a practical iterative scheme that can produce an equilibrium structure satisfying the conditions of Section 3. Separate steps of iteration are designed to bring the configuration closer to satisfaction of the conditions. To start a set of bubble centres and radii is required. These can be obtained from other software, such as Plat [10], Bubble model [14], or Surface Evolver [15].

A) Force-deformation relation.

Given a configuration and deformations of each bubble contact, the corresponding forces are found by solving Equation (10) for $F_{ik}$, for each bubble in turn. This is a nonlinear equation, hence we solve it iteratively. This difficulty is also to be found in the work of Höhler and Cohen-Addad [4], and we adopt the same method as was used by them. That is, in each iteration $n$, the forces from the previous iteration are inserted in the quadratic term, leaving a linear equation to be solved. Additionally, we apply some damping to this procedure, implemented as

[TABLE]

where we have found $a=0.9$ to be a good choice. This helps to prevent oscillations in the forces, without slowing convergence too much.

B) Deformation-displacement relations.

The deformations are updated by

[TABLE]

in order to satisfy Equation (11).

C) Action-reaction.

$F_{ij}$ and $F_{ji}$ are replaced by their average.

D) Equilibrium of forces.

Each bubble located at position $\mathbf{c}_{i}$ is moved in the direction of the net force acting on it, according to

[TABLE]

where $b=0.1R_{0}/\gamma$. As convergence speed is directly proportional to $b$, we have selected as large a $b$ as possible for which the algorithm still converges.

The flowchart of the iteration is shown in Figure 4. Note that it contains additional steps in which the contact network is, if necessary, altered.

4.2 Updating the contact network

Negative forces

A negative force indicates a spurious contact, i.e. a contact which has arisen from an overlap while the system is out of equilibrium, and this is removed from the contact network. In practice, at most one negative force is eliminated for each bubble in a given iteration to provide stability of the algorithm. This is performed after updating the forces.

Overlapping bubbles

After moving the bubble positions, nominally non-contacting bubbles may overlap with each other. To detect this we calculate

[TABLE]

For $x_{ij}>0$ bubbles $i$ and $j$ overlap. This requires an update of the contact network, which is performed before updating the forces.

4.3 Convergence

The algorithm is terminated when the foam being simulated is close to equilibrium, satisfying all of the above requirements. This is determined numerically by calculating the net force on each bubble in the foam using the left hand side of Equation (13). We deem this equation to be satisfied for all bubbles if the largest net force encountered is less than $\gamma\times 10^{-4}$.

In this case the centroid positions given by the recurrence relation, Equation (16), will have converged, leading also to a convergence of the deformations, Equation (15). Thus, solving the deformation-force relationship Equation (10) repeatedly will produce the same set of contact forces each time and all the defining equations will be satisfied.

5 Tests and Typical Results

We have run tests of the above scheme for systems of up to 200 bubbles (the run time of the program scales quadratically with the number of bubbles), in a square box with periodic boundary conditions. The computations converged satisfactorily for liquid fraction exceeding around $\phi=0.12$, a liquid fraction where at least 80% of the Plat simulations fail [11]. We have not identified the reason for non-convergence beyond that point, but it is hardly surprising in a nonlinear problem of this kind, and may be rectified in due course. In order to validate the method, we have compared it with simulations using the Plat software, as introduced in Section 1.

To begin, a system of ten bubbles with a polydispersity of $\sqrt{\langle R^{2}\rangle/\langle R\rangle^{2}-1}\simeq 0.12$ was generated using the Plat software and the liquid fraction $\phi$ increased until a hard disk packing was achieved (Figure 5(a), top). A hard disk packing corresponds to a foam in the wet limit (at $\phi=\phi_{c}$) where all of the degrees of freedom are exactly taken up by the contacts between bubbles, and there are no additional constraints. In this case the average number of contacts is $Z_{c}=4(1-1/N)=3.6$ [16]. Ten bubbles constitutes a small enough system that, despite the general failure of Plat to converge in the wet limit, the cost of repeating simulations until it is found to be successful is sufficiently small so as to make it feasible. The centre positions of the bubbles were extracted and used to create a Morse–Witten simulation of the same system. The liquid fraction of both simulations was then decremented in parallel, down to a liquid fraction of approximately 0.12. The Morse–Witten simulation produced almost the same contact changes as the Plat simulation, although at values of $\phi$ shifted by roughly $\Delta\phi\simeq 0.01$ higher. In looking at this comparison, it should be borne in mind that the Morse–Witten formalism is inherently approximate.

We next consider the excess energy of a Morse–Witten foam, defined (in dimensionless form) by

[TABLE]

where $j$ enumerates the contacts of bubble $i$.

For ordered monodisperse foam, Princen calculated $\varepsilon(\phi)$ exactly [17, 18]. This presents a good test for the Morse–Witten model, which can be solved exactly in this case. Figure 6(a) shows excellent agreement between Princen’s exact result and the analytic solution of the Morse–Witten model in the wet limit ($\Delta\phi<0.02$). Our numerical simulation results match the analytic solution of the Morse–Witten model.

Also shown is a simple approximate solution of the Morse–Witten model, which can be obtained as follows. The energy per contact is given by elementary methods as $0.5F\delta R/(\gamma R_{0})$ and using the relation

[TABLE]

from Weaire et al.[5], along with the affine compression relation $\delta R/R_{0}=\Delta\phi/2(1-\phi)$, we obtain

[TABLE]

where $\phi_{h}=\pi/2\sqrt{3}$ is the critical packing fraction for a hexagonal disk arrangement. This relation (shown in Figure 6(a)) is in excellent agreement with the result of Princen for $\Delta\phi<0.015$.

In order to study the variation of excess energy $\varepsilon$ as a function of liquid fraction of polydisperse foams, 1000 foams of 100 bubbles each were prepared with an average polydispersity of $0.21\pm 0.02$. These simulations were run for a range of liquid fraction from $0.18$ to $0.12$ in steps of $0.001$. They were started deliberately higher than the expected value of $\phi_{c}\simeq 0.16$ so that the transition from unjammed collection of disks to jammed foams will not be missed. The critical value $\phi_{c}$, marks the onset of the excess energy.

Our simulations show that, similar to results from Plat [11], close to $\phi_{c}$, the energy varies roughly quadratically with the distance $\Delta\phi=\phi_{c}-\phi$ from $\phi_{c}$. Therefore, the values for $\phi_{c}$ were calculated individually for each 100 bubble system by fitting a straight line to the lowest eight points of the square root of the energy curve that were above $10^{-4}$. The average value obtained from this procedure is $0.843\pm 0.003$, consistent with previously published values for $\phi_{c}$ [9, 16, 19, 20, 14, 21]. The energy curves for these simulations were shifted by their respective $\phi_{c}$ values, and then averaged with a bin width of $0.001$ in $\Delta\phi$ to smooth the data. Figure 6(b) shows that, based on our 1000 simulations, $\varepsilon(\Delta\phi)\propto\Delta\phi^{2.2}$.

A further quantity of interest in the context of random packings is the variation of the average coordination number, $Z$, with liquid fraction. A log-log plot of our data (Figure 7) reveals a scaling of $Z-Z_{c}=\Delta Z\propto\Delta\phi^{0.52}$, consistent with results for packings using the soft disk model [22]. Such a scaling was recently disputed based on extensive computer simulations with Plat which resulted in $\Delta Z\propto\Delta\phi$, and it was argued that this was due to the deformability of soft bubbles [16]. The results presented here appear to put some doubts on this argument. Further simulations with Plat and the Surface Evolver software [15] (currently restricted to finite contact angles in two dimensions [18]) would be required to determine whether the reported linear scaling with $\Delta\phi$ might be due to some inherent bubble-bubble attraction that arises from the algorithms.

In the study of granular matter it is common to compute the contact force network [23, 24]. Granular packings are characterised by a very slow decay of the distribution of forces greater than the mean. Whether this is exponential or faster than exponential depends on the details of the simulations/experiments, such as dimensionality, solid friction, and partial size distribution [25, 26, 27].

In Figure 8(a) we show the contact force network for an equilibrated Morse–Witten foam of 100 bubbles at a liquid fraction of $\phi=0.13$. The width of each line in the contact network is proportional to the force magnitude. In addition, the bubbles are shaded according to their individual excess energies. Also shown in Figure 8 is a preliminary normalised distribution of contact forces. This is broadly similar to that found by Höhler and Cohen-Addad [4], however, further simulations are required to analyse its shape.

6 Extension to a 3D foam

The methodology developed above for the simulation of a 2D foam based on the Morse–Witten model lends itself to application also for 3D. As in 2D the foam will be represented by the centroid of all bubbles and a network of contacts. In 3D, the profile is expressed analogously to Equation (1) and Equation (2) becomes

[TABLE]

where

[TABLE]

from [1]. The expression for the deformation of bubble $i$, equivalent to Equation (10) and derivable in the same way, is given by

[TABLE]

to lowest order in $F$. The relative change in separation (equivalent to Equation (9)) between two bubbles where $R_{1}=R_{0}+\Delta R$ and $R_{2}=R_{0}-\Delta R$ is

[TABLE]

Again, symmetry tells us not to expect any terms of odd orders of $\Delta R$ in the separation. Equation (23) would need to be expanded to order $F^{2}$ to give terms of $\Delta R^{2}$. Taking for example $\Delta R=0.1R_{0}$, $F=0.5\gamma R_{0}$, the relative error that would result from using a formula for monodisperse foam would be of order $10^{-4}$. This explains the success by Höhler and Cohen-Addad [4] in using an expression derived for the monodisperse case in treating a slightly polydisperse case.

In order to model a 3D foam, an equivalent to Equation (10) is required. This is obtained by adding a non-local term to Equation (23) (see [4]) giving

[TABLE]

Similar to the procedure of Section 2.2, we determined the separation of two 3D Morse–Witten bubbles at their point of contact, for a given force $F$. We find that Equation (24) is reasonably accurate up to $F/(R\gamma)\sim 0.5$ (corresponding to the dry limit) for low polydispersity, and $F/(R\gamma)\sim 0.05$ (corresponding to $\phi\sim 0.24$) for high polydispersity. The appearance of a non-linear $F\ln(F)$ term makes the 3D case somewhat different from the 2D one presented here: nevertheless we hope that further improvement of the 2D methods will assist in the greater computational task of implementation in 3D.

7 Conclusion

We have shown how polydispersity can be accommodated in the Morse–Witten theory, in such a way as to give satisfactory results for a typical disordered polydisperse foam that is close to the wet limit. The extension of the theory to 3D is quite natural, although the implementation becomes conceptually more difficult to visualise and check, and there is an obvious increase in computational demands. The transparency of the theory and its direct relation to a force network (Figure 8) is attractive. However, it should be noted that it has proven a computational challenge that was hardly anticipated, and is worthy of further attention.

In the polydisperse foam the bubble-bubble interfaces have pronounced curvature: this is accounted for in the present formulation, being related to differences in bubble sizes. One might well ask what is the case in a monodisperse disordered foam? (Despite some doubts in the past, this can indeed exist, even in 2D). Since the bubbles are not equivalent, surely their pressures are slightly different, hence the interfaces are curved? This is correct in principle, but the effect is surely very small, and of higher order in the forces than what is considered here.

8 Acknowledgements

DW wishes to acknowledge discussions with Reinhard Höhler and Tom Witten.

Research supported in part by a research grant from Science Foundation Ireland (SFI) under grant number 13/IA/1926 and from an Irish Research Council Postgraduate Scholarship (project ID GOIPG/2015/1998). We also acknowledge the COST action MP1305 ‘Flowing matter’ and the European Space Agency ESA MAP Metalfoam (AO-99-075) and Soft Matter Dynamics (contract: 4000115113).