Search

[LBM Study] 1. The Lattice Boltzmann Equation

Category
Study
Kewords
LBM
Kruger
3 more properties
Table of Contents

1.1 Introduction

•
If we can solve the Boltzmann equation numerically, this may also indirectly give us a solution to the Navier-Stokes Equation (NSE).
•
The numerical scheme for the Boltzmann equation somewhat paradoxically turns out to be quite simple, both to implement and to parallelise.
◦
A simple hyperbolic equation which essentially describes the advection of the distribution function ff with the particle velocity ξ\xi.
◦
The ource term Ω(f)\Omega(f) depends only on the local value of ff and not on its gradients.

1.2 The Lattice Boltzmann Equation in a Nutshell

1.2.1 Overview

•
The basic quantity of the Lattice Boltzmann Method (LBM) is the discrete-velocity distribution function fi(x,t)f_{i}(\mathbf{x},t).
◦
It represents the density of particles with velocity ci=(cix,ciy,ciz)\bold{c}_{i} = (c_{ix}, c_{iy}, c_{iz}) at position x\bold{x} and time tt.
◦
Mass density ρ\rho and momentum density ρu\rho\bold{u} can be expressed as:
ρ(x,t)=∑ifi(x,t)  ,ρu(x,t)=∑icifi(x,t).\rho(\mathbf{x},t) = \sum_{i}f_{i}(\mathbf{x},t)\;,\qquad\rho\mathbf{u}(\mathbf{x},t) = \sum_{i}\mathbf{c}_{i}f_{i}(\mathbf{x},t).
•
The time step Δt\Delta t and lattice spacing Δx\Delta x respectively represent a time resolution and a space resolution in any set of units.
◦
The most common choice in the LB (Lattice Boltzmann) literature is lattice units → Δt=1\Delta t=1 and Δx=1\Delta x=1.
◦
We can convert quantities between lattice units and physical units. → Law of Similarity: only ensure that the relevant dimensionless numbers.
•
The discrte velocities ci\bold{c}_{i} and weighting coefficients wiw_{i} form velocity sets {ci,wi}\{ \bold{c}_{i}, w_{i} \}.
◦
Different velocity sets are used for different purposes.
◦
These velocity sets are usually denoted by DdQqDdQq, where dd is the number of spatial dimensions the velocity set covers and qq is the set’s number of velocities.
◦
Commonly used velocity sets: D1Q3D1Q3, D2Q9D2Q9, D3Q15D3Q15, D3Q19D3Q19, and D3Q27D3Q27. → In 3D, the most commonly used velocity set is D3Q19D3Q19.
•
In the basic isothermal LBE (Lattice Boltzmann Equation), csc_{s} determines the relation p=cs2ρp = c_{s}^{2} \rho between pressure pp and density ρ\rho.
◦
Because of this relation, it can be shown that csc_{s} represents the isothermal model’s speed of sound.
◦
In all the velocity sets, this constant is cs2=(1/3)Δx2/Δt2c_{s}^{2}=(1/3) \Delta x^2 / \Delta t^2.
•
By discretising the Boltzmann equation in velocity space, physical space, and time, we find the LBE:
fi(x+ciΔt,t+Δt)=fi(x,t)+Ωi(x,t).f_{i}(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f_{i}(\mathbf{x}, t) + \Omega_{i}(\mathbf{x}, t).
◦
This expresses that particles fi(x,t)f_{i}(\bold{x}, t) move with velocity ci\bold{c}_{i} to a neighbouring point x+ciΔt\bold{x} + \bold{c}_{i} \Delta t at the next time step t+Δtt + \Delta t.
◦
At the same time, particles are affected by a collision operator Ωi\Omega_{i}.
•
The simplest collision operator that can be used for Navier-Stokes simulations is the BGK (Bhatnagar-Gross-Krook) operator:
Ωi(f)=−fi−fieqτΔt.\Omega_{i}(f)=-\frac{f_{i} - f_{i}^\mathrm{eq}}{\tau} \Delta t.
◦
It relaxes the populations towards an equilibrium fieqf_{i}^\mathrm{eq} at a rate determined by the relaxation time τ\tau.
◦
This equilibrium is given by:
fieq(x,t)=wi ρ (1+u⋅cics2+(u⋅ci)22 cs4−u⋅u2 cs2).f_i^{\mathrm{eq}}(\mathbf{x},t)= w_i\,\rho\,\biggl( 1 + \frac{\mathbf{u}\cdot\mathbf{c}_i}{c_s^2} + \frac{(\mathbf{u}\cdot\mathbf{c}_i)^2}{2\,c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2\,c_s^2}\biggr).
◦
The equilibrium is such that its moments are the same as those of fif_{i}, i.e. ∑ifieq=∑ifi=ρ\sum_{i} f_{i}^\mathrm{eq} = \sum_{i} f_{i}=\rho and ∑icifieq=∑icifi=ρu\sum_{i} \mathbf{c}_{i} f_{i}^\mathrm{eq} = \sum_{i} \mathbf{c}_{i} f_{i}=\rho \mathbf{u}.
◦
The equilibrium fieqf_{i}^\mathrm{eq} depends on the local quantities density ρ\rho and fluid velocity u\mathbf{u} only.
•
The link between the LBE and the NSE can be determined using the Chapman-Enskog analysis.
◦
The LBE results in macroscopic behaviour according to the NSE, with the kinematic shear viscosity ν\nu given by the relaxation time τ\tau as:
ν=cs2(τ−Δt2).\nu=c_{s}^{2} \biggl(\tau-\frac{\Delta t}{2} \biggl).
◦
The kinematic bulk viscosity given as νB=2ν/3.\nu_{B} = 2 \nu / 3.

1.2.2 The Time Step: Collision and Streaming

•
The Lattice BGK (LBGK) equation reads:
fi(x+ci Δt,  t+Δt)=fi(x,t)−Δtτ(fi(x,t)−fieq(x,t)).f_i\bigl(\mathbf{x} + \mathbf{c}_i\,\Delta t,\; t + \Delta t\bigr)= f_i(\mathbf{x}, t)- \frac{\Delta t}{\tau}\bigl(f_i(\mathbf{x}, t) - f_i^{\mathrm{eq}}(\mathbf{x}, t)\bigr).
◦
We can decompose this equation into two distinct parts: Collision (or Relaxation) & Streaming (or Propagation).
•
The first part is collision:
fi⋆(x,t)=fi(x,t) (1−Δtτ)+fieq(x,t) Δtτ .f_i^{\star}(\mathbf{x},t) = f_i(\mathbf{x},t)\,\Bigl(1 - \frac{\Delta t}{\tau}\Bigr) + f_i^{\mathrm{eq}}(\mathbf{x},t)\,\frac{\Delta t}{\tau}\,.
◦
fi⋆f_{i}^{\star} represents the distribution function after collisions
◦
fieqf_{i}^\mathrm{eq} is found from:
fieq(x,t)=wi ρ (1+u⋅cics2+(u⋅ci)22 cs4−u⋅u2 cs2).f_i^{\mathrm{eq}}(\mathbf{x},t)= w_i\,\rho\,\biggl( 1 + \frac{\mathbf{u}\cdot\mathbf{c}_i}{c_s^2} + \frac{(\mathbf{u}\cdot\mathbf{c}_i)^2}{2\,c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2\,c_s^2}\biggr).
◦
If τ/Δt=1\tau / \Delta t = 1 → fi⋆(x,t)=fieq(x,t).f_i^{\star}(\mathbf{x},t) = f_i^{\mathrm{eq}}(\mathbf{x},t).
•
The second part is streaming:
fi(x+ci Δt,t+Δt)=fi⋆(x,t).f_i\bigl(\mathbf{x} + \mathbf{c}_i\,\Delta t,t + \Delta t\bigr) = f_i^{\star}(\mathbf{x},t).
•
Overall, the LBE concept is straightforward.
◦
Collision: one calculates the density ρ\rho and the macroscopic velocity u\mathbf{u} to find the equilibrium distributions fieqf_i^{\mathrm{eq}} and the post-collision distribution fi⋆f_i^{\star}.
◦
Streaming: we stream the resulting distribution fi⋆f_i^{\star} to neighbouring nodes.

1.3 Implementation of the Lattice Boltzmann Method in a Nutshell

1.3.1 Initialization

•
The simplest approach to initializing the populations at the start of a simulation is to set them to:
fieq(x, t=0)=fieq(ρ(x, t=0), u(x, t=0))f_i^{\mathrm{eq}}\bigl(\mathbf{x},\,t=0\bigr) = f_i^{\mathrm{eq}}\bigl(\rho(\mathbf{x},\,t=0),\,\mathbf{u}(\mathbf{x},\,t=0)\bigr)
via:
fieq(x,t)=wi ρ (1+u⋅cics2+(u⋅ci)22 cs4−u⋅u2 cs2).f_i^{\mathrm{eq}}(\mathbf{x},t)= w_i\,\rho\,\biggl( 1 + \frac{\mathbf{u}\cdot\mathbf{c}_i}{c_s^2} + \frac{(\mathbf{u}\cdot\mathbf{c}_i)^2}{2\,c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2\,c_s^2}\biggr).
•
Often, the values ρ(x,t=0)=1\rho(\mathbf{x},t=0)=1 and u(x, t=0)=0\mathbf{u}(\mathbf{x},\,t=0)=\mathbf{0} are used.

1.3.2 Time Step Algorithm

•
The core LBM algorithm consists of a cyclic sequence of substeps, with each cycle corresponding to one time step.
◦
The dark grey boxes show sub-steps that are necessary for the evolution of the solution.
◦
The light grey box indicates the optional output step.
◦
The pale boxes indicate steps with details provided in later sections (Boundary Conditions & Forces).
1.
Compute the macroscopic moments ρ(x,t)\rho(\mathbf{x},t) and u(x,t)\mathbf{u}(\mathbf{x},t) from fi(x,t)f_i(\mathbf{x},t) via:
ρ(x,t)=∑ifi(x,t)  ,ρu(x,t)=∑icifi(x,t).\rho(\mathbf{x},t) = \sum_{i}f_{i}(\mathbf{x},t)\;,\qquad\rho\mathbf{u}(\mathbf{x},t) = \sum_{i}\mathbf{c}_{i}f_{i}(\mathbf{x},t).
2.
Obtain the equilibrium distribution fieq(x,t)f_i^{\mathrm{eq}}(\mathbf{x},t) from:
fieq(x,t)=wi ρ (1+u⋅cics2+(u⋅ci)22 cs4−u⋅u2 cs2).f_i^{\mathrm{eq}}(\mathbf{x},t)= w_i\,\rho\,\biggl( 1 + \frac{\mathbf{u}\cdot\mathbf{c}_i}{c_s^2} + \frac{(\mathbf{u}\cdot\mathbf{c}_i)^2}{2\,c_s^4} - \frac{\mathbf{u}\cdot\mathbf{u}}{2\,c_s^2}\biggr).
3.
If desired, write the macroscopic fields ρ(x,t)\rho(\mathbf{x},t), u(x,t)\mathbf{u}(\mathbf{x},t), and/or σ(x,t)\mathbf{\sigma}(\mathbf{x},t) to the hard disk for visualisation or post-processing. The viscous stress tensor σ\mathbf{\sigma} can be computed from:
σαβ≈−(1−Δt2τ)∑iciα ciβ fineq.\sigma_{\alpha\beta} \approx - \Bigl(1 - \frac{\Delta t}{2\tau}\Bigr) \sum_i c_{i\alpha}\,c_{i\beta}\,f_i^{\mathrm{neq}} .
4.
Perform collision (relaxation) as shown in:
fi⋆(x,t)=fi(x,t) (1−Δtτ)+fieq(x,t) Δtτ .f_i^{\star}(\mathbf{x},t) = f_i(\mathbf{x},t)\,\Bigl(1 - \frac{\Delta t}{\tau}\Bigr) + f_i^{\mathrm{eq}}(\mathbf{x},t)\,\frac{\Delta t}{\tau}\,.
5.
Perform streaming (propagation) via:
fi(x+ci Δt,t+Δt)=fi⋆(x,t).f_i\bigl(\mathbf{x} + \mathbf{c}_i\,\Delta t,t + \Delta t\bigr) = f_i^{\star}(\mathbf{x},t).
6.
Increase the time step, setting tt to t+Δtt + \Delta t, and go back to step 1 until the last time step or convergence has been reached.

1.3.3 Notes on Memory Layout and Coding Hints

•
Initialisation
◦
In a 2D simulation,
▪
Macroscopic fields (2D arrays) → ρ[Nx][Ny]\rho [N_x][N_y], ux[Nx][Ny]u_x [N_x][N_y], uy[Nx][Ny]u_y [N_x][N_y]
▪
Populations (3D array) → f[Nx][Ny][q]f[N_x][N_y][q]
▪
NxN_x and NyN_y are the number of lattice nodes in x- and y-directions and qq is the number of velocities.
•
Streaming
◦
The streaming step must be implemented carefully to prevent overwriting population data before it has been streamed from the source node.
◦
Three common strategies exist to handle this:
1.
Opposite-Direction Sweep: Use a temporary buffer for a single population while sweeping through memory in the direction opposite to streaming.
•
Helper functions like circshift (Matlab), numpy.roll (Python), and CSHIFT (Fortran) can simplify this approach.
2.
Two-Array (Swap) Method: Allocate two full population arrays (fioldf_i^{\mathrm{old}}, finewf_i^{\mathrm{new}}).
•
Read from fiold(x)f_i^{\mathrm{old}}(\mathbf{x}) and write the streamed populations to finew(x+ciΔt)f_i^{\mathrm{new}}(\mathbf{x} + \mathbf{c}_i \Delta t).
3.
Combined Stream-and-Collide: Perform streaming and collision in a single, fused step.
•
This reads required populations from adjacent nodes directly into the collision calculation at the current node.
•
Updating Macroscopic Variables
◦
Instead of using loops, unroll the summations for calculating macroscopic moments like density ρ\rho and velocity u\mathbf{u}. → Increase the CPU’s computing efficiency.
◦
For the D2Q9D2Q9 velocity set, the unrolled implementations would be:
ρ=f0+f1+f2+f3+f4+f5+f6+f7+f8,ux=[(f1+f5+f8)−(f3+f6+f7)]/ρ,uy=[(f2+f5+f6)−(f4+f7+f8)]/ρ.ρ = f_0 + f_1 + f_2 + f_3 + f_4 + f_5 + f_6 + f_7 + f_8, \\\\ u_x = [(f_1 + f_5 + f_8) - (f_3 + f_6 + f_7)] / ρ, \\\\ u_y = [(f_2 + f_5 + f_6) - (f_4 + f_7 + f_8)] / ρ.
•
Equilibrium
◦
It is recommended to unroll the loops for the equilibrium computation ( fieqf_i^{\mathrm{eq}}), writing separate expressions for each direction ii of the qq distributions.
◦
Pre-calculate constants to avoid expensive divisions inside the main loop ( cs−2c_s^{-2}, cs−4c_s^{-4}).
•
Collision
◦
To optimize the collision step fi∗=(1−Δt/τ)fi+(Δt/τ)fieqf_{i}^{*} = (1 - Δt/τ)f_{i} + (Δt/τ)f_{i}^{\mathrm{eq}}, pre-calculate the relaxation frequency.
▪
Define constants like ω=Δt/τω = Δt/τ and ω′=1−ωω' = 1 - ω once before the main loop.
▪
The collision step then becomes a more efficient calculation: fi∗=ω′fi+ωfieqf_{i}^{*} = ω'f_{i} + ωf_{i}^{\mathrm{eq}}.

1.4 Discretization in Velocity Space

1.4.1 Non-Dimensionalization

•
To simplify the derivation, the governing equations are non-dimensionalized.
•
The continuous and force-free Boltzmann equation in its non-dimensional form is:
∂f∂t+ξα∂f∂xα=Ω(f).\frac{\partial f}{\partial t} + \xi_α \frac{\partial f}{\partial x_α} = \Omega(f).
◦
ff: Particle distribution function.
◦
xαx_{α}: Components of the spatial coordinate vector x\mathbf{x}.
◦
ξα\xi_{α}: Components of the particle velocity vector ξ\mathbf{ξ}.
◦
Ω(f)\Omega(f): Collision operator.
•
The corresponding non-dimensional equilibrium distribution function is:
feq(ρ,u,θ,ξ)=ρ(2πθ)d/2e−(ξ−u)2/(2θ).f^{\mathrm{eq}}(ρ, \mathbf{u}, θ, \mathbf{ξ}) = \frac{ρ}{(2πθ)^{d/2}} e^{-(\mathbf{ξ}-\mathbf{u})^2 / (2θ)}.
◦
feqf^{\mathrm{eq}}: Equilibrium distribution function.
◦
θ\theta: Non-dimensional temperature.
◦
dd: Number of spatial dimensions.

1.4.2 Hermite Polynomials

•
Hermite Polynomials (HPs) are used for the discretization of integrals and form the mathematical basis of the Hermite series expansion.
•
They are a set of orthogonal polynomials constructed from a Gaussian weight function ω(x)\omega(\mathbf{x}).
•
The multidimensional HP of rank nn is defined as:
H(n)(x)=(−1)n1ω(x)∇(n)ω(x),ω(x)=1(2π)d/2e−x2/2.\mathbf{H}^{(n)}(\mathbf{x}) = (-1)^n \frac{1}{ω(\mathbf{x})} \mathbf{∇}^{(n)}ω(\mathbf{x}), \quad \quad ω(\mathbf{x}) = \frac{1}{(2π)^{d/2}} e^{-\mathbf{x}^2/2}.
◦
H(n)(x)\mathbf{H}^{(n)}(\mathbf{x}): Hermite Polynomial tensor of rank nn.
◦
ω(x)\omega(\mathbf{x}): Gaussian weight function.
◦
∇(n)\mathbf{∇}^{(n)}: nn th order gradient operator.

1.4.3 Hermite Series Expansion of the Equilibrium Distribution

•
The equilibrium distribution function feq(ξ)f^{\mathrm{eq}}(\mathbf{ξ}) is expanded into an infinite Hermite series.
◦
A key insight is that the expansion coefficients are directly related to the macroscopic moments of the fluid (density, momentum, energy, etc.).
•
The equilibrium distribution function feq(ξ)f^{\mathrm{eq}}(\mathbf{ξ}) has the same mathematical form as the Hermite weight function ω(ξ)ω(\mathbf{ξ}), allowing it to be expressed as:
feq(ρ,u,θ,ξ)=ρ(2πθ)d/2e−(ξ−u)2/(2θ)=ρθd/2ω(ξ−uθ).f^{\mathrm{eq}}(ρ, \mathbf{u}, θ, \mathbf{ξ}) = \frac{ρ}{(2πθ)^{d/2}}e^{{-(\mathbf{ξ}-\mathbf{u})^2} / {(2θ)}} = \frac{ρ}{θ^{d/2}} ω\left(\frac{\mathbf{ξ}-\mathbf{u}}{\sqrt{θ}}\right).
◦
feqf^{\mathrm{eq}}: The equilibrium distribution function.
◦
ρ\rho: Macroscopic density.
◦
u\mathbf{u}: Macroscopic velocity.
◦
θθ: Non-dimensional temperature.
◦
ξ\mathbf{ξ}: Particle velocity vector.
◦
dd: Number of spatial dimensions.
◦
ω\omega: The Hermite weight function (a Gaussian).
•
Key Conclusion: To recover the correct macroscopic conservation laws for hydrodynamics, it is sufficient to truncate the Hermite series expansion after the second-order term.
◦
This is a crucial simplification that makes the LBM computationally feasible.

1.4.4 Discretisation of the Equilibrium Distribution Function

•
The Gauss-Hermite quadrature rule is used to replace continuous integrals over velocity space with discrete, weighted sums. ⇒ This rule allows for the exact integration of polynomials by sampling them at a few specific points called abscissae ξi\mathbf{\xi}_i.
•
To obtain velocity sets with simple integer components, the abscissae are rescaled to define the discrete particle velocities ci\mathbf{c}_i:
ci=ξi3.\mathbf{c}_i = \frac{\mathbf{ξ}_i}{ \sqrt{3}}.
•
This procedure yields the final form of the discrete equilibrium distribution function. This equation is a cornerstone of the LBM, shown here for an isothermal (θ=1θ=1) system:
fieq=wiρ(1+ci⋅ucs2+(ci⋅u)22cs4−u⋅u2cs2)=wiρ(1+ciαuαcs2+uαuβ(ciαciβ−cs2δαβ)2cs4).f_i^{\mathrm{eq}} = w_i ρ \left( 1 + \frac{\mathbf{c}_i \cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i \cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u} \cdot \mathbf{u}}{2c_s^2} \right) = w_i ρ \left( 1 + \frac{c_{iα} u_α}{c_s^2} + \frac{u_α u_β (c_{iα} c_{iβ} - c_s^2 δ_{αβ})}{2c_s^4} \right).
◦
The Greek letters α\alpha and β\beta are indices representing the Cartesian coordinate components ( xx, yy, zz).
▪
The Einstein summation convention is applied. → ciαuα=∑αciαuα=ci⋅uc_{iα} u_α = \sum_α c_{iα} u_α = \mathbf{c}_i \cdot \mathbf{u} (the vector dot product).
▪
δαβδ_{αβ} is the Kronecker delta, which is 1 if α=βα=β and 0 if α≠βα \neq β.
◦
fieqf_i^{\mathrm{eq}}: The discrete equilibrium distribution for the ii th velocity direction.
◦
wiw_i: The lattice weight associated with direction ii.
◦
ρ\rho: The macroscopic fluid density.
◦
u\mathbf{u}: The macroscopic fluid velocity vector (uαu_\alpha and uβu_\beta are its components).
◦
ci\mathbf{c}_i: The discrete velocity vector for direction ii (ciαc_{i\alpha} and ciβc_{i\beta} are its components).
◦
csc_s: The lattice speed of sound, a constant determined by the velocity set.

1.4.5 Discretisation of the Particle Distribution Function

•
The same discretization process used for the equilibrium function is now applied to the non-equilibrium particle distribution function f(x,ξ,t)f(\mathbf{x}, \mathbf{ξ}, t).
◦
This projects the continuous function onto the discrete velocity set {ci}\{ \mathbf{c}_i \}, resulting in a finite set of qq populations fi(x,t)f_i(\mathbf{x}, t).
◦
These populations become the fundamental variables of the LBM.
•
The evolution of these discrete populations is governed by the discrete-velocity Boltzmann equation:
∂tfi+ciα∂αfi=Ωi(fi) ,     i=1,...,q.\partial_t f_i + c_{iα}\partial_α f_i = Ω_i(f_i)\:, \:\:\:\:\: i=1,...,q.
◦
fi(x,t)f_i(\mathbf{x}, t): The discrete particle distribution (or population) for direction ii at position x\mathbf{x} and time tt.
◦
ciαc_{iα}: The αα component of the discrete velocity ci\mathbf{c}_i. The repeated αα implies summation over all spatial components.
◦
Ωi(fi)Ω_i(f_i): The discrete collision operator acting on the ii th population.
•
A major consequence of this discretization is that macroscopic quantities are now computed through simple, computationally efficient summations over the discrete populations:
ρ=∑ifi=∑ifieq, ρu=∑ifici=∑ifieqci.ρ = \sum_i f_i = \sum_i f_i^{\mathrm{eq}}, \\\\ \: \\\\ ρ\mathbf{u} = \sum_i f_i \mathbf{c}_i = \sum_i f_i^{\mathrm{eq}} \mathbf{c}_i.

1.4.6 Discretisation of the Particle Distribution Function

•
A key component of the Lattice Boltzmann Method is the choice of the discrete velocity set {ci}\{ \mathbf{c}_i \} and its corresponding weights wiw_i.
•
General Comments and Definitions
◦
Velocity sets are commonly named using the DdQqDdQq notation, where dd is the number of spatial dimensions and qq is the number of discrete velocities.
◦
Implementation Warning: The numbering of velocity vectors (”i=0...q−1i=0...q-1” vs. ”i=1...qi=1...q”) is not consistent across the literature.
◦
Here follows the convention of using index i=0i=0 for the rest velocity (c0=0\mathbf{c}_0 = \mathbf{0}) if one exists.
•
Construction and Requirements of Velocity Sets
◦
To correctly solve the NSE, a velocity set and its weights must satisfy certain isotropy conditions.
◦
The key conditions are:
∑iwi=1,∑iwiciα=0,∑iwiciαciβ=cs2δαβ,∑iwiciαciβciγ=0,∑iwiciαciβciγciδ=cs4(δαβδγδ+δαγδβδ+δαδδβγ),∑iwiciαciβciγciδciε=0.\begin{aligned} \sum_i w_i = 1, \\ \sum_i w_i c_{iα} = 0, \\ \sum_i w_i c_{iα} c_{iβ} = c_s^2 δ_{αβ}, \\ \sum_i w_i c_{iα} c_{iβ} c_{iγ} = 0, \\ \sum_i w_i c_{iα} c_{iβ} c_{iγ} c_{iδ} = c_s^4 (δ_{αβ}δ_{γδ} + δ_{αγ}δ_{βδ} + δ_{αδ}δ_{βγ}), \\ \sum_i w_i c_{iα} c_{iβ} c_{iγ} c_{iδ} c_{iε} = 0. \end{aligned}
◦
Additional Constraints for Standard LBM:
▪
All weights must be non-negative: wi≥0w_i ≥ 0.
▪
For simple streaming on a uniform grid, the velocity components are typically integer multiples of Δx/Δt\Delta x / \Delta t. ⇒ In lattice units where Δx=1\Delta x=1 and Δt=1\Delta t=1, the velocity vectors ci\mathbf{c}_i have integer components.
•
Common Velocity Sets for Hydrodynamics
◦
The most common velocity sets used for hydrodynamic simulations are D1Q3D1Q3, D2Q9D2Q9, D3Q15D3Q15, D3Q19D3Q19, and D3Q27D3Q27.
◦
Practical Guidance for 3D Simulations:
▪
D3Q19D3Q19 is generally a good compromise for simulating laminar flows.
▪
D3Q27D3Q27 offers better rotational invariance and is preferred for high Reynolds number or turbulent flows, although it is more computationally expensive.
•
Equilibrium Distributions
◦
To improve computational performance, it is highly recommended to use the "unrolled" form of the equilibrium equation for each direction rather than calculating it in a loop.
◦
Example for D2Q9D2Q9 (using u2=ux2+uy2\mathbf{u^2} = u_x^2 + u_y^2 and cs2=1/3c_s^2 = 1/3 in lattice units):
f0eq=2ρ9(2−3u2),f1eq=ρ18(2+6ux+9ux2−3u2),f5eq=ρ36[1+3(ux+uy)+9uxuy+3u2],f2eq=ρ18(2+6uy+9uy2−3u2),f6eq=ρ36[1−3(ux−uy)−9uxuy+3u2],f3eq=ρ18(2−6ux+9ux2−3u2),f7eq=ρ36[1−3(ux+uy)+9uxuy+3u2],f4eq=ρ18(2−6uy+9uy2−3u2),f8eq=ρ36[1+3(ux−uy)−9uxuy+3u2]. \begin{aligned} f_0^{\mathrm{eq}} &= \frac{2ρ}{9}(2 - 3\mathbf{u^2}), \\ f_1^{\mathrm{eq}} &= \frac{ρ}{18}(2 + 6u_x + 9u_x^2 - 3\mathbf{u^2}), \quad & f_5^{\mathrm{eq}} &= \frac{ρ}{36}[1 + 3(u_x + u_y) + 9u_xu_y + 3\mathbf{u^2}], \\ f_2^{\mathrm{eq}} &= \frac{ρ}{18}(2 + 6u_y + 9u_y^2 - 3\mathbf{u^2}), \quad & f_6^{\mathrm{eq}} &= \frac{ρ}{36}[1 - 3(u_x - u_y) - 9u_xu_y + 3\mathbf{u^2}], \\ f_3^{\mathrm{eq}} &= \frac{ρ}{18}(2 - 6u_x + 9u_x^2 - 3\mathbf{u^2}), \quad & f_7^{\mathrm{eq}} &= \frac{ρ}{36}[1 - 3(u_x + u_y) + 9u_xu_y + 3\mathbf{u^2}], \\ f_4^{\mathrm{eq}} &= \frac{ρ}{18}(2 - 6u_y + 9u_y^2 - 3\mathbf{u^2}), \quad & f_8^{\mathrm{eq}} &= \frac{ρ}{36}[1 + 3(u_x - u_y) - 9u_xu_y + 3\mathbf{u^2}]. \end{aligned}
•
Macroscopic Moments
◦
The velocity sets are designed such that the moments of the equilibrium distribution can be expressed simply in terms of macroscopic variables. ⇒ These relations are fundamental to the Chapman-Enskog analysis that proves the LBM-to-NSE correspondence.
◦
The key equilibrium moments are:
Πeq=∑ifieq=ρ,Παeq=∑ifieqciα=ρuα,Παβeq=∑ifieqciαciβ=ρcs2δαβ+ρuαuβ,Παβγeq=∑ifieqciαciβciγ=ρcs2(uαδβγ+uβδαγ+uγδαβ).\begin{aligned} \Pi^{\mathrm{eq}} &= \sum_i f_i^{\mathrm{eq}} = ρ, \\ \Pi_{α}^{\mathrm{eq}} &= \sum_i f_i^{\mathrm{eq}} c_{iα} = ρu_α, \\ \Pi_{αβ}^{\mathrm{eq}} &= \sum_i f_i^{\mathrm{eq}} c_{iα} c_{iβ} = ρc_s^2 δ_{αβ} + ρu_α u_β, \\ \Pi_{αβγ}^{\mathrm{eq}} &= \sum_i f_i^{\mathrm{eq}} c_{iα} c_{iβ} c_{iγ} = ρc_s^2(u_α δ_{βγ} + u_β δ_{αγ} + u_γ δ_{αβ}). \end{aligned}

1.5 Discretization in Space and Time

1.5.1 Method of Characteristics

•
Core Idea:
◦
The method transforms the partial differential equation (PDE) for the discrete populations fif_i into a simpler ordinary differential equation (ODE) that is valid along specific trajectories in spacetime, known as "characteristics."
•
Starting Equation:
∂tfi+ciα∂αfi=Ωi.\partial_t f_i + c_{iα}\partial_α f_i = Ω_i.
◦
The starting point is the discrete-velocity Boltzmann equation.
◦
fi(x,t)f_i(\mathbf{x}, t): The discrete population for direction ii.
◦
ci\mathbf{c}_i: The discrete velocity vector for direction ii (ciαc_{iα} is its αα component).
◦
ΩiΩ_i: The discrete collision operator.
•
The Transformation to an ODE (The Characteristic Trajectory):
◦
The advection operator on the left ∂t+ci⋅∇\partial_t + \mathbf{c}_i \cdot \nabla, is converted into a total derivative dfi/ζdf_i/ \zeta along a characteristic trajectory defined by:
dtdζ=1,dxαdζ=ciα.\frac{dt}{d \zeta} = 1, \quad \frac{d\mathbf{x}_\alpha}{d \zeta} = {c}_{i \alpha}.
◦
This simplifies the original PDE to an ODE along this trajectory:
dfidζ=(∂fi∂t)dtdζ+(∂fi∂xα)dxαdζ=Ωi(x(ζ),t(ζ)).\frac{df_i}{d\zeta} = \left( \frac{\partial f_i}{\partial t} \right) \frac{dt}{d\zeta} + \left( \frac{\partial f_i}{\partial {x}_\alpha} \right) \frac{d\mathbf{x}_\alpha}{d\zeta} = \Omega_i(\mathbf{x}(\zeta), t(\zeta)).
•
The Integral Form (Key Result):
◦
By integrating the ODE along the characteristic trajectory over a single time step from ζ=0\zeta=0 to ζ=Δt\zeta=\Delta t, we obtain the integral form of the LBE.
◦
The integration of the left-hand side is exact, according to the fundamental theorem of calculus.
◦
The resulting equation, generalized for any arbitrary starting point (x,t)(\mathbf{x}, t) is:
fi(x+ciΔt,t+Δt)−fi(x,t)=∫0ΔtΩi(x+ciζ,t+ζ)dζ.f_i(\mathbf{x} + \mathbf{c}i \Delta t, t + \Delta t) - f_i(\mathbf{x}, t) = \int_{0}^{\Delta t} \Omega_i(\mathbf{x} + \mathbf{c}_i \zeta, t + \zeta) d\zeta.
◦
This equation is central because it reveals the discretization pattern:
1.
Left-hand side: An exact, discrete update rule for propagation, showing that during a time step Δt\Delta t, the population fi(x,t)f_i(\mathbf{x}, t) moves from x\mathbf{x} to x+ciΔt\mathbf{x} + \mathbf{c}_i\Delta t.
2.
Right-hand side: An integral of the collision term along the trajectory, which must be approximated in subsequent steps.

1.5.2 First- and Second-Order Discretization

•
First-Order Discretization
◦
The simplest approximation uses a first-order (rectangular) rule, which evaluates the collision operator Ωi\Omega_i at the beginning of the time step tt.
◦
This approximation of the integral ∫0ΔtΩidζ≈ΔtΩi(x,t)\int_{0}^{\Delta t} \Omega_i d\zeta \approx \Delta t\Omega_i(\mathbf{x}, t) results in the explicit discretized LBE:
fi(x+ciΔt,t+Δt)=fi(x,t)+ΔtΩi(x,t).f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f_i(\mathbf{x}, t) + \Delta t\Omega_i(\mathbf{x}, t).
•
Second-Order Discretization
◦
A more accurate approximation uses the second-order trapezoidal rule to evaluate the collision integral:
∫0ΔtΩidζ≈Δt2[Ωi(x,t)+Ωi(x+ciΔt,t+Δt)].\int_{0}^{\Delta t} \Omega_i d\zeta \approx \frac{\Delta t}{2} \left[ \Omega_i(\mathbf{x}, t) + \Omega_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) \right].
◦
This scheme is initially implicit because the collision term Ωi\Omega_i at time t+Δtt + \Delta t depends on the unknown populations fif_i at that same time.
◦
However, through a specific change of variables (fi→f~if_i \rightarrow \tilde{f}_i), this implicit equation can be transformed into an explicit equation that has the same form as the one derived from the first-order rule.
•
The fact that the second-order accurate trapezoidal scheme can be rewritten into the same form as the first-order scheme proves that the standard LBE is second-order accurate in time.

1.5.3 BGK Collision Operator

•
To complete the LBE, the collision operator Ωi\Omega_i must be specified. ⇒ The simplest and most widely used choice is the Bhatnagar-Gross-Krook (BGK) model, also known as the single-relaxation-time (SRT) model.
•
Definition
◦
The BGK operator models the complex process of particle collisions as a simple linear relaxation of each population fif_i towards its corresponding local equilibrium state fieqf_i^{\mathrm{eq}}:
Ωi=−fi−fieqτ.\Omega_i = -\frac{f_i - f_i^{\mathrm{eq}}}{\tau}.
◦
fif_i: The discrete particle population for direction ii.
◦
fieqf_i^{\mathrm{eq}}: The discrete equilibrium population.
◦
τ\tau: A characteristic time scale known as the relaxation time, which controls the rate of equilibration.
•
The Lattice BGK (LBGK) Equation
◦
Substituting the BGK operator into the discretized LBE (first-order discretization) yields the full Lattice BGK (LBGK) equation:
fi(x+ciΔt,t+Δt)=fi(x,t)−Δtτ(fi(x,t)−fieq(x,t)).f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f_i(\mathbf{x}, t) - \frac{\Delta t}{\tau} \left( f_i(\mathbf{x}, t) - f_i^{\mathrm{eq}}(\mathbf{x}, t) \right).
◦
This simple form is capable of reproducing the full Navier-Stokes equations, which is a primary reason for the popularity of the LBM.
◦
Note on Accuracy and Other Operators:
▪
Although the LBGK equation above can be derived from a first-order time integration, it is actually second-order accurate in time.
▪
More advanced operators like the Two-Relaxation-Time (TRT) and Multiple-Relaxation-Time (MRT) models exist to overcome some of its limitations.
•
Relaxation Regimes and Stability
◦
The behavior of the discrete LBGK equation depends on the dimensionless ratio τ/Δt\tau / \Delta t.
▪
Under-Relaxation (τ/Δt>1\tau / \Delta t > 1): fif_i decays exponentially towards fieqf_i^{\mathrm{eq}}, similar to the continuous-time case.
▪
Full-Relaxation (τ/Δt=1\tau / \Delta t = 1): fif_i relaxes directly to fieqf_i^{\mathrm{eq}} in a single step.
▪
Over-Relaxation (1/2<τ/Δt<11/2 < \tau / \Delta t < 1): fif_i oscillates around fieqf_i^{\mathrm{eq}} with an exponentially decreasing amplitude.
◦
This analysis reveals a necessary condition for the numerical stability of the LBGK model:
τΔt≥12.\frac{\tau}{\Delta t} \ge \frac{1}{2}.
▪
If τ/Δt<1/2\tau / \Delta t < 1/2, the populations will oscillate with an exponentially increasing amplitude, leading to instability.

1.5.4 Streaming and Collision

•
The Lattice BGK (LBGK) equation can be logically and algorithmically separated into two distinct substeps that are performed: collision (or relaxation) and streaming (or propagation).
1.
Collision (Relaxation) Step:
•
This is a purely local operation where the particle populations at each lattice site fi(x,t)f_i(\mathbf{x}, t) relax towards their local equilibrium fieq(x,t)f_i^{\mathrm{eq}}(\mathbf{x}, t).
•
The state of the population after collision fi∗(x,t)f_i^*(\mathbf{x}, t), is calculated as:
fi∗(x,t)=fi(x,t)−Δtτ(fi(x,t)−fieq(x,t)).f_i^*(\mathbf{x}, t) = f_i(\mathbf{x}, t) - \frac{\Delta t}{\tau} \left( f_i(\mathbf{x}, t) - f_i^{\mathrm{eq}}(\mathbf{x}, t) \right).
•
For computational efficiency, this is often implemented as:
fi⋆(x,t)=fi(x,t) (1−Δtτ)+fieq(x,t) Δtτ .f_i^{\star}(\mathbf{x},t) = f_i(\mathbf{x},t)\,\Bigl(1 - \frac{\Delta t}{\tau}\Bigr) + f_i^{\mathrm{eq}}(\mathbf{x},t)\,\frac{\Delta t}{\tau}\,.
◦
The choice τ=Δt\tau = \Delta t is particularly efficient as it simplifies the collision to fi∗=fieqf_i^* = f_i^{\mathrm{eq}}.
2.
Streaming (Propagation) Step:
•
This is a non-local operation; the post-collision populations fi∗(x,t)f_i^*(\mathbf{x}, t) move along their associated velocity vectors ci\mathbf{c}_i to the neighboring lattice sites.
•
This updates the populations for the next time step t+Δtt + \Delta t:
fi(x+ciΔt,t+Δt)=fi∗(x,t).f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f_i^*(\mathbf{x}, t).