Ideal wet two-dimensional foams and emulsions with finite contact angle

We present simulations that show that an ideal two-dimensional foam with a finite contact angle develops an inhomogeneity for high liquid fraction $\phi$. In liquid-liquid emulsions this inhomogeneity is known as flocculation. In the case of an ordered foam this requires a perturbation, but in a disordered foam inhomogeneity grows steadily and spontaneously with $\phi$, as demonstrated in our simulations performed with the Surface Evolver.


Introduction
In emulsions, the term flocculation refers to the (spontaneous) clustering of droplets, leading to the formation of density inhomogeneities.Here we describe the onset of flocculation in computer simulations of two-dimensional (2D) liquid foams.This only occurs in systems where the liquid-gas interfaces meet at a finite contact angle θ, as illustrated in figure 1(b).Previous simulations of 2D foams with finite liquid fraction φ have taken this contact angle to be zero [1] or, where that was not feasible for numerical reasons, as small as possible [2].We treat both liquid and gas as incompressible, and so the results apply equally to emulsions.
Two-dimensional foams have properties that are broadly similar to their three-dimensional counterparts, but are much simpler to analyse.Their study has early antecedents [3], and continues to be of interest today.The usual theoretical model is entirely two-dimensional, while the third dimension may be significant in relevant experimental systems, such as a foam trapped between two plates [3,4].
A simulated example of an ideal 2D foam with finite liquid fraction φ and finite contact angle θ is shown in figure 1(a).It was produced by the method described in Appendix A. The gas bubbles are surrounded by a network of smoothly-curved thin films connecting the liquid-filled Plateau borders, which each have three or more sides.There are therefore two types of interface: liquid-gas interfaces around each Plateau border, and gas-liquid-gas interfaces forming the bubble-bubble contacts.On each of the interfaces the Laplace-Young law relates the product of interfacial tension and curvature to the pressure difference across the interface [5,6], and consequently liquid films are represented by arcs of circles that meet at the vertices of Plateau borders.The films are considered to be infinitesimally thin, so all of the liquid in the foam is considered to be contained in the Plateau borders.A finite contact angle implies that the interfacial tension associated with the bubble-bubble interfaces is less than twice that associated with the Plateau borders (see figure 1(b)).Note that throughout this paper "interfacial tension" is used for what is really a line tension (or energy per unit line length) in such an idealised 2D model.The contact angle θ is given by where γ f is the interfacial tension of each side of a liquid film; this is equal to or smaller than the (bulk) interfacial tension γ associated with the gas-liquid interfaces of a Plateau border [10].The surface energy E(φ, θ) of the foam is the sum of the lengths of all interfaces, multiplied by their appropriate interfacial tension.
For bubbles immersed in a liquid, the presence of finite contact angles entails net attractive forces between them (see figure 2) when they are only slightly compressed together.This is similar to the attraction between droplets, when considering emulsions [7].
In the present paper we address some basic consequences of introducing finite contact angles into the standard model of 2D foams, by analysing simulations carried out with the Surface Evolver software of Ken Brakke [8].We will see that finite contact angles, even if apparently very small, can have large effects.
Henry Princen introduced the concept of a contact angle between a thin liquid film and its adjacent Plateau border [9] and he established its physical significance in foams and emulsions by analysing 2D ordered (hexagonal) monodisperse structures [10], which admit analytical solutions.The Surface Evolver and current computational resources enable the simulation of disordered foams, which is the usual practical case of interest.Princen's work was stimulated by his own measurements of contact angles for soap films in contact with bulk solution [11].He found that finite contact angles up to 17 • could be achieved in surfactant (SDS) solutions at sufficiently high concentrations of added electrolytes.
We shall begin by recapitulating Princen's 2D model, which gave a number of exact results.This turns out to be instructive when discussing our own results from simulations of disordered polydisperse foams.

Ordered hexagonal foams
Even the apparently trivial case of an ordered 2D foam proves to present some challenges to detailed understanding, so we shall examine it carefully.
Monodisperse 2D bubbles are arranged in equilibrium on an hexagonal lattice (figure 3).Princen [10] showed that the liquid fraction φ θ 0 of such a packing at zero compression, which corresponds to zero osmotic pressure [10,12,13]) and minimal surface energy [10], is given by In the case of a zero contact angle (θ = 0), φ θ 0 reduces to the familiar value which is the liquid fraction of an hexagonal close packing of circular bubbles; it is known as the wet limit of an ordered monodisperse foam with zero contact angle.However, the bubbles are never circular when the contact angle is finite.Figure 3 shows two examples of an ordered foam with finite contact angle, for different values of the liquid fraction.
Princen calculated the work per unit area, ∆W θ , required to compress (by removing liquid) a foam from the liquid fraction φ θ 0 to any given φ [10].The bubble area is written in terms of the radius R of an undeformed circular bubble of the same area, resulting in 1 where the normalizing factor γ/R is often called the Laplace pressure.
1 Princen [10] describes his calculation in terms of the deformation of columns of hypothetical cylindrical emulsion drops.We have re-written his expression in terms of liquid fraction φ, rather than gas fraction (1 − φ), and corrected one misprint (missing superscript θ in his equation ( 35)).2)), the maximum value of liquid fraction at which the bubbles still remain in contact (eqn.( 6)).
In the following we will consider the energy E(φ, θ) per bubble as a function of liquid fraction φ and contact angle θ.Using eqn.( 4) and In the dry limit (φ = 0, i.e., close-packed hexagons), this reduces to The energy E(φ, θ) in eqn.( 5) may be characterised by two different critical values of the liquid fraction, as shown in Figure 4.These are • the liquid fraction φ θ 0 (eqn.( 2)) at which the energy has a minimum; • the maximum value of liquid fraction φ θ m at which the bubbles remain in contact, given by with ϕ = π/3 − 2θ.This is the hypothetical wet limit of an ordered monodisperse foam with finite contact angle.
Figure 5 displays the variation of φ θ 0 and φ θ m with contact angle θ.Note that θ = 0 is a special case because the critical values coincide, 3) 0.093.Figure 3 shows an ordered foam at each of the critical liquid fractions φ θ m and φ θ 0 , for contact angle θ = 6.26 • .When the contact angle is finite, the non-circular bubbles at φ θ m contact their neighbours at a point.However, this situation differs from the conventional wet limit for zero contact angle, as we explore below.All points on the curves for energy as a function of liquid fraction, shown in figure 4, correspond to (possibly metastable) equilibrium structures.Since the liquid fraction φ θ 0 corresponds to the minimum energy, any homogeneous structure for which φ > φ θ 0 must be metastable, at least for an infinite sample.This is because an inhomogeneous structure can be defined, with constant energy (in the limit of infinite   sample size) close to E(φ θ 0 ).This scenario is illustrated in figure 6, where we have represented the energy of the inhomogeneous structure by a horizontal line beyond φ θ 0 .The formation of inhomogeneities in the ordered monodisperse structure can be simulated by considering a representative area of foam that contains a large number of bubbles, and increasing the liquid fraction by expanding the system while keeping the bubble areas fixed.An example of such an inhomogeneous structure is shown in figure 7. Finite contact angles were included by assigning different values for the line tension in the bubble-bubble interfaces (2γ f ) and the Plateau border sides (γ); the contact angle is given by eqn.(1).
(For further details see Appendix A.) The inhomogeneity does not arise spontaneously in these simulations, but requires a perturbation at a liquid fraction just above φ θ 0 .This is achieved by randomly displacing all of the Plateau border vertices by a small distance (less than a typical Plateau border width); this allows the system to escape from the metastable branch by undergoing topological transitions that are triggered when film lengths go to zero.These result in the formation of cracks, or liquid pools, as illustrated in figure 7. The requirement of a perturbation to trigger instability is a familiar feature of highly symmetric structures that are locally stable, for which alternative structures of lower energy are available.Subsequently, increasing or decreasing the liquid fraction keeps the foam at roughly constant energy.Small fluctuations in the energy occur because, as the liquid fraction changes in this finite sample of foam, there are short periods during which the energy increases elastically, followed by topological transitions that reduce the energy.As the sample size increases, such fluctuations become less marked.This instability is reminiscent of the observations by Abd el Kadar and Earnshaw [14] in experiments with monodisperse 2D bubble rafts in an hexagonal confinement.About 25 minutes after foam formation, cracks appeared within the perfectly ordered structure, rapidly leading to the formation of a hole in the monolayer.
No bubbles are lost in this process.The authors attribute the development of this inhomogeneity to the movement and build-up of local stresses, for example due to small variations in bubble size.It is conceivable that the existence of a finite contact angle plays a further role, although the concept cannot be applied straightforwardly to bubble rafts, as opposed to 2D foam trapped between plates.

Disordered foams with finite contact angle
Having established the effect of a finite contact angle for ordered foams we now turn to the results of Surface Evolver simulations of disordered foams, as exemplified in figure 8.The results presented below are based on five samples with N = 1500 bubbles and different values of the contact angle θ between 2.6 • and 15.9 • .The foams are polydisperse and the variation in bubble areas differs between samples.Details of the simulations are presented in Appendix A.
We again find that inhomogeneities develop as the liquid fraction φ is increased.However, these inhomogeneities do not necessarily result from metastable structures of the kind discussed in the previous section.
The disordered samples spontaneously undergo discrete topological transitions as the liquid fraction is increased; these occur when the length of a film between two neighbouring bubbles goes to zero, causing them to separate.In the simplest case, two bubbles that share a shrinking edge separate, and two three-sided Plateau borders merge to form a four-sided Plateau border, as shown in figure A.11. Another possibility occurs when non-adjacent edges on a Plateau border with four or more sides come into contact and a new edge is formed.Various combinations of these topological transitions are responsible for the highly irregular liquid regions shown in figure 8.
Figure 9 illustrates the growth of inhomogeneity as the liquid fraction is increased.The average area of the Plateau borders A pb is given by where N is the constant number of bubbles, A is the constant average bubble area, and N pb is the number of Plateau borders, which decreases when they merge.If the Plateau borders do not merge, i.e. if they all remain three-sided, the pre-factor 2N/N pb is unity.The average Plateau border area, scaled by the value of A pb that would result if they didn't merge, is shown in the inset to figure 9. Note that the average Plateau border area rises more quickly at lower contact angles; this occurs because the tendency to combine is more pronounced when the contact angle is smaller, so the film lengths are smaller, and consequently, topological transitions are more frequent.The evolving inhomogeneity exhibits an acceleration in the growth of a few, increasingly larger, Plateau borders, while most remain small.The area of the largest Plateau border, A max pb , within a sample (figure 9) increases roughly exponentially with liquid fraction.9: The development of inhomogeneities in disordered foams with finite contact angle is illustrated by the observation that, while the average area of a Plateau border A pb grows sublinearly with liquid fraction φ (inset), the area of the largest Plateau border A max pb (at each value of φ) grows approximately exponentially (note the logarithmic vertical axis).The data shown is for five simulations of disordered foams, with values of contact angle as indicated.The areas are normalized by the average area A pb that would result if Plateau borders did not merge, eqn.(7).Note that the drops in area at low liquid fraction correspond to the splitting of a Plateau border with more than three sides into three-sided Plateau borders, as explained in Appendix A.
(a)  8, is accompanied by a decrease in energy as a function of φ followed by saturation to a roughly constant value.The energy has been rescaled using the polydispersity parameter p defined in eqn.(B.1).

Appendix A. Simulations using the Surface Evolver
We carried out simulations of both ordered (hexagonal) and disordered (polydisperse) foams with finite contact angles using Ken Brakke's Surface Evolver software [8].
In principle, all the properties of an ordered hexagonal foam can be captured from a single cell, considered with periodic boundary conditions.However, since we seek an instability of this structure to an inhomogeneous state, we must instead simulate many cells.We chose to reproduce the basic hexagonal cell 256 times to form a 16x16 hexagonal lattice with periodic boundary conditions.
The disordered foams are made as dry foams, also with periodic boundary conditions, from a Voronoi construction in the usual way [22]; these have N = 1500 bubbles, with average area A close to 1, and different polydispersity.
Both ordered and disordered foams are turned into wet foams of liquid fraction φ ≈ 0.03 by adding small triangular Plateau borders at each three-fold vertex; the liquid fraction is set by the total area of all Plateau borders.They have no individual area constraints, and therefore all have the same pressure.
Each side of each Plateau border is associated with a fixed bulk interfacial tension γ ≥ 1 2 ; the interfacial tension in the thin liquid films is set to 2γ f = 1.The contact angle is then given by eqn.(1).The current version of the Surface Evolver software does not allow for the simulation of 2D foams with zero contact angle and we find that θ = 2.6 • is the minimum value that we can choose.Examination of the consequences of this limitation provided some initial motivation for this study.
In the Surface Evolver, each edge is represented as a circular arc, and a local minimum of the interfacial energy is sought using up to 2 × 10 4 iterations to achieve a relative accuracy close to 10 −6 .
The liquid fraction is changed by increasing the area of liquid and keeping the area of gas constant (so the size of the periodic box increases).Small increments in φ allows us to explore a large range of liquid fractions: at each step, the liquid fraction is increased by 0.001, up to about φ = 0.25.We used the gradient descent method for energy minimisation, with occasional Hessian iterations, to move towards a minimum of interfacial energy j γ j L j .Here L j denotes the lengths of the interfaces, and γ j is either γ or 2γ f , depending on whether the interface is associated with a Plateau border or a film.
As the foam evolves, bubble rearrangements (topological changes) occur.We chose a critical film length of 10 −4 below which these are triggered.The difficult step is in recognising when a four-sided Plateau border should split into two three-sided Plateau borders (or, in general, a many-sided Plateau border should split into two parts).To achieve this we check whether two sides of any Plateau border with more than three sides overlap, as illustrated in figure A.11.If they do, these two sides are joined, to split the Plateau border.
We briefly increase γ and perform a few iterations, and then the minimization continues.A similar process is used for Plateau borders with more than four sides.Note that all Plateau borders have the same pressure, irrespective of their number of sides, and that their total area is conserved (rather than individual areas), so that no repartition of the liquid is required.

Figure 1 :
Figure 1: (a) Typical simulation of a two-dimensional foam with liquid fraction φ = 0.1 and contact angle θ = 9.5 • .The image shows part of an equilibrated foam of 1500 bubbles.(b) Close-up of a four-sided Plateau border, showing the contact angle θ.The vectors represent the forces, due to interfacial tension, acting on the point of contact. θ

Figure 2 :
Figure 2: The (line) energy of two isolated circular bubbles is reduced when they share a common interface, if the contact angle θ is finite.

Figure 6 :
Figure6: Variation of energy E(φ, θ) of an initially ordered hexagonal foam (for contact angle θ = 10 • ).For values of liquid fraction exceeding φ θ 0 (marked by a blue square), the hexagonal bubble arrangement is metastable, since alternative inhomogeneous structures exist with lower energy.(Figure7(a) shows such an example, obtained from simulation.)The light blue triangle at critical liquid fraction φ θ m corresponds to the point where bubbles in an hexagonal arrangement no longer touch.

Figure 7 :
Figure 7: (a) Example of a Surface Evolver simulation of a monodisperse foam (φ = 0.12, contact angle θ = 6.26 • , 256 bubbles) which shows the formation of inhomogeneities upon an increase in liquid fraction above φ θ 0 = 0.084.The appearance of inhomogeneities was triggered by a random perturbation of the ordered structure.(b) Upon a further increase in liquid fraction the energy remains roughly constant.

Figure 8 :
Figure 8: Surface Evolver simulation of a disordered foam of 1500 bubbles and contact angle θ = 2.6 • .An increase in liquid fraction leads to the appearance of cracks or incipient flocculation.The images shown correspond to close-ups of the foam at values of liquid fraction φ = 0.04, 0.08, 0.12, 0.16, 0.20 and 0.24 (from top left to bottom right).

Figure A. 11 :
Figure A.11: Numerical method for performing topological changes.(a) Two three-sided Plateau borders join when the interface between them shrinks beow a critical length.The four-sided Plateau border is created by deleting this short interface and merging the Plateau borders.(b) Two sides of a four-sided Plateau border approach each other; when they overlap, a new interface is inserted, creating two three-sided Plateau borders.A similar process is used for Plateau borders with more than four sides.Note that all Plateau borders have the same pressure, irrespective of their number of sides, and that their total area is conserved (rather than individual areas), so that no repartition of the liquid is required.
Variation of energy with liquid fraction φ for five disordered foams with finite contact angles.The gradual development of inhomogeneities, as shown in figure