3D simulations of wet foam coarsening evidence a self similar growth regime
In wet liquid foams, slow diffusion of gas through bubble walls changes bubble pressure, volume and wall curvature. Large bubbles grow at the expenses of smaller ones. The smaller the bubble, the faster it shrinks. As the number of bubbles in a given volume decreases in time, the average bubble size increases: i.e. the foam coarsens. During coarsening, bubbles also move relative to each other, changing bubble topology and shape, while liquid moves within the regions separating the bubbles. Analyzing the combined effects of these mechanisms requires examining a volume with enough bubbles to provide appropriate statistics throughout coarsening. Using a Cellular Potts model, we simulate these mechanisms during the evolution of three-dimensional foams with wetnesses of \xcf\x95 = 0.00, 0.05 and 0.20. We represent the liquid phase as an ensemble of many small fluid particles, which allows us to monitor liquid flow in the region between bubbles. The simulations begin with 2 \xc3\x97 105 bubbles for \xcf\x95 = 0.00 and 1.25 \xc3\x97 105 bubbles for \xcf\x95 = 0.05 and 0.20, allowing us to track the distribution functions for bubble size, topology and growth rate over two and a half decades of volume change. All simulations eventually reach a self-similar growth regime, with the distribution functions time independent and the number of bubbles decreasing with time as a power law whose exponent depends on the wetness.
1. Introduction<!>2.1. Foam structure<!>2.2. Foam dynamics<!>3. Results<!>4. Discussion and conclusions
<p>Liquid foams include soap froths and food and industrial foams, as well as many liquid-phase pollutants. Foams consist of gas bubbles surrounded by a continuous liquid phase. Bubble walls are very thin; their thickness strongly depends on the intensity of mutual repulsion between surfactant-loaded gas-liquid interfaces and only weakly on foam wetness (ϕ=liquid volume/total volume) below a critical wetness at which the foam transforms into a bubbly liquid. When three bubble walls meet, they form edges, known as Plateau borders, which are tube-like structures with a cross section which is roughly a triangle with its straight sides replaced by concave circular arcs. Wet foams and dry foams both have bubbles separated by thin walls, the difference being that wetter foams have larger Plateau borders [1, 2].</p><p>Slow diffusion of gas through bubble walls and Plateau borders, due to the Laplace pressure difference, changes bubble pressure, volume and wall curvature. Bubbles, Plateau borders and walls move and relax towards mechanical equilibrium, while liquid flows inside the Plateau borders and walls as well. Bubbles larger than the mean bubble size (or bubbles with more than the mean number of faces for dry foams) grow at the expenses of smaller (fewer-faced) ones. Smaller (fewer-faced) bubbles shrink more and more quickly, and eventually disappear. The number of bubbles decreases in time, and the average bubble size increases: the foam coarsens [1, 2, 3] (figure 1).</p><p>X-ray tomography experiments [5] and numerical simulations [6, 7] have found that dry foams reach a self-similar growth regime where the topology and relative bubble size distributions are steady while the mean bubble radius increases as t1/2, as expected from dimensional analysis [8]. Hence, the mean bubble volume increases as t3/2 and the number of bubbles decreases as t−3/2. In the opposite limit, in a bubbly liquid, gas flows from small bubbles to the liquid phase, where it diffuses, to eventually flow back into larger bubbles due to Laplace pressure difference, now between the bubbles' internal pressures and the gas saturation pressure in the liquid. This Ostwald ripening dynamics also generates a self-similar growth regime, where the relative bubble size distribution is constant, and the mean bubble radius grows as t1/3 [9]. Hence, the mean bubble volume grows as t and the number of bubbles decreases as t−1.</p><p>To investigate coarsening for foams of intermediate wetnesses, Isert, Maret and Aegerter [10] levitated an aqueous foam in a magnetic field: the water's diamagnetism counterbalanced gravity, so the wetness remained uniform throughout the foam. At long times, for foams with wetness less than ϕ = 0.25 the mean bubble radius grew as t1/2 as in a dry foam, for foams with wetness greater than ϕ = 0.35 the mean bubble radius grew as t1/3 as in a bubbly liquid, and the mean bubble radius growth exponent decreased from 1/2 to 1/3 between ϕ = 0.25 and ϕ = 0.35. Determining whether power law growth in bubble size at these wetnesses reflected self-similar growth with constant relative bubble size distributions would require determining individual bubble properties within the levitated foam. Such a measurement would be quite difficult, though it might be possible using e.g. X-ray tomography [4].</p><p>Numerical simulations provide access to properties of individual bubbles (volume, face-number and edge-number distributions, and their spatial correlations) with arbitrary spatial and temporal resolution. Complicating such simulations are the multiple mechanisms, length and time scales affecting coarsening. Slow volume changes in bubbles cause the gradual increase of Plateau border cross section. Both bubble and Plateau-border growth can induce topological transformations in bubble and Plateau-border adjacency, such as neighbors swapping (T1), bubble disappearance (T2), or Plateau-border fusion, which can cause further rapid cascades of rearrangement. These rearrangements, in turn, change the slower dynamics of gas exchange (between neighboring bubbles, and between bubbles and the liquid phase) and flow in Plateau borders. Length scales include the thickness of the thin walls (which can range from hundreds of nanometers to many microns), the size of the Plateau borders (which can range from hundreds of microns to millimeters), the lengths of Plateau borders and the radii of bubbles (which can range from sub-millimeter to centimeters) depending on the foam, wetness and degree of coarsening. The close feedback among all these mechanisms and scales makes simulation of wet foam dynamics challenging [11].</p><p>The Cellular Potts model (CPM) can simulate the structure of foams in mechanical equilibrium, with a resolution (number of pixels per bubble) which matches or exceeds that of experimental images in two or three dimensions. The CPM can also simulate the quasi-static dynamics of foam coarsening, as slow gas diffusion leads to a succession of mechanically equilibrated foam structures. The CPM correctly reproduces bubble rearrangements and disappearance.</p><p>The CPM developed from Potts model simulations of coarsening in metallic polycrystals [12], because the mass flux between bubbles and crystals both obey Fick's law (a local dynamics which produces the same growth dynamics as the Young-Laplace law for bubbles). Such CPM simulations agree in almost all respects with experimental observations of dry liquid foams [13]. However, in the CPM, mass flow between bubbles and surface-tension-driven relaxation of wall shape occur over similar time scales, while the former is much slower than the latter in almost all liquid foams. Extending the CPM method by interleaving steps in which mass is free to diffuse between bubbles with steps in which we constrain the bubble volumes to their current values but allow the bubble walls to relax allows us to reproduce the experimental hierarchy of time scales. The more relaxation steps per diffusion step, the faster the relaxation rate relative to the diffusion rate [6].</p><p>CPM simulations of 3D dry foam coarsening [6, 7] and 2D wet foam coarsening [14] both reach self-similar growth regimes. Here we compare 3D CPM simulations (with additional shape relaxation steps) of foams with wetnesses of ϕ = 0.00, 0.05 and 0.20. All three reach self-similar growth regimes with enough bubbles to yield statistically significant distribution functions.</p><!><p>Like experimental images, CPM 3D foam simulations represent space as a cubic grid of voxels, and each bubble as a connected set of voxels with particular properties. We employ a space of VT = 2003 voxels with periodic boundary conditions, initially filled with 2 × 105 bubbles for ϕ = 0.00 and 1.25 × 105 bubbles for ϕ = 0.05 and 0.20. Each voxel, located at a position r⃗, has a label S, corresponding to a bubble or portion of liquid. Defining the liquid phase as a single domain would allow liquid to move instantaneously within the foam, effectively giving it zero inertia and viscosity. To ensure local conservation of liquid volume and to keep the liquid velocity finite, we follow Fortuna et al. [14], and subdivide the liquid phase into numerous small fluid particle subdomains (see figure S1 in the supplementary materials).</p><p>The CPM energy function includes a surface energy and a volume constraint: 
(1)E=∑r→∑v→(r→)J(Sr→,Sv→)+∑SJ(S)(VS-VStarget)2 where the sum over v⃗(r⃗) extends through the 5th neighbors around the voxel at r⃗ (56 neighbors), to reduce lattice-anisotropy effects [15], and is only evaluated for neighboring voxels with distinct labels S. Figure 2 illustrates the grid neighborhood in the energy calculation. The second term in the r.h.s. of Eq. (1) constrains the volume of the labeled region S, VS, around its target volume, 
VStarget, when λ(S ), which is the compression modulus (inverse compressibility) of the material, is > 0.</p><p>While an "ideal" fluid particle would be negligibly small and perfectly incompressible, we assign our fluid particles a Vtarget = 7 voxels so they are much smaller than typical bubbles, but large enough to move freely (see figure S2 in the supplementary materials for the effect of varying the target volume of the fluid particles, which affects coarsening dynamics). The compression modulus λliq = 10 is not critical, since coarsening dynamics are robust under moderate variation of this parameter (see figure S3 in the supplementary materials). However, making λliq too large (i.e. the compressibility too small) leads to non-physical "freezing" of the liquid in the Plateau borders.</p><p>Simulations begin with a number of fluid particles randomly distributed over the grid. The foam wetness ϕ is simply the number of voxels in fluid particles (controlled by their number and target volume) divided by the lattice size. We fill the rest of the grid with bubbles with initial target volumes randomly drawn from a log-normal, normal or bidisperse probability distribution.</p><p>The first term of the r.h.s. of Eq. (1) defines the interfacial energy between the different types of regions:
(2)J(liq,liq)=0.10,J(liq,gas)=J(gas,liq)=1.00,J(gas,gas)=1.99.</p><p>J(liq, liq) would be zero in a real liquid, but a small positive value helps to maintain the fluid particles as connected domains without significantly affecting their dynamics. J(liq, gas) is the energy per unit area of the gas-liquid interface and J(gas, gas) is of the order of 2J(liq, gas) because a wall separating two bubbles consists of two gas-liquid interfaces. If J(gas, gas) were strictly equal to 2J(liq, gas), the fluid particles would spread in the walls as in a bubbly liquid. We set J(gas, gas) slightly lower than 2J(liq, gas) to account for the positive disjoining pressure (2J(liq, gas) − J(gas, gas)) between the soap films in the walls, which serves to keep the walls thin, as observed in experiment. The greater the disjoining pressure, the greater the fraction of liquid that localizes in the Plateau borders and the less that penetrates into the bubbles walls. As a result, we expect coarsening to be faster for higher disjoining pressures. For the range 2J(liq, gas) − J(gas, gas) = 0.1, 0.01 and 0.001, the simulation results are robust and the liquid's spatial distribution is essentially the same (see figure S4 in the supplementary materials).</p><p>The wall thickness and Plateau-border shape depend on the interaction range (v⃗(r⃗) in Eq. (1)) and the size of the fluid particles.</p><!><p>Configurations change stochastically in the CPM. We randomly pick a pair of first neighbor voxels with different labels and calculate the current energy, Ei, and the energy after copying the label S from the first to the second voxel, Ef. If ΔE = Ef − Ei ≤ 0 we always make the copy. If ΔE > 0 we make the copy with probability exp (−ΔE/T), where T is a Boltzmann-like fluctuation amplitude which we choose just large enough to overcome lattice pinning effects (the macroscopic fluctuation amplitude in a real foam would be close to 0, but real foams do not experience lattice pinning). The simulation unit which corresponds to time is the Monte Carlo Step (MCS): 1 MCS consists of N voxel-copy attempts, where N is the total number of voxels in the grid.</p><p>Like all models, the CPM's simplification of reality may introduce artifacts. The CPM models mass transfer as voxel transfer. The transfer of mass between two bubbles in contact across a thin wall obeys Fick's law in both CPM and experiment, so the key question is the rate of transfer as a function of parameters. In both experiment and simulation, we expect the transfer of mass across the Plateau borders to be slower than across thin walls. However, the dependence of the rate of transfer of bubble mass between bubbles via the Plateau borders depends less obviously both on model parameters and geometry and is not simple to separate from transfer across thin walls, so both the rate and functional form of the transfer might differ from experiment. Prior work has shown that 2D CPM simulations of wet foams using the current method reproduce the bulk growth behaviors and distribution functions seen in experiments, including their dependence on foam wetness [14]. While such agreement shows that the functional form of the mass transfer via Plateau borders must be similar in CPM and wet foams, the agreement does not demonstrate that the functional forms are identical, nor how the rate depends on model parameters. In particular, for foams composed of a gas with fairly low solubility in the liquid phase, we expect gas transport via walls to be much faster than via Plateau borders. In our current CPM simulations, the interaction range is a significant fraction of the typical Plateau-border size, thus we expect the transport rates via walls and Plateau borders to be more similar and the difference in growth rate between dry and wet foams to be less pronounced in the simulations than in an experiment with a foam with a low solubility gas phase. In the near future, we plan to examine the explicit relationship between Plateau border shape, size and model parameters, which will allow us to control the relative transport rates explicitly.</p><p>We first relax the starting configuration for 200 MCS, to allow the bubble and liquid volumes and geometries to adjust to generate a foam with the initially specified area distribution, then begin the simulation proper. To ensure that the equilibration of wall and Plateau-border shapes is much faster than the rate of mass diffusion, after each 20 MCS of coarsening steps we run at least 10 MCS of relaxation steps, continuing until the number of accepted voxel copies during a relaxation MCS is zero or less than 1% of the number of accepted voxel copies during the previous relaxation MCS. Because mass transfer controls the time-scale of coarsening in experiment, with relaxation occuring simultaneously, the simulation time t in our graphs reflects the number of coarsening steps only. During coarsening steps we set λgas = 0 so bubble volumes are free to vary and T = 8 for ϕ = 0, and T = 25 for ϕ> 0 since the volume constraint on the fluid particles increases their pinning to the grid. During relaxation steps, we set T = 0, λgas = 7, and 
VStarget for each bubble equal to the bubble's volume at the beginning of the step. The volume constraint prevents mass transfer between bubbles, while allowing boundaries to move to minimize their interface energy.</p><p>The simulations use the CompuCell3D environment [16], publicly available at http://www.compucell3d.org/. We ran the simulations on an Intel Core I7-3700K processor. We ran 10 replicas for each set of simulation parameters, starting from a log-normal volume distribution for each set of parameters, with ϕ = 0.00, 0.05, and 0.20. Each simulation replica took from 2 to 7 days to run, depending on the number of simulated objects (bubbles plus fluid particles), which varies with wetness. Because CompuCell3D limits the total number of objects in a simulation, we chose the size of the grid to respect that limit.</p><!><p>Figure 3 shows 2D sections of a simulated 3D foam before (left) and after (right) a set of relaxation steps. Relaxation smoothes the boundaries while preserving bubble volumes. Figure 4 shows 3D images of a larger simulation for ϕ = 0.20, run on a 3003 grid: boundaries are smooth and the foam texture closely resembles that in the experiments in Fig. 1.</p><p>To check whether the foam reached a self-similar growth regime we compared the bubble growth rate, G, as a function of the normalized bubble volume V/〈V〉 (figure 5) at different times. The growth rates are identical within error for t > 600 MCS for every wetness value (upper panel for ϕ = 0.20, data not shown for ϕ = 0.00, 0.05). The lower panel shows the self-similar growth rates (averaged over 10 runs and 31 times) for ϕ = 0.00, 0.05, and 0.20. Figure S5 in the supplementary materials shows a detail of this figure: the growth rates for ϕ = 0.20 significantly differ from those in the drier foam simulations.</p><p>Figures 6 and 7 show the distributions of normalized bubble volumes and face number. Both distributions are time invariant for t > 600 MCS for all wetnesses. The lower panels compare the average distributions for simulations with ϕ = 0.00, 0.05, and 0.20. The normalized volume distributions are similar, though statistically distinct (see figure S6 in the supplementary materials). The face-number distributions differ dramatically for foams with different wetnesses.</p><p>The upper panel in figure 8 shows log-log plots of the decrease in number of bubbles starting from initial configurations with a log-normal distribution of bubble volumes. The straight line behavior observed for all wetnesses for t > 600 MCS over roughly one decade is compatible with a power law, which is a requirement for a self-similar growth regime [8]. Fits for 600 MCS ≤ t ≤ 3500 MCS yield volume scaling exponents of 1.50, 1.33, and 1.20 ± 0.01, for ϕ = 0.00, 0.05 and 0.20 respectively (corresponding to radius scaling exponents of 0.5, 0.44 and 0.4). For t > 600 MCS, the number of bubble rearrangements NT1 per MCS also decreases as a power of time (figure 8, inset in the upper panel). The lower panel in figure 8 shows log-log plots of the decrease of the number of bubbles for ϕ = 0.20, starting from initial configurations with normal, log-normal and bidisperse bubble volume distributions. The growth exponent is independent of the initial volume distribution. Figure S7 in the supplementary materials shows the self-similar distributions for number of faces for ϕ = 0.20 and different initial distributions: the overlapping of these functions is a further evidence for the scaling regime.</p><!><p>So far, experiments and simulations have investigated self-similar growth primarily in dry foams [5, 6, 7]. Together, the growth rates in figure 5 and distribution functions in figures 6 and 7, and the power laws seen in figure 8, show that the foams reach a self-similar growth regime after a transient of ≈ 600 MCS for all wetnesses and all three initial volume distributions.</p><p>Radius scaling exponents in experimental foams depend on the wetness, dropping from 1/2 to 1/3 as ϕ goes from zero to one. In our simulations, the radius scaling exponent for ϕ = 0.20 is significantly smaller than in the experiments of Isert [10] but similar to that in Fortuna et al. [14]. This difference between experiments, as well as between experiment and simulation, may depend on the relative transport rates in thin walls and Plateau borders, which we will investigate in future work.</p>