Search

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