Search
moon
sun

[LBM Study] 4. Forces

Category
Study
Kewords
LBM
Kruger
3 more properties
Table of Contents

4.1 Motivation and Background

β€’
Forces play a central role in many hydrodynamic problems.
β€’
Gravitational acceleration g\mathbf{g} can be cast into force density: Fg=ρg\mathbf{F}_g = \rho \mathbf{g}.
β€’
In hydrodynamics, we encounter force densities rather than forces since momentum equation is a PDE for momentum density. β‡’ Forces are momentum density source terms in the Cauchy equation.

Key Applications of Forces in LBM

β€’
Gravity effects:
β—¦
Rayleigh-BΓ©nard instability - convection patterns when warmed fluid rises from hot surface.
β—¦
Rayleigh-Taylor instability - denser fluid descends as lower-density fluid rises.
β—¦
Gravity waves at free water surface.
β€’
Other physical problems:
β—¦
Rotating reference frames subject to radial and Coriolis forces.
β—¦
Charged or magnetic particles in fluids with electromagnetic fields.
β—¦
Electrical Double Layer (EDL) effects in electrolytes near charged surfaces.
β€’
Computational advantages:
β—¦
In incompressible flows, pressure gradients can be replaced by divergence-free body forces.
β—¦
Useful in periodic flow configurations (e.g., porous media flows).
β—¦
Helps avoid accuracy loss from compressibility errors in pressure fields.
β€’
Additional applications:
β—¦
Modeling multiphase or multi-component flows.
β—¦
Fluid-structure interactions via immersed boundary method.

4.2 LBM with Forces in a Nutshell

Order of Operations in a Single Time Step (with BGK collision operator)

β€’
Step 1: Determine force density F\mathbf{F}.
β—¦
Example: gravity force Fg\mathbf{F}_g
β€’
Step 2: Compute fluid density and velocity:
ρ=βˆ‘ifi,u=1Οβˆ‘ifici+FΞ”t2ρ.\rho = \sum_i f_i, \quad\quad \mathbf{u} = \frac{1}{\rho} \sum_i f_i \mathbf{c}_i + \frac{\mathbf{F} \Delta t}{2\rho}.
β€’
Step 3: Compute equilibrium populations fieq(ρ,u)f_i^\mathrm{eq}(\rho, \mathbf{u}) to construct the collision operator:
Ξ©i=βˆ’1Ο„(fiβˆ’fieq).\Omega_i = -\frac{1}{\tau}(f_i - f_i^\mathrm{eq}).
β€’
Step 4: Output macroscopic quantities (optional).
β—¦
Deviatoric stress:
ΟƒΞ±Ξ²β‰ˆβˆ’(1βˆ’Ξ”t2Ο„)βˆ‘ifineqciΞ±ciΞ²βˆ’Ξ”t2(1βˆ’Ξ”t2Ο„)(FΞ±uΞ²+uΞ±FΞ²).\sigma_{\alpha\beta} \approx -\left(1 - \frac{\Delta t}{2\tau}\right) \sum_i f_i^\mathrm{neq} c_{i\alpha} c_{i\beta} - \frac{\Delta t}{2}\left(1 - \frac{\Delta t}{2\tau}\right)(F_\alpha u_\beta + u_\alpha F_\beta).
β€’
Step 5: Compute source term:
Si=(1βˆ’Ξ”t2Ο„)wi(ciΞ±cs2+(ciΞ±ciΞ²βˆ’cs2δαβ)uΞ²cs4)FΞ±.S_i = \left(1 - \frac{\Delta t}{2\tau}\right) w_i \left(\frac{c_{i\alpha}}{c_s^2} + \frac{(c_{i\alpha}c_{i\beta} - c_s^2 \delta_{\alpha\beta})u_\beta}{c_s^4}\right) F_\alpha.
β—¦
Where the source SiS_i and forcing FiF_i terms are related as Si=(1βˆ’12Ο„)FiS_i = (1 - \frac{1}{2\tau})F_i.
β€’
Step 6: Apply collision and source to find the post-collision populations:
fiβˆ—=fi+(Ξ©i+Si)Ξ”t.f_i^* = f_i + (\Omega_i + S_i)\Delta t.
β€’
Step 7: Propagate populations.
β€’
Step 8: Increment the time step and go back to Step 1.

Important Remarks

β€’
The form of the force F\mathbf{F} depends on the underlying physics - not given by LB algorithm.
β€’
Velocity u\mathbf{u} contains the so-called half-force correction.
β—¦
Ensures second-order space-time accuracy.
β—¦
The velocity u\mathbf{u} can be interpreted as the average velocity during the time step, i.e. the average of pre- and post-collision values.
β€’
Forcing scheme presented here is based on a Hermite expansion (same as Guo et al.). β†’ Alternative forcing schemes exist.
β€’
Any cyclic permutation of steps permitted with proper initialization.

4.3 Discretisation

4.3.1 Discretisation in Velocity Space

β€’
Main Question: what is the equivalent polynomial representation in velocity space of the forcing term in the Boltzmann equation?
β€’
The continuous Boltzmann equation with a forcing term:
βˆ‚fβˆ‚t+ΞΎΞ±βˆ‚fβˆ‚xΞ±+FΞ±Οβˆ‚fβˆ‚ΞΎΞ±=Ξ©(f).\frac{\partial f}{\partial t} + \xi_{\alpha} \frac{\partial f}{\partial x_{\alpha}} + \frac{F_{\alpha}}{\rho} \frac{\partial f}{\partial \xi_{\alpha}} = \Omega(f).
β—¦
Goal: Find discrete velocity structure of forcing term aligned with velocity space discretization of feqf^\mathrm{eq}.
β—¦
Challenge: FΞ±\mathbf{F}_\alpha doesn't appear as isolated term but as:
FΞ±Οβˆ‚fβˆ‚ΞΎΞ±\frac{F_{\alpha}}{\rho} \frac{\partial f}{\partial \xi_{\alpha}}

Mathematical Approach

1.
The Hermite series expansion of the distribution function f(ΞΎ)f(\boldsymbol{\xi}) is:
f(x,ΞΎ,t)β‰ˆΟ‰(ΞΎ)βˆ‘n=0N1n!a(n)(x,t)β‹…H(n)(ΞΎ).f(\mathbf{x}, \boldsymbol{\xi}, t) \approx \omega(\boldsymbol{\xi}) \sum_{n=0}^{N} \frac{1}{n!} a^{(n)}(\mathbf{x}, t) \cdot \mathbf{H}^{(n)}(\boldsymbol{\xi}).
2.
The derivative property of Hermite polynomials reads:
Ο‰(ΞΎ)H(n)=(βˆ’1)nβˆ‡ΞΎnΟ‰(ΞΎ).\omega(\boldsymbol{\xi})\mathbf{H}^{(n)} = (-1)^n \nabla_{\boldsymbol{\xi}}^n \omega(\boldsymbol{\xi}).
3.
We can rewrite the Hermite expansion of f(ΞΎi)f(\boldsymbol{\xi}_i) as follows:
fβ‰ˆβˆ‘n=0N(βˆ’1)nn!a(n)β‹…βˆ‡ΞΎnΟ‰.f \approx \sum_{n=0}^{N} \frac{(-1)^n}{n!} a^{(n)} \cdot \nabla_{\boldsymbol{\xi}}^n \omega.
4.
The forcing contribution can be simplified as follows:
FΟβ‹…βˆ‡ΞΎfβ‰ˆFΟβ‹…βˆ‘n=0N(βˆ’1)nn!a(n)β‹…βˆ‡ΞΎn+1Ο‰β‰ˆβˆ’FΟβ‹…Ο‰βˆ‘n=1N1n!na(nβˆ’1)β‹…H(n).\begin{aligned}\frac{\mathbf{F}}{\rho} \cdot \nabla_{\boldsymbol{\xi}} f &\approx \frac{\mathbf{F}}{\rho} \cdot \sum_{n=0}^{N} \frac{(-1)^n}{n!} a^{(n)} \cdot \nabla_{\boldsymbol{\xi}}^{n+1} \omega \\&\approx -\frac{\mathbf{F}}{\rho} \cdot \omega \sum_{n=1}^{N} \frac{1}{n!} n a^{(n-1)} \cdot \mathbf{H}^{(n)}.\end{aligned}

Discrete Form

β€’
Replace continuous ΞΎ\boldsymbol{\xi} with discrete ci\mathbf{c}_i.
β€’
Rescale velocities: ci=ΞΎi/3\mathbf{c}_i = \boldsymbol{\xi}_i/\sqrt{3}.
β€’
Renormalize by lattice weights wiw_i.
β€’
The discrete form of the forcing term:
Fi(x,t)=βˆ’wiΟ‰(ΞΎ)FΟβ‹…βˆ‡ΞΎfβˆ£ΞΎβ†’3ci.F_i(\mathbf{x}, t) = -\frac{w_i}{\omega(\boldsymbol{\xi})} \frac{\mathbf{F}}{\rho} \cdot \nabla_{\boldsymbol{\xi}} f \bigg|_{\boldsymbol{\xi} \rightarrow \sqrt{3} c_i}.
β€’
The discrete velocity Boltzmann equation with a forcing term:
βˆ‚tfi+ciΞ±βˆ‚Ξ±fi=Ξ©i+Fi,i=0,…,qβˆ’1.\partial_t f_i + c_{i\alpha} \partial_\alpha f_i = \Omega_i + F_i, \quad i = 0, \ldots, q - 1.

Second-Order Truncation

β€’
The truncation of the forcing term up to second velocity order (N=2N=2), corresponding to the expansion of feqf^\mathrm{eq}:
Fi=wi(ciΞ±cs2+(ciΞ±ciΞ²βˆ’cs2δαβ)uΞ²cs4)FΞ±.F_i = w_i \left( \frac{c_{i\alpha}}{c_s^2} + \frac{(c_{i\alpha} c_{i\beta} - c_s^2 \delta_{\alpha\beta}) u_\beta}{c_s^4} \right) F_\alpha.
β€’
Velocity moments:
β—¦
Zeroth: βˆ‘iFi=0\sum_i F_i = 0 (no mass source).
β—¦
First: βˆ‘iFiciΞ±=FΞ±\sum_i F_i c_{i\alpha} = F_\alpha (momentum source). β†’ It appears as a body force in the NSE.
β—¦
Second: βˆ‘iFiciΞ±ciΞ²=FΞ±uΞ²+uΞ±FΞ²\sum_i F_i c_{i\alpha} c_{i\beta} = F_\alpha u_\beta + u_\alpha F_\beta (energy source). β†’ a purposefully designed correction to precisely cancel out a spurious error term that arises in the momentum equation of weakly compressible Lattice Boltzmann Methods.
β–ͺ
For a incompressible flow:
Fi=wiciΞ±cs2FΞ±.F_i = w_i \frac{c_{i\alpha}}{c_s^2} F_\alpha.

4.3.2 Discretisation in Space and Time

β€’
Task consists of two parts:
β—¦
Advection: Exact via method of characteristics ( fi=fi(x(ΞΆ),t(ΞΆ))f_i = f_i (\mathbf{x}(\zeta), t(\zeta))):
∫tt+Ξ”tdfidΞΆdΞΆ=fi(x+ciΞ”t,t+Ξ”t)βˆ’fi(x,t).\int_t^{t+\Delta t} \frac{df_i}{d\zeta} d\zeta = f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) - f_i(\mathbf{x}, t).
β—¦
Collision with forces: Requires approximation:
∫tt+Ξ”t(Ξ©i+Fi)dΞΆ.\int_t^{t+\Delta t} (\Omega_i + F_i) d\zeta.

First-Order Integration

β€’
Rectangular discretization:
∫tt+Ξ”t(Ξ©i+Fi)dΞΆ=[Ξ©i(x,t)+Fi(x,t)]Ξ”t+O(Ξ”t2).\int_t^{t+\Delta t} (\Omega_i + F_i) d\zeta = [\Omega_i(\mathbf{x}, t) + F_i(\mathbf{x}, t)] \Delta t + O(\Delta t^2).
β—¦
The integral of collision and forcing terms is approximated by just one point.
β€’
LBE with force (BGK):
fi(x+ciΞ”t,t+Ξ”t)βˆ’fi(x,t)=βˆ’Ξ”tΟ„(fiβˆ’fieq)+FiΞ”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 - f_i^{\mathrm{eq}} \right) + F_i \Delta t.
β€’
Issues:
β—¦
Only first-order accurate in time.
β—¦
Cannot absorb errors into viscosity when forces present.
β—¦
Leads to discrete lattice artifacts.

Second-Order Integration

β€’
Trapezoidal discretization:
∫tt+Ξ”t(Ξ©i+Fi)dΞΆ=(Ξ©i(x,t)+Ξ©i(x+ciΞ”t,t+Ξ”t)2+Fi(x,t)+Fi(x+ciΞ”t,t+Ξ”t)2)Ξ”t+O(Ξ”t3).\int_t^{t+\Delta t} (\Omega_i + F_i) d\zeta = \left( \frac{\Omega_i(\mathbf{x}, t) + \Omega_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t)}{2} + \frac{F_i(\mathbf{x}, t) + F_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t)}{2} \right) \Delta t + O(\Delta t^3).
β—¦
More accurate but time-implicit.
β€’
Change of variables to recover explicit form:
f~i=fiβˆ’(Ξ©i+Fi)Ξ”t2.\tilde{f}_i = f_i - \frac{(\Omega_i + F_i) \Delta t}{2}.
β€’
Second-order accurate LBGK with forcing:
f~i(x+ciΞ”t,t+Ξ”t)βˆ’f~i(x,t)=βˆ’Ξ”tΟ„(f~iβˆ’fieq)+(1βˆ’Ξ”t2Ο„)FiΞ”t\tilde{f}_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) - \tilde{f}_i(\mathbf{x}, t) = -\frac{\Delta t}{\tau} \left( \tilde{f}_i - f_i^{\mathrm{eq}} \right) + \left( 1 - \frac{\Delta t}{2\tau} \right) F_i \Delta t
β—¦
Where Ο„~=Ο„+Ξ”t/2\tilde{\tau} = \tau + \Delta t/2.

Redefined Macroscopic Moments

β€’
Density:
ρ=βˆ‘if~i+Ξ”t2βˆ‘iFi.\rho = \sum_i \tilde{f}_i + \frac{\Delta t}{2} \sum_i F_i.
β€’
Velocity:
ρu=βˆ‘if~ici+Ξ”t2βˆ‘iFiciΞ±.\rho \mathbf{u} = \sum_i \tilde{f}i \mathbf{c}_i + \frac{\Delta t}{2} \sum_i F_i \mathbf{c}_{i \alpha}.
β€’
Momentum flux:
Ξ =(1βˆ’Ξ”t2Ο„)βˆ‘if~icici+Ξ”t2Ο„βˆ‘ifieqcici+Ξ”t2Ο„(1βˆ’Ξ”t2Ο„)βˆ‘iFicici.\boldsymbol{\Pi} = \left( 1 - \frac{\Delta t}{2\tau} \right) \sum_i \tilde{f}_i \mathbf{c}_i \mathbf{c}_i + \frac{\Delta t}{2\tau} \sum_i f_i^{\mathrm{eq}} \mathbf{c}_i \mathbf{c}_i + \frac{\Delta t}{2\tau} \left( 1 - \frac{\Delta t}{2\tau} \right) \sum_i F_i \mathbf{c}_i \mathbf{c}_i.

4.4 Alternative Forcing Schemes

4.4.1 General Observations

β€’
Based on the second-order velocity and space-time discretizations, the LBE with a force can be expressed as:
fi(x+ciΞ”t,t+Ξ”t)βˆ’fi(x,t)=[Ξ©i(x,t)+Si(x,t)]Ξ”t.f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) - f_i(\mathbf{x}, t) = [\Omega_i(\mathbf{x}, t) + S_i(\mathbf{x}, t)] \Delta t.
β—¦
Where Ξ©i\Omega_i is the BGK collision operator and Si=(1βˆ’12Ο„)FiS_i = (1 - \frac{1}{2\tau})F_i denotes a source, with the forcing FiF_i given by:
Fi=wi(ciΞ±cs2+(ciΞ±ciΞ²βˆ’cs2δαβ)uΞ²cs4)FΞ±.F_i = w_i \left( \frac{c_{i\alpha}}{c_s^2} + \frac{(c_{i\alpha} c_{i\beta} - c_s^2 \delta_{\alpha\beta}) u_\beta}{c_s^4} \right) F_\alpha.
β—¦
The fluid velocity in the presence of a force is redefined to guarantee the second-order space-time accuracy:
u=1Οβˆ‘icifi+FΞ”t2ρ.\mathbf{u} = \frac{1}{\rho} \sum_i \mathbf{c}_i f_i + \frac{\mathbf{F} \Delta t}{2\rho}.
β€’
The complexity in the LB literature is caused by the fact that: β‡’ There exist different force algorithms that decompose Ξ©i\Omega_i and SiS_i differently but lead to essentially the same results on the Navier-Stokes level.
β€’
To generalize the forcing method, the equilibrium velocity can be written as :
ueq=1Οβˆ‘ifici+AFΞ”tρ.\mathbf{u}^{eq} = \frac{1}{\rho} \sum_i f_i \mathbf{c}_i + A \frac{\mathbf{F}\Delta t}{\rho}.
β—¦
AA: A model-dependent parameter (A=1/2A = 1/2 for Guo forcing).
β€’
Equivalence condition:
β—¦
Schemes equivalent if Ξ©i+Si\Omega_i + S_i same up to O(F2)O(F^2) or O(u3)O(u^3).
β—¦
Valid only for sufficiently small FF and uu.

4.4.2 Forcing Schemes

Guo et al. (2002)

β€’
Based on the Chapman-Enskog analysis.
β€’
A=1/2A = 1/2.
Si=(1βˆ’Ξ”t2Ο„)wi(ciβˆ’ucs2+(ciβ‹…u)cics4)β‹…FS_i = \left( 1 - \frac{\Delta t}{2\tau} \right) w_i \left( \frac{\mathbf{c}_i - \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i \cdot \mathbf{u}) \mathbf{c}_i}{c_s^4} \right) \cdot \mathbf{F}
β—¦
Removes undesired derivatives in continuity and momentum equations.

Shan and Chen (1993, 1994)

β€’
Originally for multi-phase fluids but applicable to single-phase fluids.
β€’
A=Ο„/Ξ”tA = \tau/\Delta t.
Si=0.S_i = 0.

He et al. (1998)

β€’
Based on near-equilibrium approximation:
Fβ‹…βˆ‡cfβ‰ˆFβ‹…βˆ‡cfeq=βˆ’Fβ‹…cβˆ’ucs2feq.\mathbf{F} \cdot \nabla_{\mathbf{c}} f \approx \mathbf{F} \cdot \nabla_{\mathbf{c}} f^{\mathrm{eq}} = -\mathbf{F} \cdot \frac{\mathbf{c} - \mathbf{u}}{c_s^2} f^{\mathrm{eq}}.
β€’
A=1/2A = 1/2.
Si=(1βˆ’Ξ”t2Ο„)fieqρciβˆ’ucs2β‹…F.S_i = \left( 1 - \frac{\Delta t}{2\tau} \right) \frac{f_i^{\mathrm{eq}}}{\rho} \frac{\mathbf{c}_i - \mathbf{u}}{c_s^2} \cdot \mathbf{F}.

Kupershtokh (2004)

β€’
Exact difference method β†’ To include the force density F\mathbf{F} in such a way that it merely shifts fif_i in velocity space.
β€’
A=0A = 0.
Si=fieq(ρ,uβˆ—+Ξ”u)βˆ’fieq(ρ,uβˆ—).S_i = f_i^{eq}(\rho, u^* + \Delta u) - f_i^{eq}(\rho, u^*).
β—¦
Where uβˆ—=βˆ‘ifici/ρ\mathbf{u}^* = \sum_i f_i \mathbf{c}_i/\rho and Ξ”uu=FΞ”t/ρ\Delta \mathbf{u}u = \mathbf{F} \Delta t/\rho.
β—¦
The equilibrium for a velocity uβˆ—\mathbf{u}^* is directly replaced by the equilibrium for a velocity uβˆ—+Ξ”u\mathbf{u}^* + \Delta \mathbf{u}.

4.5 Chapman-Enskog and Error Analysis in the Presence of Forces

4.5.1 Chapman-Enskog Analysis with Forces

β€’
In order to be consistent with the remaining terms in the LBE, the forcing term must scale as Fi=O(Ο΅)F_i = O(\epsilon). β‡’ We should at least have Fi=Ο΅Fi(1)F_i = \epsilon F_i^{(1)}.
β€’
A hierarchy of Ο΅\epsilon-perturbed equations:
O(Ο΅):(βˆ‚t(1)+ciΞ±βˆ‚Ξ±(1))fieqβˆ’(1βˆ’Ξ”t2Ο„)Fi(1)=1Ο„fi(1),O(Ο΅2):βˆ‚t(2)fieq+(βˆ‚t(1)+ciΞ±βˆ‚Ξ±(1))(1βˆ’Ξ”t2Ο„)(fi(1)+Ξ”t2Fi(1))=βˆ’1Ο„fi(2).\begin{aligned} O(\epsilon): \quad & \left( \partial_t^{(1)} + c_{i\alpha} \partial_\alpha^{(1)} \right) f_i^{\mathrm{eq}} - \left( 1 - \frac{\Delta t}{2\tau} \right) F_i^{(1)} = \frac{1}{\tau} f_i^{(1)}, \\ O(\epsilon^2): \quad & \partial_t^{(2)} f_i^{\mathrm{eq}} + \left( \partial_t^{(1)} + c_{i\alpha} \partial_\alpha^{(1)} \right) \left( 1 - \frac{\Delta t}{2\tau} \right) \left( f_i^{(1)} + \frac{\Delta t}{2} F_i^{(1)} \right) = -\frac{1}{\tau} f_i^{(2)}. \end{aligned}
β€’
In the presence of an external force, the hydrodynamic moments are no longer conserved. β‡’ This leads to a redefinition of the solvability conditions for mass and momentum:
βˆ‘ifineq=βˆ’Ξ”t2βˆ‘iFi(1),βˆ‘icifineq=βˆ’Ξ”t2βˆ‘iciFi(1).\begin{aligned} \sum_i f_i^{\mathrm{neq}} &= -\frac{\Delta t}{2} \sum_i F_i^{(1)}, \\ \sum_i \mathbf{c}_i f_i^{\mathrm{neq}} &= -\frac{\Delta t}{2} \sum_i \mathbf{c}_i F_i^{(1)}. \end{aligned}
β—¦
The extension to β€œstrengthened” order-by-order solvability conditions reads:
βˆ‘ifi(1)=βˆ’Ξ”t2βˆ‘iFi(1)andβˆ‘ifi(k)=0,βˆ‘icifi(1)=βˆ’Ξ”t2βˆ‘iciFi(1)andβˆ‘icifi(k)=0.\begin{aligned} \sum_i f_i^{(1)} &= -\frac{\Delta t}{2} \sum_i F_i^{(1)} \quad \text{and} \quad \sum_i f_i^{(k)} = 0, \\ \sum_i \mathbf{c}_i f_i^{(1)} &= -\frac{\Delta t}{2} \sum_i \mathbf{c}_i F_i^{(1)} \quad \text{and} \quad \sum_i \mathbf{c}_i f_i^{(k)} = \mathbf{0}. \end{aligned}
β–ͺ
With kβ‰₯2k \ge 2, which results from Fi(1)∼O(Ο΅)F_i^{(1)} \sim O(\epsilon).
β€’
By taking the zeroth and first moments of the O(Ο΅)O(\epsilon) equation:
βˆ‚t(1)ρ+βˆ‚Ξ³(1)(ρuΞ³)=0,βˆ‚t(1)(ρuΞ±)+βˆ‚Ξ²(1)Ξ Ξ±Ξ²eq=FΞ±.\begin{aligned}\partial_t^{(1)} \rho + \partial_\gamma^{(1)} (\rho u_\gamma) &= 0, \\\partial_t^{(1)} (\rho u_\alpha) + \partial_\beta^{(1)} \Pi_{\alpha\beta}^{\mathrm{eq}} &= F_\alpha.\end{aligned}
β—¦
Here, Ξ Ξ±Ξ²eq=βˆ‘iciΞ±ciΞ²fieq=ρuΞ±uΞ²+ρcs2δαβ\Pi_{\alpha\beta}^{\mathrm{eq}} = \sum_i c_{i\alpha} c_{i\beta} f_i^{\mathrm{eq}} = \rho u_\alpha u_\beta + \rho c_s^2 \delta_{\alpha\beta}.
β€’
By taking the zeroth and first moments of the O(Ο΅2)O(\epsilon^2) equations :
βˆ‚t(2)ρ=0,βˆ‚t(2)(ρuΞ±)+βˆ‚Ξ²(1)(1βˆ’Ξ”t2Ο„)Ξ Ξ±Ξ²(1)=0.\begin{aligned} \partial_t^{(2)} \rho &= 0, \\ \partial_t^{(2)} (\rho u_\alpha) + \partial_\beta^{(1)} \left( 1 - \frac{\Delta t}{2\tau} \right) \Pi_{\alpha\beta}^{(1)} &= 0. \end{aligned}
β€’
By combining the mass and momentum equations in the O(Ο΅)O(\epsilon) and O(Ο΅)O(\epsilon) equations:
(Ο΅βˆ‚t(1)+Ο΅2βˆ‚t(2))ρ+Ο΅βˆ‚Ξ³(1)(ρuΞ³)=0,(Ο΅βˆ‚t(1)+Ο΅2βˆ‚t(2))(ρuΞ±)+Ο΅βˆ‚Ξ²(1)Ξ Ξ±Ξ²eq=Ο΅FΞ±(1)βˆ’Ο΅2βˆ‚Ξ²(1)(1βˆ’Ξ”t2Ο„)Ξ Ξ±Ξ²(1).\begin{aligned} \left( \epsilon \partial_t^{(1)} + \epsilon^2 \partial_t^{(2)} \right) \rho + \epsilon \partial_\gamma^{(1)} (\rho u_\gamma) &= 0, \\ \left( \epsilon \partial_t^{(1)} + \epsilon^2 \partial_t^{(2)} \right) (\rho u_\alpha) + \epsilon \partial_\beta^{(1)} \Pi_{\alpha\beta}^{\mathrm{eq}} &= \epsilon F_\alpha^{(1)} - \epsilon^2 \partial_\beta^{(1)} \left( 1 - \frac{\Delta t}{2\tau} \right) \Pi_{\alpha\beta}^{(1)}. \end{aligned}
β€’
Ξ Ξ±Ξ²(1)\Pi_{\alpha\beta}^{(1)} is the contribution responsible for the viscous stress at macroscopic level. β‡’ Therefore, the role of βˆ‘iFi(1)ciΞ±ciΞ²\sum_i F_i^{(1)} c_{i\alpha} c_{i\beta} is to remove spurious forcing terms possibly appearing in Ξ Ξ±Ξ²(1)\Pi_{\alpha\beta}^{(1)} so that its form is the same as for the force-free case:
Ξ Ξ±Ξ²(1)=βˆ’Οcs2Ο„(βˆ‚Ξ²(1)uΞ±+βˆ‚Ξ±(1)uΞ²)+O(u3).\Pi_{\alpha\beta}^{(1)} = -\rho c_s^2 \tau \left(\partial_\beta^{(1)} u_\alpha + \partial_\alpha^{(1)} u_\beta\right) + O(u^3).
β—¦
The viscous stress is still given by σαβ=βˆ’(1βˆ’Ξ”t2Ο„)Ξ Ξ±Ξ²(1)\sigma_{\alpha\beta} = -\left( 1 - \frac{\Delta t}{2\tau} \right) \Pi_{\alpha\beta}^{(1)}.
β€’
Finally, we can re-assemble βˆ‚t=Ο΅βˆ‚t(1)+Ο΅2βˆ‚t(2)\partial_t = \epsilon \partial_t^{(1)} + \epsilon^2 \partial_t^{(2)} and use Ξ Ξ±Ξ²eq\Pi_{\alpha\beta}^{\mathrm{eq}} and Ξ Ξ±Ξ²(1)\Pi_{\alpha\beta}^{(1)} to obtain the correct form of the unsteady NSE with forcing term (up to O(u3)O(u^3) error terms):
βˆ‚tρ+βˆ‚Ξ³(ρuΞ³)=0,βˆ‚t(ρuΞ±)+βˆ‚Ξ²(ρuΞ±uΞ²+ρcs2δαβ)=βˆ‚Ξ²[Ξ·(βˆ‚Ξ²uΞ±+βˆ‚Ξ±uΞ²)]+FΞ±.\begin{aligned}\partial_t \rho + \partial_\gamma (\rho u_\gamma) &= 0, \\\partial_t (\rho u_\alpha) + \partial_\beta \left( \rho u_\alpha u_\beta + \rho c_s^2 \delta_{\alpha\beta} \right) &= \partial_\beta \left[ \eta \left( \partial_\beta u_\alpha + \partial_\alpha u_\beta \right) \right] + F_\alpha.\end{aligned}
β—¦
As usual, the dynamic shear and bulk viscosities are Ξ·=ρcs2(Ο„βˆ’Ξ”t2)\eta = \rho c_s^2 \left( \tau - \frac{\Delta t}{2} \right) and Ξ·B=2Ξ·/3\eta_{\mathrm{B}} = 2 \eta / 3, respectively.

4.5.2 Errors Caused by an Incorrect Force Model

Discretization of Velocity Space: The Issue of Unsteady and Steady Cases

β€’
Unsteady state:
β—¦
Term βˆ‚t(1)Ξ Ξ±Ξ²eq\partial_t^{(1)} \Pi_{\alpha\beta}^{eq} contains the contribution FΞ±uΞ²+uΞ±FΞ²F_\alpha u_\beta + u_\alpha F_\beta. β‡’ This contribution can be exactly cancelled by βˆ‘iFiciΞ±ciΞ²\sum_i F_i c_{i\alpha} c_{i\beta}, providing the force term FiF_i is expanded up to the second velocity order.
β€’
Steady state with standard equilibrium:
β—¦
Same spurious term βˆ‘iFiciΞ±ciΞ²\sum_i F_i c_{i\alpha} c_{i\beta} is still required as a correction due to the gradient of the velocity u\mathbf{u}.
β—¦
Still needs second-order expansion.
β€’
Steady state with incompressible equilibrium:
β—¦
The steady incompressible NSE is recovered with no spurious terms. β‡’ We must set βˆ‘iFiciΞ±ciΞ²=0\sum_i F_i c_{i\alpha} c_{i\beta} = 0.

Discretization of Space and Time: The Issue of Discrete Lattice Effects

β€’
Let us assume a time-dependent process and a forcing term with second-order velocity discretization. β‡’ The macroscopic equations reproduced in this case have the following incorrect form:
β—¦
First-order time integration errors:
β–ͺ
Continuity:
βˆ‚tρ+βˆ‚Ξ³(ρuΞ³)=βˆ’Ξ”t2βˆ‚Ξ³FΞ³.\partial_t \rho + \partial_\gamma (\rho u_\gamma) = -\frac{\Delta t}{2} \partial_\gamma F_\gamma.
β–ͺ
Momentum:
βˆ‚t(ρuΞ±)+βˆ‚Ξ²(ρuΞ±uΞ²+ρcs2δαβ)=βˆ‚Ξ²[Ξ·(βˆ‚Ξ²uΞ±+βˆ‚Ξ±uΞ²)]+FΞ±βˆ’Ξ”t2[βˆ‚tFΞ±+βˆ‚Ξ²(uΞ±FΞ²+FΞ±uΞ²)].\partial_t (\rho u_\alpha) + \partial_\beta \left( \rho u_\alpha u_\beta + \rho c_s^2 \delta_{\alpha\beta} \right) = \partial_\beta \left[ \eta \left( \partial_\beta u_\alpha + \partial_\alpha u_\beta \right) \right] + F_\alpha - \frac{\Delta t}{2} \left[ \partial_t F_\alpha + \partial_\beta \left( u_\alpha F_\beta + F_\alpha u_\beta \right) \right].
β€’
Discrete lattice artifacts:
β—¦
Act on same scale as viscous term (O(Ξ”t)O(\Delta t)).
β—¦
Corrupt both mass and momentum equations.
β—¦
More problematic than velocity discretization errors.

4.6 Boundary and Initial Conditions with Forces

4.6.1 Initial Conditions

β€’
Equilibrium initialization with forces:
fi(x,t=0)=fieq(ρ0(x),u~0(x)),u~0=u0βˆ’FΞ”t2ρ0.f_i(\mathbf{x}, t = 0) = f_i^{\mathrm{eq}} \left( \rho_0(\mathbf{x}), \tilde{\mathbf{u}}_0(\mathbf{x}) \right), \quad \tilde{\mathbf{u}}_0 = \mathbf{u}_0 - \frac{\mathbf{F} \Delta t}{2\rho_0}.
β—¦
For low-order forcing schemes, where the macroscopic velocity is computed from ρu=βˆ‘ifici\rho \mathbf{u} = \sum_i f_i \mathbf{c}_i, the equilibrium initialization is the same as in the force-free case, i.e. u~0=u0\tilde{\mathbf{u}}_0 = \mathbf{u}_0.
β€’
Non-equilibrium initialization with forces (adding the modified non-equilibrium term):
fineqβ‰ˆβˆ’wiΟ„cs2ρQiΞ±Ξ²βˆ‚Ξ±uΞ²βˆ’wiΞ”t2cs2(ciΞ±FΞ±+QiΞ±Ξ²2cs2(uΞ±FΞ²+FΞ±uΞ²)).f_i^{\mathrm{neq}} \approx -\frac{w_i \tau}{c_s^2} \rho Q_{i\alpha\beta} \partial_\alpha u_\beta - \frac{w_i \Delta t}{2c_s^2} \left( c_{i\alpha} F_\alpha + \frac{Q_{i\alpha\beta}}{2c_s^2} \left( u_\alpha F_\beta + F_\alpha u_\beta \right) \right).
β—¦
Where QiΞ±Ξ²=ciΞ±ciΞ²βˆ’cs2δαβQ_{i\alpha\beta} = c_{i\alpha}c_{i\beta} - c_s^2 \delta_{\alpha\beta}

4.6.2 Boundary Conditions

Bounce-Back

β€’
The principle of the bounce-back rule is not changed by the inclusion of forces. β‡’ Its accuracy does depend on the force implementation.
β€’
Simple example: a hydrostatic equilibrium where a constant force is balanced by a pressure gradient.
β—¦
The hydrostatic solution established by the bounce-back rule at boundary node xb\mathbf{x}_\mathrm{b}.
β—¦
Second-order space-time discretization for the bulk dynamics (βˆ‚tfi=0\partial_t f_i = 0):
(Ο„βˆ’Ξ”t2)(cs2βˆ‚Ξ±Οβˆ’FΞ±)∣xb=0.\left(\tau - \frac{\Delta t}{2}\right) \left(c_s^2 \partial_\alpha \rho - F_\alpha\right)\bigg|_{\mathbf{x}_\mathrm{b}} = 0.
β–ͺ
The first factor is positive due to the stability requirement Ο„>Ξ”t/2\tau > \Delta t / 2 and can be cancelled. β†’ cs2βˆ‚Ξ±Ο=FΞ±c_s^2 \partial_\alpha \rho = F_\alpha.
β—¦
First-order space-time discretization for the bulk dynamics (βˆ‚tfi=0\partial_t f_i = 0):
(Ο„βˆ’Ξ”t2)(cs2βˆ‚Ξ±Οβˆ’FΞ±)∣xb=Ξ”t2FΞ±(xb).\left(\tau - \frac{\Delta t}{2}\right) \left(c_s^2 \partial_\alpha \rho - F_\alpha\right)\bigg|_{\mathbf{x}_\mathrm{b}} = \frac{\Delta t}{2} F\alpha(\mathbf{x}_\mathrm{b}).
β–ͺ
The first-orderdiscretization retains discrete lattice artifacts even for constant forces.

Non-Equilibrium Bounce-Back (NEBB)

β€’
Modified density calculation at wall:
ρw=cc+uyw(f0+f1+f3+2(f2+f5+f6)+FywΞ”t2c).\rho_{\mathrm{w}} = \frac{c}{c + u_y^{\mathrm{w}}} \left( f_0 + f_1 + f_3 + 2(f_2 + f_5 + f_6) + \frac{F_y^{\mathrm{w}} \Delta t}{2c} \right).
β—¦
The unknown boundary populations still have to be determined by the bounce-back of their non-equilibrium components.
β€’
Momentum corrections:
β—¦
Tangential:
Nx=βˆ’12(f1βˆ’f3)+ρwuxw3cβˆ’FxwΞ”t4c.N_x = -\frac{1}{2}(f_1 - f_3) + \frac{\rho_{\mathrm{w}} u_x^{\mathrm{w}}}{3c} - \frac{F_x^{\mathrm{w}} \Delta t}{4c}.
β—¦
Normal:
Ny=βˆ’FywΞ”t6c.N_y = -\frac{F_y^{\mathrm{w}} \Delta t}{6c}.
β€’
Unknown populations with forces (top wall example):
f4=f2βˆ’2ρwuyw3c+FywΞ”t6c,f7=f5+12(f1βˆ’f3)βˆ’Οwuxw2cβˆ’Οwuyw6c+FxwΞ”t4c+FywΞ”t6c,f8=f6βˆ’12(f1βˆ’f3)+ρwuxw2cβˆ’Οwuyw6cβˆ’FxwΞ”t4c+FywΞ”t6c.\begin{aligned}f_4 &= f_2 - \frac{2\rho_{\mathrm{w}} u_y^{\mathrm{w}}}{3c} + \frac{F_y^{\mathrm{w}} \Delta t}{6c}, \\f_7 &= f_5 + \frac{1}{2}(f_1 - f_3) - \frac{\rho_{\mathrm{w}} u_x^{\mathrm{w}}}{2c} - \frac{\rho_{\mathrm{w}} u_y^{\mathrm{w}}}{6c} + \frac{F_x^{\mathrm{w}} \Delta t}{4c} + \frac{F_y^{\mathrm{w}} \Delta t}{6c}, \\f_8 &= f_6 - \frac{1}{2}(f_1 - f_3) + \frac{\rho_{\mathrm{w}} u_x^{\mathrm{w}}}{2c} - \frac{\rho_{\mathrm{w}} u_y^{\mathrm{w}}}{6c} - \frac{F_x^{\mathrm{w}} \Delta t}{4c} + \frac{F_y^{\mathrm{w}} \Delta t}{6c}.\end{aligned}

4.7 Benchmark Problems

β€’
Velocity Order:
1st:Fi=wiciΞ±cs2FΞ±.\text{1st:} \quad F_i = w_i \frac{c_{i\alpha}}{c_s^2} F_\alpha.
2nd:Fi=wi(ciΞ±cs2+(ciΞ±ciΞ²βˆ’cs2δαβ)uΞ²cs4)FΞ±.\text{2nd:} \quad F_i = w_i \left( \frac{c_{i\alpha}}{c_s^2} + \frac{(c_{i\alpha} c_{i\beta} - c_s^2 \delta_{\alpha\beta}) u_\beta}{c_s^4} \right) F_\alpha.
β€’
Space-Time Order:
1st:fi(x+ciΞ”t,t+Ξ”t)βˆ’fi(x,t)=βˆ’Ξ”tΟ„(fiβˆ’fieq)+FiΞ”t.\text{1st:} \quad f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) - f_i(\mathbf{x}, t) = -\frac{\Delta t}{\tau} \left( f_i - f_i^{\mathrm{eq}} \right) + F_i \Delta t.
2nd:f~i(x+ciΞ”t,t+Ξ”t)βˆ’f~i(x,t)=βˆ’Ξ”tΟ„(f~iβˆ’fieq)+(1βˆ’Ξ”t2Ο„)FiΞ”t\text{2nd:} \quad \tilde{f}_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) - \tilde{f}_i(\mathbf{x}, t) = -\frac{\Delta t}{\tau} \left( \tilde{f}_i - f_i^{\mathrm{eq}} \right) + \left( 1 - \frac{\Delta t}{2\tau} \right) F_i \Delta t
β€’
Comparing four possible forcing strategies:
Scheme
Velocity Order
Space-Time Order
I
1st
1st
II
2nd
1st
III
1st
2nd
IV
2nd
2nd

4.7.1 Problem Description

β€’
A 2D Poiseuille channel flow driven by a combined pressure gradient βˆ‚p/βˆ‚x\partial p/\partial x and body force FxF_x:
ΟΞ½βˆ‚uxβˆ‚y=βˆ‚pβˆ‚xβˆ’Fx.\rho \nu \frac{\partial u_x}{\partial y} = \frac{\partial p}{\partial x} - F_x.
β€’
Analytical velocity solution:
ux(y)=12ρν(βˆ‚pβˆ‚xβˆ’Fx)[y2βˆ’(H2)2]u_x(y) = \frac{1}{2\rho\nu} \left( \frac{\partial p}{\partial x} - F_x \right) \left[ y^2 - \left( \frac{H}{2} \right)^2 \right]
β—¦
Where the no-slip condition (ux=0u_x = 0) holds at the bottom and top walls (y=Β±H/2y = \pm H/2).

4.7.2 Numerical Procedure

β€’
The BGK collision operator with incompressible equilibrium.
β€’
The bounce-back and the non-equilibrium bounce-back (NEBB) boundary conditions are tested.
β€’
Initialization: fi(x,t=0)=fieq(ρ=1,u=0)f_i(\mathbf{x}, t = 0) = f_i^{\mathrm{eq}}(\rho = 1, \mathbf{u} = \mathbf{0}).
β€’
Steady-state criterion: L2≀10βˆ’10L_2 \leq 10^{-10} between 100 consecutive time steps.
β€’
Grid: NxΓ—Ny=5Γ—5N_x \times N_y = 5 \times 5.
MATLAB Code of Poiseuille flow with bounce-back
MATLAB Code of Poiseuille flow with NEBB

4.7.3 Constant Force

β€’
A purely force-driven Poiseuille flow: βˆ‚p/βˆ‚x=0\partial p/\partial x = 0.
β€’
Using periodic boundary conditions at the inlet and outlet.
β€’
The force magnitude is Fx=10βˆ’3F_x = 10^{-3} (in simulation units).

Results:

β€’
Bounce-back: Exact solution at specific Ο„\tau values.
β—¦
Schemes I & II: Ο„=(13/64+5/8)Ξ”t\tau = (\sqrt{13}/64 + 5/8)\Delta t.
β—¦
Schemes III & IV: Ο„=(3/16+1/2)Ξ”t\tau = (\sqrt{3}/16 + 1/2)\Delta t.
β€’
NEBB: Exact for all schemes (no bulk errors for constant force).

4.7.4 Constant Force and Pressure Gradient

β€’
Combined driving: Fxβˆ’βˆ‚p/βˆ‚x=2Γ—10βˆ’3F_x - \partial p/\partial x = 2 \times 10^{-3} (in simulation units) β†’ Considering a 50/50 contribution from each term.
β€’
Using pressure periodic boundary conditions at the inlet and outlet.

Results:

β€’
Bounce-back: Exact solution at specific Ο„\tau values.
β—¦
Schemes I & II: Ο„=Ξ”t\tau = \Delta t.
β—¦
Schemes III & IV: Ο„=(3/16+1/2)Ξ”t\tau = (\sqrt{3}/16 + 1/2)\Delta t.
β–ͺ
Only second-order space-time discretization (Schemes III & IV) maintains physical equivalence.
β€’
NEBB: Exact for all schemes (no bulk errors for constant force).

4.7.5 Linear Force and Pressure Gradient

β€’
The force increases linearly along the streamwise direction.
β—¦
The total contribution remains constant so that the overall magnitude remains locally Fxβˆ’βˆ‚p/βˆ‚x=2Γ—10βˆ’3F_x - \partial p/\partial x = 2 \times 10^{-3} (in simulation units) β†’ Considering a 50/50 contribution from each term.
β€’
The force bulk errors do not vanish any more. β†’ Both bulk and boundary errors can now interfere with the LB solution.

Results:

β€’
Bounce-back: Only Scheme III achieves exact solution at specific Ο„\tau values.
β—¦
Schemes III: Ο„=(3/16+1/2)Ξ”t\tau = (\sqrt{3}/16 + 1/2)\Delta t.
β—¦
First-order velocity discretization required for a steady incompressible flow.
β€’
NEBB: Only Scheme III achieves exact solution for all Ο„\tau values.
β—¦
Schemes III: the bulk solution free from errors.
β—¦
Ο„\tau-independent boundary treatment.

4.7.6 Role of Compressibility

β€’
The standard equilibrium recovers the compressible NSE, which approximates incompressible hydrodynamics in the limit of slow flows and small density (pressure) variations. β‡’ Compressibility Errors.
β€’
The compressibility errors typically have a secondary impact. β†’ They always contaminate the solutions.
β—¦
For spatially varying forces:
β–ͺ
Incompressible equilibrium: First-order velocity discretization optimal.
β–ͺ
Standard equilibrium: Second-order slightly better (but differences small).

Key Takeaways

β€’
Most accurate forcing scheme for general flows: Second-order in both velocity and space-time (Scheme IV).
β€’
For steady incompressible flows: First-order velocity with second-order space-time (Scheme III).
β€’
Critical implementation details:
β—¦
Half-force correction in velocity calculation essential.
β—¦
Boundary conditions must account for forces.
β—¦
Initial conditions need force corrections.
β€’
Error hierarchy: Space-time discretization errors > velocity discretization errors > compressibility errors.