LBfoam: An open-source software package for the simulation of foaming using the Lattice Boltzmann Method

TL;DR

LBfoam is an open-source simulation tool that models foam formation and behavior using the Lattice Boltzmann Method, incorporating physics like gas diffusion, bubble dynamics, and surface tension for large-scale parallel computations.

Contribution

The paper introduces LBfoam, a comprehensive open-source software package for simulating foaming processes with detailed physics and scalable parallel performance.

Findings

Successfully simulates bubble growth and interactions

Predicts final foam structures accurately

Enables large-scale parallel foam simulations

Abstract

This paper presents a 2D/3D Free Surface Lattice Boltzmann Method simulation package called LBfoam for the simulation of foaming processes. The model incorporates the essential physics of foaming phenomena: gas diffusion into nucleated bubbles, bubble dynamics and coalescence, surface tension, the stabilizing disjoining pressure between bubbles, and Newtonian and non-Newtonian rheological models. The software can simulate the growth and interaction of bubbles, and predict final foam structures. The implementation is based on the Palabos library (in C++), which enables large-scale parallel simulations. The software is freely available under the GNU Affero General Public License version 3 at: https://github.com/mehdiataei/LBfoam

Tables

Table 1. Table 1 : Advection-diffusion simulation parameters.

$v_x$	$0.2$
$q$	$1e^{-6}$
$\nu$	$0.25$
$D$	$0.06$
domain size	$200\times 20$
$k_H$	$0.01$
$p_g$	$\frac{1}{3}$
$L$	$40$

Equations

$$f_i(\vec{x}+{\vec{e}}_i{\delta }_t,t+{\delta }_t)=f_i(\vec{x},t)+\hat{C}(\vec{x},t)$$

$${\vec{e}}_i^{\top }=\{\begin{array}{ll}(0,0,0),& i=0\\ (\pm 1,0,0),(0,\pm 1,0),(0,0,\pm 1),& i=1,2,3,4,5,6\\ (\pm 1,\pm 1,0),(\pm 1,0,\pm 1),(0,\pm 1,\pm 1)& i=7,\dots ,18\end{array}$$

$$\text{Collision}:f_i^c(\vec{x},t)=f_i(\vec{x},t)+\hat{C}(\vec{x},t)$$

$$\text{Streaming}:f_i(\vec{x}+{\vec{e}}_i{\delta }_t,t+{\delta }_t)=f_i^c(\vec{x},t)$$

$$\ \hat{C}(\vec{x},t)=-\frac{\delta_{t}}{\tau}\left(f_{i}(\vec{x},t)-f_{i}^{eq}(\vec{x},t)\right)-F_{i}\$$

$$f_i^{eq}(\rho ,\vec{u})=\rho w_i\left(1+\frac{{\vec{e}}_i^{\top }\vec{u}}{c_s^2}+\frac{({\vec{e}}_i^{\top }\vec{u}{)}^2}{2c_s^4}-\frac{{\vec{u}}^2}{2c_s^2}\right)$$

$$F_i=\rho w_i\left(\frac{{\vec{e}}_i^{\top }\vec{u}}{2c_s^4}-\frac{\vec{u}-{\vec{e}}_i}{c_s^2}\right)\cdot \vec{a}$$

$$\rho (\vec{x},t)=\sum _{i=0}^N{f}_i$$

$$\vec{u}(\vec{x},t)=\frac{1}{\rho }\sum _{i=0}^N{f}_i{\vec{e}}_i$$

$$p(\vec{x},t)=\rho (\vec{x},t)c_s^2$$

$$\frac{\partial \alpha }{\partial t}+\mathbf{u}\cdot \nabla \alpha =0$$

$$\alpha =\{\begin{array}{ll}0,& \forall \vec{x}\in G\\ 1,& \forall \vec{x}\in L\\ 0<\alpha <1& \forall \vec{x}\in I\end{array}$$

$$\alpha =\frac{M}{\rho {\delta }_l^3}$$

$$\alpha (\vec{x},t+1)=\alpha (\vec{x},t)+\frac{1}{\rho (\vec{x},t)}\sum _{i=0}^N\Theta (f_{\bar{i}}(\vec{x}+{\vec{e}}_i,t)-f_i(\vec{x},t))$$

$$\Theta =\{\begin{array}{ll}0,& \vec{x}+{\vec{e}}_i\in G\\ 1,& \vec{x}+{\vec{e}}_i\in L\\ \frac{1}{2}\left(\alpha (\vec{x},t)+\alpha (\vec{x}+{\vec{e}}_i,t)\right)& \vec{x}+{\vec{e}}_i\in I\end{array}$$

$$f_i^u(\vec{x}-{\vec{e}}_i,t+1)=f_i^{eq}({\rho }_g,{\vec{u}}_I)+f_{\bar{i}}^{eq}({\rho }_g,{\vec{u}}_I)-f_{\bar{i}}(\vec{x}-{\vec{e}}_i,t+1)$$

$${\rho }_g=\frac{p_g-2\gamma \kappa (\vec{x})}{c_s^2}\vec{x}\in I$$

$$p_g^i=\frac{m_g^i\mathcal{R}T}{V^i}$$

$$\nu =c_s^2\left(\tau -\frac{1}{2}\right)$$

$$\frac{\mu -{\mu }_{\infty }}{{\mu }_o-{\mu }_{\infty }}={\left(1+(\lambda \dot{\gamma }{)}^a\right)}^{\frac{n-1}{a}}$$

$$\frac{\tau -{\tau }_{\infty }}{{\tau }_o-{\tau }_{\infty }}={\left(1+(\lambda \dot{\gamma }{)}^a\right)}^{\frac{n-1}{a}}$$

$$\frac{\partial c}{\partial t}+\nabla \cdot (c\mathbf{u})=\nabla \cdot (D\nabla c)+q$$

$$c(\vec{x},t)=\sum _i^N{g}_i(\vec{x},t)$$

$$g_i(\vec{x}+{\vec{e}}_i,t+1)=g_i(\vec{x},t)+\frac{1}{{\tau }_g}\left(g_i^{eq}(\vec{x},t)-g_i(\vec{x},t)\right)+w_{i}q$$

$$D=c_s^2\left({\tau }_g-\frac{1}{2}\right)$$

$$g_i^{eq}(\vec{x},t)=w_{i}c(\vec{x},t)\left(1+\frac{{\vec{e}}_i^{\top }\vec{u}}{c_s^2}\right)$$

$$\Delta m_g^j=\rho \sum _{\vec{x}\in I^j}\left(\sum _{\vec{x}+{\vec{e}}_i\in F}[g_{\bar{i}}(\vec{x}+{\vec{e}}_i,t)-g_i(\vec{x},t)]-c(\vec{x},t)[\alpha (\vec{x},t)-\alpha (\vec{x},t-1)]\right)$$

$$c(\vec{x},t)=k_H{p}_g\forall \vec{x}\in I$$

$$g_i^u(\vec{x},t+1)=g_i^{eq}(c,{\vec{u}}_I)+g_{\bar{i}}^{eq}(c,{\vec{u}}_I)-g_{\bar{i}}(\vec{x},t+1)$$

$$\Pi =\{\begin{array}{ll}0& d>d_{max}\\ k_{\Pi }(d_{max}-d)& d<d_{max}\end{array}$$

Full text

https://github.com/mehdiataei/LBfoam

\geometry

top=2.5cm,bottom=4.5cm, footskip=1.5cm

1 Introduction

Foams are used in a wide range of products in the automotive, aerospace, furniture, packaging and insulation industries due to desired mechanical characteristics such as weight reduction, improved strength per weight, and energy absorption. Foams are formed by nucleation and growth of gas bubbles in a liquid medium. Initially either, a metered amount of gas is dissolved in the liquid, or the gas is produced by a foaming agent that is mixed with the liquid. A rapid drop in solubility (e.g. induced by a pressure drop) gives rise to nucleated gas bubbles, which grow due to the diffusion of gas molecules from the liquid. Ultimately, bubbles are separated by thin liquid regions (lamellae), forming a cellular structure. Simulating this foaming process is essential for developing foam structures with desired properties as, for example, homogeneous foams tend to improve the mechanical strength per weight of a final structure.

Developing a numerical model for foam evolution is challenging due to the complexity of the physical phenomena. There have been many attempts in the literature to model the foaming process by introducing simplifications. Various single bubble models (i.e., the “Cell Model”) estimate the growth rate of an isolated bubble as a 1D diffusion problem, but they neglect bubble interactions and assume that the foam is a periodic array of spherical bubbles (e.g. [1, 2, 3, 4, 5, 6]). A comparison of these models can be found in [7]. In a number of 2D models, the transport equations are solved for only one section of a lamella separating adjacent bubbles, by assuming that the bubbles are organized in known arrangements such as hexagonal arrays [8, 9, 10]. The arbitrary Lagrangian-Eulerian sharp interface algorithm has also been used to simulate bubble growth in 2D [11], although this model is prone to errors for large bubble deformations and cannot capture bubble coalescence and rupture. The Shan-Chen model [12] based on the Lattice Boltzmann Method (LBM) has shown promising results for the simulation of bubble growth and interaction for small liquid-gas density ratios [13].

There is another class of foaming simulations over much larger length scales, where the interest is not individual bubble growth and interactions, but rather to model foam flow into a given mold geometry, i.e. mold-filling. These models omit details of foaming at the bubble-scale, but use bubble scale models to extract certain constants and information that are then incorporated into particle tracers, surrogate models, or population balance equations, to approximate the effects of bubble growth on macroscale fluid properties such as fluid density and viscosity (e.g. [14, 15]).

Here we present a foaming simulation software, which is a modified version of the Free Surface Lattice Boltzmann Method (FSLBM) that is originally developed to model fluid flows when a large fraction of a domain is occupied by gas [16, 17, 18, 19, 20, 21, 22, 23, 24, 25]. The pressure and gas content of each bubble is monitored throughout the simulation, and sudden topological variations due to bubble coalesce, splitting, merging, and bursting are accounted for. The numerical model couples the gas advection-diffusion with the fluid flow model, for which bubble interfaces are arbitrarily-shaped immersed boundaries at which gas concentration must be specified. A short-range force between bubble interfaces (i.e. disjoining pressure) is considered to allow formation of stable lamella between bubbles.

This simulation model encompasses all of the essential physics of foaming, and can be used to analyze 2D and 3D foaming processes such as polymer foam injection molding and metal foaming for prediction of the final foam structure, and for understanding phenomena such as topological changes (e.g. bubble coalescence), the occurrence of foam drainage, and the dissolution of bubbles.

The most closely related work is that of Körner et al. [18] on a 2D FSLBM foaming model. They simulated the expansion of gas bubbles in molten aluminum, and studied the influence of viscosity, surface tension, and mold constraints on the final structure of a foam. Since their source code is not publicly available, the scope and capabilities of the model are not known. And although a 3D version of the model has been developed [19], advection-diffusion of gas is not included, and bubble growth is modeled by increasing bubble gas content in proportion to bubble surface area.

The software presented here extends the Palabos library [26] in a non-intrusive way, so that the user can easily access useful features of Palabos including efficient MPI parallelization for large-scale simulations, access to various collision operators, the use of complex geometries for the simulation domain using STL meshes, different viscosity models such as Newtonian, Smagorinsky, and Carreau-Yasuda, and various boundary conditions.

In the remainder of the paper, we present our formulation of FSLBM for foaming, followed by several validation and simulation results of different foaming phenomena.

2 Free Surface Lattice Boltzmann Method

The Lattice Boltzmann Method (LBM) is widely used to perform CFD simulations. The Chapman–Enskog expansion [27] shows that a discretized Boltzmann equation with a collision kernel ($\hat{C}$) can solve the incompressible Navier-Stokes equations for fluid flow. The discretized lattice Boltzmann equation reads

[TABLE]

The particle distribution function (PDF) for each lattice direction at time $t$ is given by $f_{i}(\vec{x},t)$ where $\delta_{t}$. In LBM, space is discretized with square (in 2D) or cubic lattices (in 3D) with a finite set of $N$ discrete lattice velocities $\vec{e}_{i}$ ($i=0,\dots,N-1$). For example, in 2D and 3D, one can use D2Q9 and D3Q19 lattices, respectively. The velocity sets for the D3Q19 lattice are:

[TABLE]

Fig. 1 shows the lattices and their respective weights $w_{i}$. The spatial resolution of a domain in 3D corresponds to the lattice width $\delta_{l}$, so that a box of length $l_{x}$, $l_{y}$ and $l_{z}$ is comprised of $\frac{l_{x}}{\delta_{l}}$, $\frac{l_{y}}{\delta_{l}}$ and $\frac{l_{z}}{\delta_{l}}$ lattices (also referred to as cells in this paper).

Equation (1) can be decomposed into collision and streaming steps:

[TABLE]

[TABLE]

where $f^{c}_{i}$ denotes the post-collision PDF. The collision operator $\hat{C}$ is often replaced by the Bhatnagar-Gross-Krook (BGK) model for isothermal systems with an external force term $F_{i}$ [28]:

[TABLE]

where $\tau$ is the time interval between particle collisions. The equilibrium distribution function $f_{i}^{eq}$ is the truncated second-order Taylor expansion of the Maxwell distribution function with respect to flow velocity $\vec{u}$:

[TABLE]

$\rho$ is the density and $c_{s}$ is the speed of sound, that has a value of $1/\sqrt{3}$ for the D2Q9 and D3Q19 lattices. Both $\rho$ and $\vec{u}$ are functions of $\vec{x}$ and $t$. The force term $F_{i}$ in Eq. 5 accounts for external forces such as gravity $g$. For an external force generating an acceleration $\vec{a}$, we have [28]:

[TABLE]

Palabos provides access to several other collision operators including the Multiple-Relaxation-Time (MRT), Two-Relaxation-Time (TRT), Regularized Lattice Boltzmann (RLB), and the Entropic Lattice Boltzmann (ELB) models, that can be used instead of the BGK model. A comparative study of these models can be found in [29].

The macroscopic quantities $\rho$ and $\vec{u}$ are calculated by summation of PDFs:

[TABLE]

[TABLE]

and pressure is proportional to the density:

[TABLE]

The bounce-back scheme [30] is used at wall nodes, meaning that particles colliding with a wall have their momentum reversed, which translates to a no-slip boundary condition. Hereafter, we use $\delta_{l}=\delta_{t}=1$ for simplicity.

2.1 Interface Capturing

Interfaces between gas and liquid phases in a foam are captured using a volume-of-fluid (VOF) method, where a dimensionless scalar variable $\alpha$ is introduced to track the motion of interfaces throughout the computational domain by solving a pure advection transport equation:

[TABLE]

$\alpha$ is equal to one in liquid cells, zero in gas cells, and it has a value between zero and one in all interface cells (see Fig. 2):

[TABLE]

$G$, $L$, and $I$ represent gas, liquid, and interface cells, respectively.

Here, instead of solving Eq. 11 explicitly, we use a fast mass tracking algorithm [19] that is commonly used in FSLBM simulations, which takes advantage of the evolution of PDFs to track $\alpha$ in the domain. In this algorithm, $\alpha$ is defined as:

[TABLE]

where $M$ is the liquid mass in a cell.

The value of $\alpha(\vec{x},t+1)$ is found by calculating the mass exchange between neighboring cells $\vec{x}$ and $\vec{x}+\vec{e}_{i}$:

[TABLE]

where $\bar{i}=-i$. The parameter $\Theta$ in Eq. 14 weights the mass exchange between two interface cells by their average volume fraction:

[TABLE]

Equation (14) conserves mass locally since the mass that leaves one cell is distributed to neighboring cells. An interface cell becomes a liquid cell when $\alpha(\vec{x},t)\geq 1$, and an interface cell becomes an empty cell when $\alpha(\vec{x},t)\leq 0$. Where there is an excess or shortage of mass in a cell (i.e. $\alpha(\vec{x},t)>1$ or $\alpha(\vec{x},t)<0$), the surplus or shortage is uniformly distributed to neighboring interface cells. Other interface capturing and tracking methods are also available e.g. [31, 32].

2.2 Free Surface Formulation

The liquid-gas interaction is modelled as a free surface. Since the ratios of liquid to gas density and viscosity are often very large, the hydrodynamics of the gas phase are ignored for simplicity. Only the influence of the pressure of the gas phase is considered.

The free surface simplification means that the lattice Boltzmann equation is not solved for in the gas cells. As a result, in the streaming step (Eq. 4), the PDFs arriving from gas cells are undefined in interface cells. These undefined PDFs ($f_{i}^{u}$) are imposed such that the pressure boundary condition at the free surface is enforced, and so that the velocity of gas cells is equal to the adjacent liquid velocity $\vec{u}_{I}$. For this purpose, the unknown distribution functions at the interface are set as [19]:

[TABLE]

The effect of surface tension is considered in the model by modifying the gas density $\rho_{g}$ in Eq. 16 [19]:

[TABLE]

$\gamma$ is the surface tension, and $\kappa(\vec{x})$ is the local curvature. Palabos provides different methods for curvature calculation, including using the gradient of a smoothed $\alpha$ scalar field [33], and a height function [34]. In Eq. 17, gas pressure $p_{g}$ for each individual bubble $i$ is given by the ideal gas law:

[TABLE]

where $m^{i}_{g}$, $\mathscr{R}$, $T$, and $V^{i}$ are the mass of a bubble, gas constant, temperature, and volume of a bubble.

2.3 Viscosity

In the collision operator (Eq. 5), $\tau$ is related to kinematic viscosity $\nu$ by [28]

[TABLE]

For a Newtonian fluid, $\tau$ is a constant. For shear thinning polymers, the Carreau-Yasuda model [35] can be used. In this model, the apparent viscosity $\mu$ is given by

[TABLE]

where $n$, $a$, and $\lambda$ are empirically-determined material coefficients, and $\mu_{o}$ and $\mu_{\infty}$ are zero-shear viscosity and the viscosity at infinite shear-rate. In terms of $\tau$, Eq. 20 can be written as:

[TABLE]

3 Foaming Model

3.1 Advection-Diffusion

The advection and diffusion of the dissolved gas $c$ within the liquid and into bubbles is governed by an advection-diffusion equation:

[TABLE]

where $q$ is a source term that takes into account gas generation due to, for example, chemical reactions. This equation is solved using another distribution function $g_{i}(\vec{x},t)$, such that the summation of $g_{i}$ gives the gas concentration (in terms of mass fraction) at location $\vec{x}$ at time $t$:

[TABLE]

The lattice Boltzmann equation for the evolution of $g_{i}(\vec{x},t)$ is given by:

[TABLE]

We solve the advection-diffusion equation using lattices with fewer velocity vectors [36]: on D2Q5 (in 2D) or D3Q7 (in 3D) lattice topologies, that coincide with the fluid flow lattices. This calculation requires less computational time because of the reduced number of velocities. $\tau_{g}$ is the relaxation time for the advection-diffusion equation, which relates to the diffusion constant $D$ of the gas in the liquid according to the following equation:

[TABLE]

where $c_{s}^{2}$ is equal to $1/3$ and $1/4$ for the D2Q5 and D3Q7 lattices, respectively. The equilibrium distribution in the advection-diffusion equation is given by:

[TABLE]

where $\vec{u}$ is calculated from Eq. 9.

The amount of gas diffusing into bubble $i$, $\Delta m_{g}^{j}$, is calculated based on the concentration gradient at the liquid-bubble interface:

[TABLE]

where $I^{j}$ is the liquid-gas interface of bubble $j$. The last term in Eq. 27 deducts the pure advection portion of gas transport, which is not diffused into the bubble.

3.2 Henry’s Law Boundary Condition

In many foam media such as polymers, the concentration of gas at a bubble interface obeys Henry’s law, which expresses a linear relationship between pressure and gas concentration:

[TABLE]

Similar to Eq. 16, this boundary condition is enforced by setting the unknown distribution coming from the gas cells $g^{u}_{i}$ to

[TABLE]

where $c$ is given by Eq. 28 [19].

3.3 Disjoining Pressure

To stabilize the lamella between bubbles, a disjoining pressure $\Pi$ is assumed to be active within a maximum distance $d_{max}$ between two different bubble interfaces, with a linear dependence on the distance between the interface cells $d$ [19]:

[TABLE]

where $k_{\Pi}$ is a constant. In this work, the disjoining pressure is active up to four lattice cells i.e., $d_{max}=4\delta_{l}$.

This disjoining pressure is a repulsive force that resists bubble coalescence and can stabilize a lamellae. Physically, $\Pi$ originates from the variation of Gibbs free energy with the distance $d$ between two interfaces that correspond to different bubbles [37, 38]. Similar to the surface tension implementation, the disjoining pressure is added to $\rho_{g}$:

[TABLE]

Calculating $\Pi$ is a two-step process as shown in Fig. 3. First, for each interface cell, a traversal ray tracing algorithm [39] is used to move along the bubble $i$ interface normal $\vec{n}$ as far as $d_{max}$, in search of an interface belonging to another bubble $j$. If an interface cell is detected, $d^{{}^{\prime}}$ is calculated as the width of the traversed cells minus $1-\alpha$, and projecting that distance along the normal vector to find $d+d^{*}$. The distance $d$ is then determined by subtracting the distance from the center of the interface cell to the interface $d^{*}$, which is calculated by geometrically reconstructing the interface of the bubble using a piecewise linear interface construction (PLIC) [40].

A full description of the implementation of the ray tracing algorithm is given in A. Since we use cubic or square lattices, we use the explicit analytical expressions developed by Scardovelli and Zaleski [41] for PLIC calculations based on the implementations provided in Ref. [42], which are faster than iterative methods.

Fig. 4 shows the effect of $\Pi$ on lamella stabilization between two growing bubbles. In all the three cases, the initial parameters are identical, except for the value of $k_{\Pi}$ that regulates the strength of $\Pi$. The growing of the bubbles draw the bubbles together until at $t=2800$, in the absence of a disjoining pressure ($k_{\Pi}=0$), the bubbles coalesce upon contact (the exact timestamp is not shown in Fig. 4), and no stable lamella forms. At $k_{\Pi}=0.005$ other forces overcome the disjoining pressure and the bubbles still coalesce, but the disjoining pressure delays the coalescence until $t=3400$. At $k_{\Pi}=0.08$, the disjoining pressure is strong enough to prevent coalescence, and so a stable lamella forms between the bubbles.

The formation of lamella is more apparent in a system with multiple bubbles. Fig. 5 shows a number of bubbles ascending in a bubble column due to buoyancy. Although 18 bubbles coalesce, the disjoining pressure eventually stabilizes the lamella between the remaining bubbles and enables the formation of a cellular structure near the liquid surface.

3.4 Nucleation

The location and radius of nucleated bubbles must be specified at the beginning of each simulation. A nucleation probability field (derived on the basis of a nucleation model) may be used to specify the distribution and initial radius of the nucleated bubbles, as shown in Fig. 6. If we pick random points in the domain, the nucleated bubbles may overlap, as the bubbles are disks (2D) or spheres (3D), not points. For this purpose, a Poisson Disk Sampling algorithm is implemented based on the work of [43], which enforces that the distance between selected nucleation sites be at least equal to the radius of the bubbles plus one $\delta_{l}$. Bubbles can also be distributed in an “organized” fashion (see Fig. 6).

The initial configuration of the bubbles have a major effect on the final structure of the foam. Section 4.5 presents an example of the effect of different initial bubble distributions.

3.5 Software Structure

Fig. 7 shows the directory structure of the LBfoam library. Many of the LBfoam functions and classes either call or extend the Palabos library classes. The Palabos open source framework provides basic tools for parallelization of LBM simulations using Message Passing Interface (MPI) and Lattice Boltmzann solvers. LBfoam is written in a non-intrusive manner with respect to Palabos, except for a few exceptions, to allow researchers to use it along with future versions of Palabos. Currently LBfoam uses Palabos v2.0, which is the latest version.

The algorithms folder contains the PLIC, ray-tracing, Poisson disk-sampling algorithm, and other related functions. The $dynamics$ folder contains the classes for advection-diffusion coupling and applying the Henry’s law boundary condition. The $bubble$ folder includes the classes for calculating gas diffusion into each bubble, updating bubble pressures, tracking bubble coalescence, calculation of disjoining pressure, and updating bubbles gas content. The $models$ folder includes a 2D FSLBM model which is not available in Palabos. The directory $examples/lbfoam$ include a variety of demo cases. With the ongoing development of the library, this directory structure may change in the future.

$lbfoam2D.h$ and $lbfoam3D.h$ header files provide access to LBfoam classes for 2D and 3D simulations, respectively. Each LBfoam function/class must be called using the namespace called lbfoam (e.g. lbfoam::PLIC2D).

Palabos uses a data-structure called “multi-block” that is composed of a number of Cartesian meshes, that combined create a computational domain. These cartesian meshes are distributed to a number of processing units. At each iteration, every processor tags the cells on its sub-domain that belong to a bubble with a unique ID using a flood-fill algorithm, and attaches bubble information such as volume and gas content to the ID. The ID of each bubble is unique and global among processors. In the case of foaming, there are two situations that require special communication (data transfer) between processors. The first relates to the calculation of the disjoining pressure: when a bubble interface is closer than $d_{max}$ from a sub-domain boundary, an adjacent bubble could be in another sub-domain. In this case, the communication envelope between the two processors must be expanded as much as $d_{max}$. Second, when a bubble occupies multiple sub-domains, bubble gas content must be shared among all sub-domain processors. Since gas can diffuse into a bubble from different sub-domains, this communication is necessary to calculate the correct value of gas content at each iteration. Bubble merging or splitting is carefully monitored at each iteration as the gas content must be distributed to the resulting bubble(s) in proportion to their volume(s).

While bubble tracking involves multiple communications between distributed CPUs, we evaluated the performance of our model on up to 100 nodes (with 40 processors each) on the SciNet Niagara supercomputer [44] at the University of Toronto, and were able to achieve close to linear parallel scaling.

4 Validation and Sample results

4.1 Advection-Diffusion

We begin by examining the solution of the advection-diffusion equation to an available analytical solution. The schematic of the problem is shown in Fig. 8. A 1D section of liquid (in red) of length $L$ moves at a constant velocity ($v_{x}$) within a rectangular domain, and is surrounded by gas on both sides. Boundaries are periodic, and gas is generated at a constant rate within the liquid ($q$ is the source term), and diffuses into the gas phase across the interfaces. The steady-state analytical solution for the concentration $c$ is given by

[TABLE]

where $x^{\prime}=2x/L$.

The simulation parameters are given in Table 1. Fig. 9 shows that the simulation result matches the analytical solution, demonstrating that the advection-diffusion equation is being correctly solved, and that Henry’s law on the moving interface boundaries is being properly applied. We confirmed this result using other parameter sets, and successfully conducted the grid independence test using domain sizes with $0.5$ and $2x$ resolution of the domain size presented here.

4.2 3D Bubble Growth

The following is the analytical solution to the 1D “Cell Model” for a single bubble growing in an infinite reservoir (i.e. the gas concentration at infinity $c_{\infty}$ is constant [45]):

[TABLE]

$R(t)$ is bubble radius at time $t$, $R_{0}$ is the initial bubble radius, $V_{m}=RT/p$ is the molar volume and $\Delta c=c_{\infty}-c_{R}$ where $c_{R}=k_{H}p_{g}$.

To solve the same problem in LBfoam, a bubble nucleus is placed in a reservoir 30 times larger than the initial radius of the bubble, so that for the duration of the simulation $c_{\infty}$ remains almost constant (see Fig. 10). In the simulation, $R_{0}=3\delta_{l}$, $k_{H}=0.001$, $D=0.03$, $V_{m}=3$, $p_{0}=1/3$, $\rho=1$, $\gamma=g=0$, and $\nu=0.25$, and we varied $c_{\infty}$ so that $\Delta c=0.1,0.3$ and $0.5$. The radius of bubble was calculated as $R=\sqrt[3]{3V_{b}/4\pi}$ where $V_{b}$ is the bubble volume (although the bubble remained almost spherical throughout the simulations). The comparison between Eq. 33 and LBfoam is given in Fig. 11 and shows very good agreement. The small discrepancy between the results can be attributed to the assumption that $c_{\infty}$ remains constant throughout the simulation, and the related effect of the domain boundaries on the growth of the bubble.

4.3 Bubble Dissolution

In foam injection molding (FIM), bubbles nucleate due to a pressure drop at the inlet gate during filling. Due to large shear stresses during the mold-filling process, these gate-nucleated bubbles become undesirably elongated [46]. To promote foam uniformity, in high-pressure FIM, these bubbles are re-dissolved into the polymer by increasing the cavity pressure to promote the solubility of the gas, where the gas diffuses from the bubbles back into the polymer melt. Uniform bubbles are then nucleated again by a second pressure drop (due to melt shrinkage or mold opening) [47].

LBfoam can be used to study bubble dissolution. In Fig. 12, a bubble is placed in the middle of a liquid with no initial gas content, and with the source term set to zero. The non-zero concentration of gas at the bubble-liquid interface results in the diffusion of gas from the bubble to the liquid, until the bubble is fully dissolved.

4.4 Foaming

In the following, we parameterize simulations based on lattice units. The conversion between lattice units and SI units is straightforward and can be done by defining a length scale and a reference density for the simulation. Examples of unit conversions for foaming simulations are given in [18, 19]. The system parameters for each simulation are given in the figure captions.

2D foams have been investigated in fields such as rheology and microfluidics and many others (e.g. [48]). LBfoam is capable of simulating 2D foams as shown in Fig. 13, in a rectangular container initially containing 300 nuclei. The simulation domain consists of $450\times 650$ cells. Initially a portion of the domain is filled with the liquid medium, gravity $g=0$ and, the nuclei are randomly distributed, each with a radius of $3$ cells. The initial gas concentration $c_{0}$ is zero; instead the gas is generated by a source term. Note that some of the bubbles coalesce with other bubbles during the simulation.

Under the influence of gravity the behavior is different, as liquid gradually moves to the bottom of the container, leaving the top of the foam dry of liquid. This effect is known as foam drainage, and has been extensively studied for Newtonian and non-Newtonian fluids (e.g. [49, 50, 51, 52, 53, 54]), previous foam drainage studies have been largely experimental; little mathematical and numerical modeling has been done. As shown in Fig. 14, foam drainage can be observed in the simulation results, in the presence of gravity. The liquid accumulates at the bottom of the container. It can also be seen that after drainage large sections of the foam evolve into honeycomb structures forming dry and wet foam sections (see Fig. 14) similar to experimental observations [55]. The natural emergence of foam drainage illustrates the capability of the model for foam drainage studies of complex foaming phenomena.

An example of a 3D simulation is given in Fig. 15. 40 bubble nuclei are randomly distributed in a liquid medium, each with an initial radius of 3 cells. The domain is made up of $300\times 200\times 200$ cells, or 12 million in total.

In Fig. 16 we illustrate a comparison between a cross-section of the 3D simulation result and an aluminium foam, that shows good qualitative agreement.

4.5 The effect of Nuclei Distribution

The initial nuclei distribution can affect the foam structure. Fig. 17 demonstrates two cases, one in which the bubbles are initially distributed evenly, and one in which the bubbles are nucleated near the center of the liquid. This can happen in conditions where, for example, the liquid center is hotter and promotes nucleation. As shown in Fig. 17, the different initial distributions result in different structure and expansion of the foam.

5 Conclusion

We have developed a simulation software based on the lattice Boltzmann Method that can simulate foaming in 2D and 3D. An advection-diffusion equation is coupled to a fluid flow solver to account for dissolved gas diffusion into bubbles, subject to a Henry’s law boundary condition at bubble-liquid interfaces. A volume-of-fluid method is used to track and locate bubble interfaces. The model accounts for disjoining pressure between bubbles, surface tension, film drainage, and bubble dynamics, including bubble growth, deformation, coalescence, bursting and splitting. To incorporate the disjoining pressure into the model, a fast traversal algorithm is used to calculate the distance between adjacent bubbles, and the bubble interfaces are reconstructed via a PLIC algorithm. The library is based on the Palabos library, which enables large scale parallel simulations. Simulation results demonstrate the potential of the model to provide a better understanding of the evolution, structure, and rheology of foams.

The software is publicly available under the AGPL v3 license at the following Github repository, and is being continuously improved:

https://github.com/mehdiataei/LBfoam

Acknowledgements

We thank the Natural Sciences and Engineering Research Council of Canada (NSERC) and Autodesk Inc. for their financial support. Computations were performed on the Niagara supercomputer at the SciNet HPC Consortium. SciNet is funded by: the Canada Foundation for Innovation; the Government of Ontario; Ontario Research Fund - Research Excellence; and the University of Toronto.

Appendix A Ray Tracing Traversal Algorithm

The ray tracing algorithm for detection of adjacent bubbles is as follows:

For each interface cell, the algorithm marches along the interface normal $\vec{n}$ away from the bubble, until it finds an interface that belongs to another bubble. Referencing the traversal algorithm in 2D described in Algorithm 1, for an interface cell $i$, the line along the negative direction of the interface normal is defined as:

[TABLE]

where $n_{x}$ and $n_{y}$ are the components of the normal vector. First, the algorithm calculates the values of $\beta$ required for $\underline{n}$ to cross the vertical ($betaMaxX)$) and horizontal ($betaMaxY$) cell boundaries of cell $i$ ($betaMaxX$). The first neighboring cell along $\underline{n}$ is determined by the minimum of $betaMaxX$ and $betaMaxY$. In order to locate the subsequent neighboring cells along $\underline{n}$, the minimum of $betaMaxX$ and $betaMaxY$ is found after incrementing by the value of $\beta$ required to move along $\underline{n}$ for one cell length. Extending the algorithm to 3D is trivial and only requires calculating $betaMaxZ$ in the $z$ direction to calculate the minimum of $betaMaxX$, $betaMaxY$ and $betaMaxZ$.