# Lattice Boltzmann Model (LBM)

Most numerical codes are based on the Navier-Stokes equations, which treats a fluid as a continuous field. It is well known that a fluid is made of a discrete number of particles or molecules. Since the number of molecules is extremely large (Avogadro’s number = 6.022×1023 atoms/mole) for almost all practically sized systems, it may never be computationally viable to track each particle and its interactions with other particles. The number of molecules in a given region and the molecular interaction are described through the fluid’s density and transport coefficients (i.e., viscosity) in the continuous model.

Modeling the individual molecules for a small system over a small period of time has been achieved by molecular dynamic simulations (MDS). The computational requirements needed in these simulations can be greatly reduced if the degrees of freedom of the system are reduced. Also, instead of considering individual molecules, groups of molecules can be considered. The degrees of freedom can be reduced by restricting the movement of the molecules to a lattice. A lattice is simply a predefined direction in which a molecule can move. The first of such models is the lattice gas automata (LGA; Hardy et al., 1973). The evolution equation of the LGA is:

${{n}_{\alpha }}\left( \mathbf{x}+\Delta {{\mathbf{x}}_{\alpha }},t+1 \right)-{{n}_{\alpha }}\left( \mathbf{x},t \right)={{\Omega }_{\alpha }}\left( n\left( \mathbf{x},t \right) \right)\begin{matrix} {} & {} \\ \end{matrix}\left( \alpha =0,1\cdots ,M \right)\qquad \qquad( 1 )$

Each node has a Boolean variable, ${{n}_{\alpha }}\left( \mathbf{x},t \right)$, for each connecting lattice; if there is a molecule present it is unity, and if none is present it is zero. The number of connecting lattices is M, and the collision operator, Ωα, is dictated by collision rules. On each node, only one particle can be found moving in a given direction at a given time. This condition is imposed for memory efficiency and leads to a Fermi-Dirac local equilibrium distribution (Frisch et al. 1987). The two subsequent steps of the LGA simulations are streaming and collision. In the first step, a particle can move in one lattice site corresponding to the link direction. In the second step, collision rules are applied on sites where two or more particles meet. The rules are: the particles change configuration when two or more particles meet. Because of the Boolean variables, the LGA method is quite noisy. The Lattice Boltzmann Method (LBM) is introduced by modeling the LGA equation by single particle distribution functions, ${{f}_{i}}\left( \mathbf{x},t \right)$, instead of a Boolean variable (McNamara and Zanetti, 1988). The Lattice Boltzmann Equation (LBE) is:

${{f}_{\alpha }}\left( \mathbf{r}+{{\mathbf{c}}_{\alpha }}\Delta t,t+\Delta t \right)-{{f}_{\alpha }}\left( \mathbf{r},t \right)={{\Omega }_{\alpha }}\left( f\left( \mathbf{x},t \right) \right)\Delta t\begin{matrix} {} & {} \\ \end{matrix}\left( \alpha =0,1,\cdots ,M \right)\qquad \qquad( 2 )$

A triangular lattice structure should be used because it has sufficient symmetry to make a solution independent of how the lattice is laid out on the geometry (Frisch et al., 1986). The velocity vectors surrounding a node as well as the lattice structure are presented in Figure 2 from Boltzmann Equation. The velocity vector, ${{\mathbf{c}}_{0}}$, symbolizes a particle at rest. The discretization of the LGA and the LBE is the same when $\Delta {{\mathbf{x}}_{\alpha }}/\Delta t={{\mathbf{c}}_{\alpha }}$, therefore the lattice lengths $\Delta {{\mathbf{x}}_{\alpha }}$ are equal to ${{\mathbf{c}}_{\alpha }}\Delta t$. The LBE is a discrete velocity model of the Boltzmann equation. In the LBE, the macroscopic flow variables, density and velocity, are defined as:

$\rho =\sum\limits_{\alpha }^{{}}{{{f}_{\alpha }}}\qquad \qquad( 3 )$

$\rho \mathbf{V}=\sum\limits_{\alpha }^{{}}{{{f}_{\alpha }}{{\mathbf{c}}_{\alpha }}}\qquad \qquad( 4 )$

The collision operator can also be discretized. The collision operator for the αth component of the distribution, Ωα, contains a summation over all the velocity components, as well as all the collision cross sections, ${{\vartheta }_{\beta \alpha }}$, that the Βth component of the distribution function makes with the αth component. For a more in-depth description of the collision operator, refer to the discussion on the Boltzmann equation. The scattering angle is predetermined by the discretization of f.

${{\Omega }_{\alpha }}\left( f \right)=\sum\limits_{\beta }^{{}}{\left( {{{{f}'}}_{\alpha }}{{{{f}'}}_{\beta }}-{{f}_{\alpha }}{{f}_{\beta }} \right){{\vartheta }_{\beta \alpha }}\left| {{\mathbf{c}}_{\alpha }}-{{\mathbf{c}}_{\beta }} \right|}\qquad \qquad( 5 )$

The collision operator can be linearized, and the dimension of the matrix can be greatly reduced by expanding the distribution function into an equilibrium, $\overline{{{f}_{\alpha }}}$, and non-equilibrium, $f_{\alpha }^{*}$, component.

${{f}_{\alpha }}=\overline{{{f}_{\alpha }}}+f_{\alpha }^{*}\qquad \qquad( 6 )$

The equilibrium distribution function depends on the local macroscopic variables and should satisfy the following constraints:

$\sum\limits_{\alpha }^{{}}{\overline{{{f}_{\alpha }}}}=\rho \qquad \qquad( 7 )$

$\sum\limits_{\alpha }^{{}}{\overline{{{f}_{\alpha }}}{{\mathbf{c}}_{\alpha }}}=\rho \mathbf{V}\qquad \qquad( 8 )$

If $f_{\alpha }^{*}\ll \overline{{{f}_{\alpha }}}$, then the collision operator can be linearized.

${{\Omega }_{\alpha }}\left( f \right)={{\Omega }_{\alpha }}\left( \overline{{{f}_{\alpha }}} \right)+\sum\limits_{\beta }^{{}}{{{\Omega }_{\alpha \beta }}\left( f_{\alpha }^{*} \right)}+O\left( f{{_{\alpha }^{*}}^{2}} \right)\qquad \qquad( 9 )$

The linearized collision matrix is defined as: ${{\Omega }_{\alpha \beta }}\equiv \frac{\partial {{\Omega }_{\alpha }}\left( \overline{{{f}_{\alpha }}} \right)}{\partial {{f}_{\beta }}}$. Also, the collision operator of an equilibrium distribution is zero, ${{\Omega }_{\alpha }}\left( \overline{{{f}_{\alpha }}} \right)=0$, the LBE is:

${{f}_{\alpha }}\left( \mathbf{r}+{{\mathbf{c}}_{\alpha }}\Delta t,t+\Delta t \right)-{{f}_{\alpha }}\left( \mathbf{r},t \right)=\sum\limits_{\beta }{{{\Omega }_{\alpha \beta }}\left( {{f}_{\beta }}-\overline{{{f}_{\beta }}} \right)}\Delta t\text{ }\left( \alpha =0,1,\cdots ,M \right)\qquad \qquad( 10 )$

Since the linearized collision matrix must conserve mass and momentum, it must satisfy the following constraints:

$\sum\limits_{\beta }^{{}}{{{\Omega }_{\alpha \beta }}}=0\qquad \qquad( 11 )$

$\sum\limits_{b}^{{}}{{{\Omega }_{\alpha \beta }}{{\mathbf{c}}_{\beta }}}=0\qquad \qquad( 12 )$

A further assumption of the linearized collision matrix is that the local particle distribution relaxes to an equilibrium state at a single rate, τ.

${{\Omega }_{\alpha \beta }}=-\frac{1}{\tau }{{\delta }_{\alpha \beta }}\qquad \qquad( 13 )$

The delta function is unity when α = β, and zero when $\alpha \ne \beta$. This approximation is called the BGK collision term (Bhatnagar et al., 1954). With a single rate relaxation period, the Lattice Boltzmann BGK (LBGK) equation is:

${{f}_{\alpha }}\left( \mathbf{r}+{{\mathbf{c}}_{\alpha }}\Delta t,t+\Delta t \right)-{{f}_{\alpha }}\left( \mathbf{r},t \right)=-\frac{1}{\tau }\left( {{f}_{\alpha }}-\overline{{{f}_{\alpha }}} \right)\Delta t\qquad \qquad( 14 )$

Now, the relaxation period, as well as the equilibrium distribution function must be specified. The equilibrium distribution function can be found by taking a second-order expansion of the local velocity of the Maxwell-Boltzmann distribution function (Házi et al., 2002).

$\overline{{{f}_{\alpha }}}={{w}_{\alpha }}\left[ A+B{{\mathbf{c}}_{\alpha }}\cdot \mathbf{V}+C{{\left( {{\mathbf{c}}_{\alpha }}\cdot \mathbf{V} \right)}^{2}}+D{{\mathbf{V}}^{2}} \right]\qquad \qquad( 15 )$

For a two dimensional triangular grid as shown in Fig. 2 from Boltzmann Equation, the constants are:

${{w}_{\alpha }}=\rho ,\,\,A={{d}_{0}},\,\,B=C=0,\,\,D=\frac{1}{{{c}^{2}}}\text{ for }\alpha =0\qquad \qquad( 16 )$

${{w}_{\alpha }}=\rho ,\,\,A=\frac{\left( 1-{{d}_{0}} \right)}{6},\,\,B=\frac{D}{6{{c}^{2}}},\ \,C=\frac{8}{12{{c}^{4}}},\,\,D=-\frac{1}{6{{c}^{2}}}\text{ for }\alpha =1,...,6\qquad \qquad( 17 )$

where c2 is:

$\sum\limits_{\alpha }^{{}}{\left( {{\mathbf{c}}_{\alpha }}\cdot {{\mathbf{x}}_{1}} \right)\left( {{\mathbf{c}}_{\alpha }}\cdot {{\mathbf{x}}_{1}} \right)=3{{c}^{2}}}\qquad \qquad( 18 )$

where x1 is Cartesian directions, and α is a single component of a lattice vector. From these expressions, the macroscopic momentum equations can be written. The pressure and kinematic viscosity are:

$p=\frac{\rho }{2}\left( 1-{{d}_{0}} \right){{c}^{2}}\qquad \qquad( 19 )$

$\nu =\frac{{{c}^{2}}}{4}\left( \tau -\frac{1}{2} \right)\qquad \qquad( 20 )$

Therefore the relaxation time is related to the kinematic viscosity, and the constant d0 is related to the pressure. The Lattice Boltzmann equations are solved for the distribution function fα at every node. It is important to note that the Navier-Stokes equations are matched with these criteria only if the flow has a nearly constant density. Also, note that the above derivations apply only for isothermal systems, because the conservation of energy was never considered. A brief review of models that have attempted to capture thermal effects is discussed in Házi et al. (2002). Even though these models can theoretically capture arbitrary Prandtl number fluids, they are strictly limited by stability constraints.

There are two major forms of boundary conditions for hydrodynamic problems: wall boundary conditions and flow boundary conditions, in which flow enters or exits a computational domain. The no-slip wall boundary condition in LBM can be implemented by a bounce-back method. In this method, the distribution functions that arrive at a boundary are bounced back with the same distribution function in which they arrived. In a flow boundary condition, the distribution function needs to be specified during the streaming phase. This specification is driven by the macroscopic velocity or pressure. The closure to this open set of equations was achieved through the work of Zou and He (1997).

## References

Bhatnagar, P.L., Gross, E.P., and Krook, M., 1954, “A Model for Collision Processes in Gases. I Small Amplitude Processes in Charged and Neutral One-Component Systems,” Physical Review, Vol. 94, pp. 511-525.

Frisch, U, d’Humiéres, D., Hasslacher, B., Lallemand, P., Pomeau, Y. and Rivet, P.P.,1987, “Lattice Gas Hydrodynamics in Two and Three Dimensions,” Complex System, Vol. 1, pp. 649-707.

Hardy, J., Pomeau, Y., and Pazzis, O., 1973, “Time Evolution of a Two-Dimensional Model System. I. Invariant States and Time Correlation Functions,” Journal of Mathematical Physics, Vol. 14, pp. 1746-1759.

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.

McNamara, G.R. and Zanetti, G., 1988, “Use of the Boltzmann Equation to Simulate Lattice-Gas Automata,” Physical Review Letters, Vol. 61, pp. 2332-2335.

Zou, Q.S. and He, X.Y., 1997, “On pressure and velocity boundary conditions for the lattice Boltzmann BGK model,” Physics of Fluids, Vol. 9, pp. 1591-1598.