# LBM for multiphase flows

One of the most important aspects of modeling multiphase flow is capturing interfacial dynamic effects such as surface tension. The first group of models is based on the color model of Gunstensen et al. (1991). In this model, red and blue particle distribution functions (fR and fB, respectively) were introduced to represent two different fluids. The total distribution function is the sum of the partial distribution functions:

${{f}_{\alpha }}=f_{\alpha }^{R}+f_{\alpha }^{B}\qquad \qquad( 1 )$

The LBE can be written for each phase, k.

$f_{\alpha }^{k}\left( \mathbf{x}+{{\mathbf{c}}_{i}}\Delta t,t+\Delta t \right)-f_{\alpha }^{k}\left( \mathbf{x},t \right)=\Omega _{\alpha }^{k}\left( \mathbf{x},t \right)\Delta t\qquad \qquad( 2 )$

The collision operator can be split into two components: the collisions that lead to the local equilibrium similar to the LBGK model and the collisions that contribute to the dynamics of the interfaces between the two different particles.

$\Omega _{\alpha }^{k}=\Omega _{\alpha }^{k,1}+\Omega _{\alpha }^{k,2}\qquad \qquad( 3 )$

The first collision operator is the same as the BGK operator, but the relaxation parameter only pertains to phase k.

$\Omega _{\alpha }^{k,1}\Delta t=-\frac{1}{{{\tau }_{k}}}\left( f_{\alpha }^{k}-\overline{f_{\alpha }^{k}} \right)\qquad \qquad( 4 )$

The density of each phase is:

${{\rho }_{k}}=\sum\limits_{\alpha }^{{}}{f_{\alpha }^{k}}\qquad \qquad( 5 )$

The total mixture density, ρ, is the summation of the phase densities, $\rho =\sum\limits_{k}^{{}}{{{\rho }_{k}}}$. Similarly, the phase velocity is:

${{\rho }_{k}}{{\mathbf{V}}_{k}}=\sum\limits_{\alpha }^{{}}{f_{\alpha }^{k}{{\mathbf{c}}_{\alpha }}}\qquad \qquad( 6 )$

The total mixture velocity, V, is the summation of the phase velocities, $\rho \mathbf{V}=\sum\limits_{k}^{{}}{{{\rho }_{k}}{{\mathbf{V}}_{k}}}$. The second collision operator, which contributes to the dynamics in the interfaces and generates surface tension, is

$\Omega _{\alpha }^{k,2}\Delta t=\frac{{{A}_{k}}}{2}\left| \mathbf{F} \right|\left[ \frac{{{\left( {{\mathbf{e}}_{\alpha }}\cdot \mathbf{F} \right)}^{2}}}{{{\left| \mathbf{F} \right|}^{2}}}-\frac{1}{2} \right]\qquad \qquad( 7 )$

where F is the local color gradient, defined as:

$\mathbf{F}\left( \mathbf{x} \right)=\sum\limits_{\alpha }^{{}}{{{\mathbf{e}}_{\alpha }}\left[ {{\rho }_{r}}\left( \mathbf{x}+{{\mathbf{e}}_{\alpha }} \right)-{{\rho }_{b}}\left( \mathbf{x}+{{\mathbf{e}}_{\alpha }} \right) \right]}\qquad \qquad( 8 )$

Note that F will vanish in a single phase fluid where the density is constant. The direction and length of each lattice link is represented by the vector ${{\mathbf{e}}_{\alpha }}$. The parameter Ak determines the surface tension. Defining a second collision operator alone does not keep different phases separated. To keep the phases separated, the local color momentum j, must be aligned with the direction of the local color gradient after collision.

$\mathbf{j}=\sum\limits_{\alpha }^{{}}{\left( f_{\alpha }^{r}-f_{\alpha }^{b} \right){{\mathbf{e}}_{\alpha }}}\qquad \qquad( 9 )$

$\frac{\mathbf{j}}{\left| \mathbf{j} \right|}=\frac{\mathbf{F}}{\left| \mathbf{F} \right|}\qquad \qquad( 10 )$

One of the major drawbacks of the color model is that it does not have any thermodynamic background; therefore it is limited to hydrodynamic studies.

Another two-phase flow model was originated by Shan and Chen (1993), in which the surface interaction can be maintained automatically. The second collision term is different from the color model. It is:

$\Omega _{\alpha }^{k,2}\Delta t={{\mathbf{e}}_{\alpha }}\cdot {{\mathbf{F}}^{k}}\qquad \qquad( 11 )$

Here, the effective force on the kth phase comes from a pairwise interaction between the different phases.

${{\mathbf{F}}^{k}}\left( \mathbf{x} \right)=-\sum\limits_{{{k}'}}^{{}}{\sum\limits_{\alpha }^{{}}{{{V}_{k{k}'}}\left( \mathbf{x},\mathbf{x}+{{\mathbf{e}}_{\alpha }} \right){{\mathbf{e}}_{\alpha }}}}\qquad \qquad( 12 )$

The interaction potential between the phases is Vkk', and it is defined as:

${{V}_{k{k}'}}={{G}_{k{k}'}}\left( \mathbf{x},\mathbf{{x}'} \right){{\psi }^{k}}\left( \mathbf{x} \right){{\psi }^{{{k}'}}}\left( {\mathbf{{x}'}} \right)\qquad \qquad( 13 )$

The strength of the interaction is Gkk' and ${{\psi }^{k}}\left( \mathbf{x} \right)$ is a density function for phase k at x. A discussion of the form of ψ is included in Chen and Doolen (1998) and gives a nonideal equation of state, which is needed for a two-phase flow model. The strength of the interaction only accounts for the nearest neighbor on a lattice link, $\mathbf{x}-\mathbf{{x}'}={{\mathbf{e}}_{i}}$, and is zero for any other combination. The surface tension, ${{\sigma }_{\ell v}}$, is equal to:

${{\sigma }_{\ell v}}\approx \frac{M{{G}_{k{k}'}}}{2D\left( D+1 \right)}\qquad \qquad( 14 )$

The constant D is the constant in front of the V2 term in the equilibrium distribution function. The order-parameter, M, is

$M = \frac{1}{P} {{[{{ \left( p - p \right)}^2} ]}^{1/2}} \qquad \qquad( 15 )$

The rounded overbar represents the spatial average of a value over the entire lattice. It should be noted that the collision operator, $\Omega _{\alpha }^{k,2}$, for this model does not satisfy local momentum conservation.

Another approach to solving multiphase flow problems using the LBM method is a technique that uses time splitting (Premnath et al., 2005). In this method, the interface is directly tracked using a kinematic equation for the level set function, φ. The level set function varies between -1 and 1, and the interface corresponds to the location where the level set function equals 0.

$\frac{\partial \phi }{\partial t}+\mathbf{V}\cdot \nabla \phi =0\qquad \qquad( 16 )$

The first collision term is the BGK collision term. The second collision term is:

${{\Omega }_{2}}=\left( 1-\frac{1}{2\tau } \right){{S}_{i}}\left( \mathbf{x},t \right)\qquad \qquad( 17 )$

The effect of force interactions is include in the source term, Sα. This term includes the mean-field force, F, as well as an external body force, Fext, representing gravity or a Lorentz force.

${{S}_{\alpha }}=\frac{\left( {{\mathbf{c}}_{\alpha }}-\mathbf{u} \right)\cdot \left( \mathbf{F}+{{\mathbf{F}}_{ext}} \right)}{\rho c_{s}^{2}}\overline{{{f}_{\alpha }}}\qquad \qquad( 18 )$

The mean-field force can be broken into two components, the force from the non-ideal part of the equation of state, ψ, and the surface tension force, ${{\mathbf{F}}_{s}}$.

$\mathbf{F}=-\nabla \psi +{{\mathbf{F}}_{s}}\qquad \qquad( 19 )$

The nonideal part of the equation of state was modeled by the Carnahan-Starling-van der Waals equation of state. The surface tension force is related to the density and its gradient by:

${{\mathbf{F}}_{s}}=\kappa \rho \nabla {{\nabla }^{2}}\rho \qquad \qquad( 20 )$

where the surface tension parameter,κ, is related to the surface tension, ${{\sigma }_{\ell g}}$, by:

${{\sigma }_{\ell g}}=\kappa \int{{{\left( \frac{\partial \rho }{\partial \mathbf{n}} \right)}^{2}}d\mathbf{n}}\qquad \qquad( 21 )$

where n is the normal direction of the interface. Finally, the macroscopic mass and momentum balance can be incorporated into the distribution function in nodes near the interface by writing the distribution function in terms of the gradient of the momentum flux. This is done by expanding the distribution function around its local equilibrium in terms of the Knudson number and applying the Chapman-Enskog analysis to the Taylor series expansion of the LBM equation. In this method, the interface is directly tracked; therefore, there are nodes that lie near the interface. If a node lies near an interface, the distribution function for the phase that exists on the node is known. However, the distribution on the nodes in the other phase is not known, but is calculated through the interfacial mass and momentum balances; these balances are related to the distribution function from the results of the Chapman-Enskog analysis.

A final approach to multiphase flow modeling using the LBM is the free energy approach, in which the equilibrium distribution function is defined consistently with thermodynamics. This model is introduced by Swift et al. (1995, 1996). The models mentioned here do not take into account thermal transport. Other models that capture heat transfer of two-phase systems are reviewed by Házi et al. (2002) and include modeling the energy transfer with an added distribution function.

Recently work has been done to incorporate solid/liquid/vapor interaction with the LBM method. Latva-Kokko and Rothman (2005) used the LBM to simulate the capillary rise between two horizontal plates as well as in a rectangular tube. Their model used the color-gradient LBM to distinguish the phases and showed that it gave appropriate static contact angles for both imbibition and drainage.

The development of the LBM method has the potential to be very useful in computational flow modeling of multiphase systems, however, some issues first need to be resolved. The first issue is that the models that conserve energy are not stable across a wide variety of problems. Since many multiphase systems utilize the latent heat associated with phase change, the LBM method must have the capability of conserving energy to make it a useful model for practical applications.

## References

Chen, S. and Doolen, G.D., 1998, “Lattice Boltzmann Method for Fluid Flows,” Annual Review of Fluid Mechanics, Vol. 30, pp. 329-364.

Gunstensen, A.K., Rothman, D.H., Zaleski, S., and Zanetti, G., 1991, “Lattice Boltzmann Model of Immiscible Fluids,” Physical Review A., Vol. 43, pp. 4320-4327.

Házi, G., Imre, A.R., Mayer, G., and Farkas, I., 2002, “Lattic Boltzmann methods for two-phase flow modeling,” Annals of Nuclear Energy, Vol., 29, pp. 1421-1453.

Latva-Kokko, M. and Rothman, D.H. 2005, “Static Contact Angle in Lattice Boltzmann Models of Immiscible Fluids,” Physical Review E, Vol. 72, pp. 046701(1-7).

Premnath, K.N., Nave, J-C., and Banerjee, S., 2005, “Computation of Multiphase Flows with Lattice Boltzmann Methods,” Proc. 2005 ASME IMECE, Nov. 5-11, Orlando, FL, USA.

Shan, X. and Chen, H., 1993, “Lattice Boltzmann Model for Simulating Flows with Multiple Phases and Components,” Physical Review E., Vol. 47, pp. 1815-1819.

Swift, M.R., Osborn, W.R. and Yeoman, J.M, 1995, “Lattice Boltzmann simulation of nonideal fluids,” Physical Review Letters, Vol. 75, pp. 830-833.

Swift, M.R., Orlandini, S.E., Osborn, W.R. and Yeoman, J.M, 1996, “Lattice Boltzmann simulations of liquid-gas and binary-fluid systems,” Physical Review Letters E, Vol. 54, pp.5041-5052.