Numerical simulations of three-dimensional foam by the immersed boundary method

https://doi.org/10.1016/j.jcp.2014.03.016Get rights and content

Abstract

In this paper, we extend (Kim et al., 2010 [13]) to the three-dimensional dry foam case, i.e., a foam in which most of the volume is attributed to its gas phase. Dry foam dynamics involves the interaction between a gas and a collection of thin liquid–film internal boundaries that partitions the gas into discrete cells or bubbles. The liquid–film boundaries are flexible, contract under the influence of surface tension, and are permeable to the gas which moves across them by diffusion at a rate proportional to the local pressure difference across the boundary. Such problems are conventionally studied by assuming that the pressure is uniform within each bubble. Here, we introduce instead an immersed boundary method that takes into account the non-equilibrium fluid mechanics of the gas. To model gas diffusion across the internal liquid–film boundaries, we allow normal slip between the boundary and the gas at a velocity proportional to the (normal) force generated by the boundary surface tension. We implement this method in the three-dimensional framework, and test it by verifying the 3D generalization of the von Neumann relation, which governs the coarsening of a three-dimensional dry foam.

Introduction

We use an immersed boundary (IB) method to simulate a “dry” foam in 3D, i.e., a gas–liquid mixture in which most of the volume is attributed to its gas phase. One interesting phenomenon, which is called diffusive coarsening, is the evolution in bubble size and topological structure that occurs as a result of gas exchange between bubbles [1]. This gas exchange occurs by diffusion through the thin liquid films that separate one bubble from another. The diffusive flux of gas through such a film is proportional to the pressure difference between the two bubbles that are separated by that film. In a dry foam, those thin films forming the faces of the roughly polyhedral bubbles are called “lamella”, and the junctions of these films are called “Plateau borders” which form the tubes in which liquid may flow. The vertices, at which typically four Plateau borders intersect, are called “nodes”.

In 1952, von Neumann [2], [3] showed that the rate of change of the area of a given bubble (a curved polygon) in a two-dimensional dry foam is independent of bubble size and solely dependent on the number of walls (or edges). The derivation of this remarkable result is based on the fact that the net rate of outward diffusion of gas per unit length through a wall of the (two-dimensional) bubble is proportional to the pressure difference across that wall, which in turn is equal to the product of the surface tension γ and the curvature κ. The rate of change of the area of a bubble isdAdt=MγΓκdl=2πMγ(n/61), where M is a permeability coefficient, the curve Γ is the closed boundary of the bubble, dl is the arc length along Γ, and n is the number of walls. The second equality can be proved by using the fact that the boundary of the bubble is a closed curve and that each of the exterior (turning) angles at the vertex (triple junction) must be equal to 2π/6. Thus, the two-dimensional von Neumann law is purely topological, since the rate of change of area of any particular bubble depends only on its number of edges (the constants M and γ being the same for all of the bubbles in a given foam).

In 3D, the net rate of outward diffusion of gas per unit area through a wall of the bubble is proportional to γH where H is the so-called mean curvature, that is the sum of the two principal curvatures of the wall surface. Let S(t) be the surface of a bubble, and let V(t) be the enclosed volume, then the rate of change of volume of the bubble isdVdt=MγS(t)HdS, where dS is the area element on S(t). The sign convention for H is that inward curvature is considered to be negative. Note that the integral in (2) is to be understood as the sum of the corresponding integrals over each of the smooth faces of the bubble; there is no contribution from the Plateau borders or from the nodes, where the curvature is infinite, but where gas transport is considered negligible.

It has been a long-standing mathematical problem to establish an explicit formula for the surface integral of the mean curvature along a given surface S(t). Based on the assumption that each of the dihedral angles along junction edge is equal to 2π/6 and on the ingeniously defined mean width of a bubble domain, MacPherson and Srolovitz [4] have derived a general formula for the 3D von Neumann relation, see [4] and its supplementary material.

In the present paper, we derive a discretized version of the 3D von Neumann relation. Let the foam boundary S(t) be modeled as a faceted surface, with triangular facets. Then the discretized version of the 3D von Neumann relation can be written asdVdt=MγeE0LeθeMγ2eE1Leθe, where Le is the length of edge e of the faceted surface, and θe is the angle between the two facets with the same edge e. Here, θe>0 if edge e is convex when viewed from outside the surface. The set Ei, where i=0, 1, or 2, is the collection of edges which have exactly i of their two vertices on a junction edge, which is defined as an edge on which three or more facets meet, see Fig. 1. Thus a junction edge is part of a Plateau border (see above). Such an edge belongs to three or more faces of the foam. It is member of the set E2 (bold lines), and, as such, it does not appear at all on the right-hand side of (3). An edge in E1 (dotted lines) share only one vertex with a Plateau border, and an edge in E0 (solid lines) has both vertices in the interior of some face. The detailed derivation of the formula (3) is given in Appendix B, and it is used for comparison with our computational results. A similar formula can also be found in [5].

In [1], [6], [7], [8], the authors simulated the evolution of a two-dimensional dry foam within the framework of the following assumptions: Laplace–Young condition, Plateau's rule, and the von Neumann relation (1). There are other foam simulations that take the fluid dynamics into account. In [9], a numerical study based on a boundary integral formulation is presented to simulate two-dimensional, doubly periodic, diluted and concentrated emulsions or foams structures in a simple shear Stokes flow. Recent work in [10] used a nonsingular boundary integral method to simulate the three-dimensional wet foam drop formation and its dynamics in simple shear flow. More numerical works on the two- and three-dimensional dry or wet foams can also be found in [11], [12] and the reference therein.

In [13], three of the present authors have introduced an immersed boundary (IB) method to simulate the fluid dynamics of a two-dimensional dry foam by modeling the gas phase of the foam as a viscous incompressible fluid, and the liquid phase as a massless network of permeable internal boundaries under surface tension. The gas diffusion through the liquid phase of the foam was modeled by allowing the internal boundaries to slip relative to the fluid, at a velocity (speed and direction) proportional to the boundary force [13], [14], [15]. The paper [13] has shown that the IB method is a proper tool to simulate the 2D foam dynamics by verifying the 2D von-Neumann relation (1) and simulating foams with an arbitrary shape.

Here we extend the methodology introduced in [13] and present the IB method to study the three-dimensional dry foam dynamics. Although we basically use the same idea as in [13], the mathematical derivation and numerical implementation are complicated enough to require a separate detailed description. In 3D, the bubble boundaries are surfaces which we model by triangulation, and the discrete force density has to be evaluated at the vertices of the triangles. The derivation of this discrete force density from the surface energy defined on the triangulated surface is given in Appendix A. Moreover, an algorithm is needed to maintain the resolution of each internal boundary within predetermined bounds despite large changes in the lengths of edges and in the areas of the triangular facets.

The rest of the paper is organized as follows. In Section 2, we describe the equations of motion of the foam in immersed boundary formulation. These are the typical IB equations of motion, generalized to handle a permeable boundary under surface tension. The numerical implementation including external boundary conditions is described in Section 3. In Section 4, we validate the method by doing a convergence study and by showing numerically that the discretized 3D von Neumann relation (3) is satisfied. Then we consider two more complicated cases: one is the case in which a dynamic non-equilibrium flow of gas interacts with the foam boundaries, and the other case involves a foam with multiple bubbles. Conclusions and future work are discussed in Section 5.

Section snippets

Equations of motion

In this section, we state the equations of motion of a three-dimensional dry foam, in which the boundaries between bubbles are idealized as massless surfaces under a constant surface tension γ. These boundaries are assumed permeable to the gas phase of the foam, with permeability coefficient M. The gas phase is modeled as a viscous incompressible fluid [13]. The constant density and viscosity of the gas are denoted by ρ and μ, respectively.

In the following formulation, the parameter pair (r,s)

Computational procedure

What has been stated so far is the mathematical formulation of the problem in immersed boundary form (i.e., with delta-function forces instead of explicit boundary conditions). For the numerical implementation, we use a first-order IB method, generalized to take a permeable foam boundary into account [13], [14]. We use a superscript to denote the time level; thus, Xn(r,s) is a shorthand for X(r,s,nΔt), where Δt is the duration of the time step, and similarly for all other variables.

Throughout

Initial setting

Consider a computational domain Ω filled with an incompressible fluid in which foam boundaries are immersed and divide the fluid into several bubbles, see Fig. 4. Throughout this paper, we fix the computational domain Ω=[1,1]×[1,1]×[1,1] mm3, and the no-slip (rigid) boundary conditions are imposed along the circular cylinder x2+y2=1 mm2, unless otherwise stated. Note that the rigid cylinder extends from the top to the bottom of the periodic domain, and that its ends are open, i.e., there is

Summary and conclusions

We have presented an immersed boundary method to simulate the fluid dynamics of a three-dimensional dry foam. We model the gas phase of the foam as a viscous incompressible fluid, and the liquid phase as a massless network of permeable internal boundaries under surface tension. The internal boundary force, generated by the surface tension, is everywhere normal to the internal boundaries. Permeability is modeled by allowing the internal boundaries to slip relative to the fluid, at a velocity

Acknowledgements

The first author was supported by National Research Foundation of Korea Grant (2012R1A1A2043238). The second author was supported in part by National Science Council of Taiwan, Republic of China under grants NSC97-2628-M-009-007-MY3, NSC98-2115-M-009-014-MY3.

References (32)

  • J. von Neumann
  • W.W. Mullins
  • R.D. MacPherson et al.

    The von Neumann relation generalized to coarsening of three-dimensional microstructures

    Nature

    (2007)
  • S. Hilgenfeldt et al.

    An accurate von Neumann's law for three-dimensional foams

    Phys. Rev. Lett.

    (2001)
  • D. Weaire et al.

    Computer simulation of a two-dimensional soap froth II. Analysis of results

    Philos. Mag. B

    (1984)
  • T. Herdtle et al.

    Numerical experiments on two-dimensional foam

    J. Fluid Mech.

    (1992)
  • Cited by (13)

    • High-resolution method for evolving complex interface networks

      2018, Computer Physics Communications
      Citation Excerpt :

      With Lagrangian methods, such as front-tracking [20], immersed-boundary [21], or arbitrary Lagrangian–Eulerian (ALE) [22] methods, the interface is represented explicitly by conforming discretization elements. Although these methods have been extended to multi-region systems [23–27,10], it is difficult to handle complex topological changes during interface-network evolution, especially in three dimensions. With Eulerian methods, such as volume-of-fluid (VOF) [28] and level-set methods [29], the interface is reconstructed from scalar fields, i.e., volume fraction or level-set field.

    • An Immersed Boundary method with divergence-free velocity interpolation and force spreading

      2017, Journal of Computational Physics
      Citation Excerpt :

      Their IB method with modified finite-difference operators (herein referred to as IBModified) was applied to a two-dimensional model of the heart, and it achieved improvement in volume conservation by one-to-two orders of magnitude compared to IBCollocated. Nevertheless, a major drawback of IBModified that limits its use in applications is its complex, non-standard finite-difference operators that uses coefficients derived from the regularized delta function (but see [21,22] for applications). To address the issue of spurious currents across immersed structure supporting extremely large pressure differences, Guy and Strychalski [39] developed a different extension of the IB method that uses non-uniform Fast Fourier Transform [8,12] (NUFFT) to generate “spectral” approximations to the delta function, which also has superior volume conservation.

    • Simulation of a Soap Film Möbius Strip Transformation

      2017, East Asian Journal on Applied Mathematics
    View all citing articles on Scopus
    View full text