Lattice Boltzmann Theory

Overview

The Lattice Boltzmann Method (LBM) is a mesoscopic simulation technique based on the discrete Boltzmann equation. Rather than solving the Navier–Stokes equations directly, LBM evolves a set of particle distribution functions on a regular lattice. Macroscopic quantities (density, momentum, pressure) are obtained as moments of these distribution functions.

LBM is particularly well suited to:

• Complex and moving boundary geometries
• Multi-phase and multi-component flows
• Coupling with phase-field models
• Efficient parallelization

---

Discrete Boltzmann equation

The continuous Boltzmann equation with BGK (Bhatnagar–Gross–Krook) collision approximation is:

$$\begin{array}{}\text{(1)}& \frac{\partial f}{\partial t}+\mathbf{e}\cdot \mathrm{\nabla }f=-\frac{1}{\tau }(f-f^{\text{eq}})+F\end{array}$$

where:

• $f(\mathbf{x},\mathbf{e},t)$ — particle distribution function
• $\mathbf{e}$ — particle velocity
• $\tau$ — relaxation time (controls viscosity)
• $f^{\text{eq}}$ — Maxwell–Boltzmann equilibrium distribution
• $F$ — external force source term

In LBM this is discretized on a lattice with a finite set of velocities $\{{\mathbf{e}}_i\}$, yielding a set of populations $f_i(\mathbf{x},t)$.

---

D3Q27 lattice

OpenPhase uses the D3Q27 velocity set: 3 spatial dimensions, 27 discrete velocities. Each velocity ${\mathbf{e}}_i=(e_{ix},e_{iy},e_{iz})$ with components in $\{-1,0,1\}$.

The three stencils (1D, 2D, 3D) use different weight sets but are stored in the same $3\times 3\times 3$ array structure:

D1Q3 weights (1D, active dimension only):

$$\begin{array}{}\text{(2)}& w_0=\frac{4}{6},w_{\pm 1}=\frac{1}{6}\end{array}$$

D2Q9 weights (2D):

$$\begin{array}{}\text{(3)}& w_{00}=\frac{4}{9},w_{\pm 10}=w_{0\pm 1}=\frac{1}{9},w_{\pm 1\pm 1}=\frac{1}{36}\end{array}$$

D3Q27 weights (3D):

$$\begin{array}{}\text{(4)}& w_{000}=\frac{8}{27},w_{\pm 100}=\frac{2}{27},w_{\pm 1\pm 10}=\frac{1}{54},w_{\pm 1\pm 1\pm 1}=\frac{1}{216}\end{array}$$

The stencil is automatically selected during initialization based on the number of active spatial dimensions in Settings::Grid.

---

Equilibrium distribution

The local equilibrium distribution for component $n$ at lattice velocity $(x,y,z)$ is:

$$\begin{array}{}\text{(5)}& f_{xyz}^{\text{eq}}=w_{xyz}{\rho }_n[1+\frac{{\mathbf{e}}_{xyz}\cdot \mathbf{u}}{c_s^2}+\frac{({\mathbf{e}}_{xyz}\cdot \mathbf{u}{)}^2}{2c_s^4}-\frac{|\mathbf{u}{|}^2}{2c_s^2}]\end{array}$$

where:

• ${\rho }_n$ — lattice density of component $n$
• $\mathbf{u}$ — local macroscopic fluid velocity
• $c_s^2=\frac{1}{3}$ — lattice speed of sound squared (in lattice units)
• $w_{xyz}$ — D3Q27 stencil weight for direction $(x,y,z)$

---

Macroscopic quantities

Density and momentum are obtained as moments of the populations:

$$\begin{array}{}\text{(6)}& \rho =\sum _{x,y,z}{f}_{xyz}\end{array}$$

$$\begin{array}{}\text{(7)}& \rho \mathbf{u}=\sum _{x,y,z}{f}_{xyz}{\mathbf{e}}_{xyz}\end{array}$$

The lattice pressure (isothermal EOS) is:

$$\begin{array}{}\text{(8)}& p=c_s^2\rho =\frac{\rho }{3}\end{array}$$

---

BGK collision step

The BGK collision operator relaxes populations toward equilibrium:

$$\begin{array}{}\text{(9)}& f_{xyz}^{\ast }=f_{xyz}-\frac{1}{\tau }(f_{xyz}-f_{xyz}^{\text{eq}})+\mathrm{\Delta }{f}_{xyz}^{\text{force}}\end{array}$$

where $\tau$ is the relaxation time and $\mathrm{\Delta }{f}^{\text{force}}$ is the force contribution (Guo or exact-difference forcing scheme).

The kinematic viscosity is related to $\tau$ by:

$$\begin{array}{}\text{(10)}& \nu =c_s^2(\tau -\frac{1}{2})\mathrm{\Delta }t\end{array}$$

In lattice units $\mathrm{\Delta }t=1$, so choosing $\tau$ directly sets the viscosity.

---

Streaming step

After collision, populations are propagated to neighboring nodes along their respective velocity directions:

$$\begin{array}{}\text{(11)}& f_{xyz}(\mathbf{x}+{\mathbf{e}}_{xyz},t+1)=f_{xyz}^{\ast }(\mathbf{x},t)\end{array}$$

This is the exact advection step — no numerical diffusion is introduced.

---

Chapman–Enskog expansion

A Chapman–Enskog expansion of the LBM equations recovers the incompressible Navier–Stokes equations at second order:

$$\begin{array}{}\text{(12)}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }\cdot (\rho \mathbf{u})=0\end{array}$$

$$\begin{array}{}\text{(13)}& \rho (\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \mathrm{\nabla }\mathbf{u})=-\mathrm{\nabla }p+\mu {\mathrm{\nabla }}^2\mathbf{u}+\mathbf{F}\end{array}$$

with dynamic viscosity $\mu =\rho \nu$.

---

Force incorporation

The module supports two force-incorporation schemes:

Guo forcing scheme

Modifies both the equilibrium (via corrected velocity) and the collision term:

$$\begin{array}{}\text{(14)}& \mathrm{\Delta }{f}_{xyz}^{\text{force}}=w_{xyz}(1-\frac{1}{2\tau })[\frac{{\mathbf{e}}_{xyz}-\mathbf{u}}{c_s^2}+\frac{({\mathbf{e}}_{xyz}\cdot \mathbf{u}){\mathbf{e}}_{xyz}}{c_s^4}]\cdot \mathbf{F}\end{array}$$

Exact difference method (EDM)

Computes the force contribution as the difference of equilibrium distributions evaluated at shifted momenta:

$$\begin{array}{}\text{(15)}& \mathrm{\Delta }{f}_{xyz}^{\text{force}}=f_{xyz}^{\text{eq}}(\rho ,\mathbf{u}+\mathrm{\Delta }\mathbf{u})-f_{xyz}^{\text{eq}}(\rho ,\mathbf{u})\end{array}$$

where $\mathrm{\Delta }\mathbf{u}=\mathbf{F}\mathrm{\Delta }t/\rho$.

EDM conserves momentum exactly and is preferred for multi-phase flows.

---

Lattice units and physical scaling

LBM operates in dimensionless lattice units. Physical quantities are recovered via scaling coefficients defined in FlowSolverLBM:

Field	Coefficient	Unit
Density	dRho	kg/m³
Pressure	dP	Pa
Momentum density	dm	kg/(m² s)
Force density	df	kg/(m² s²)
Kinematic viscosity	dnu	m²/s
Time	dt	s

Proper choice of these coefficients ensures the simulation is in the physically relevant regime (low Mach number for incompressible flow).

---

Stability considerations

• The relaxation time $\tau$ must satisfy $\tau >0.5$ for stability.
• Typical values: $\tau \in [0.55,2.0]$.
• Very low viscosity ($\tau \to 0.5$) leads to instability.
• Use Do_FixPopulations = true to correct negative populations that can arise near interfaces or obstacles.
• For two-phase flows, the density ratio ${\rho }_{\text{liquid}}/{\rho }_{\text{vapor}}$ is limited by the pseudo-potential model; excessively high ratios cause instability.