Search
moon
sun

[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).