Search

[LBM Study] 3. Boundary and Initial Conditions

Category
Study
Kewords
LBM
Kruger
3 more properties
Table of Contents

3.1 Boundary and Initial Conditions in LBM in a Nutshell

LB Initialization

•
Steady-state solutions: Set initial populations to equilibrium.
◦
fi(x,t=0)=fieq(ρ(x,t=0),u(x,t=0))f_i(\mathbf{x}, t=0) = f_i^\mathrm{eq}(\rho(\mathbf{x}, t=0), u(\mathbf{x}, t=0)).
◦
Typical choice: ρ(x,t=0)=1\rho(\mathbf{x}, t=0) = 1 (simulation units), u(x,t=0)=0u(\mathbf{x}, t=0) = 0.
•
Time-dependent solutions: Require both equilibrium and non-equilibrium components.
◦
fi(x,t=0)=fieq(x,t=0)+fineq(x,t=0)f_i(\mathbf{x}, t=0) = f_i^\mathrm{eq}(\mathbf{x}, t=0) + f_i^\mathrm{neq}(\mathbf{x}, t=0).

Main LB Algorithm Steps

1.
Compute fluid density and velocity via moment equations.
2.
Compute equilibrium populations fieq(ρ,u)f_i^\mathrm{eq}(\rho, \mathbf{u}) from Maxwell-Boltzmann distribution.
3.
Output macroscopic fields if desired.
4.
Apply collision to find post-collision populations fi∗f_i^*.
5.
Propagate populations following streaming step.
6.
Propagate populations at boundaries according to local boundary conditions.
7.
If wet-node boundary conditions used, compute their populations as neccessary.
8.
Increment time step and repeat.

3.1.1 Boundary Conditions

•
Boundary conditions apply at boundary nodes xb\mathbf{x}_b - sites with at least one link to both solid and fluid nodes.
•
LB boundary conditions prescribe mesoscopic populations fif_i rather than macroscopic variables.
•
Two main families of boundary schemes:
◦
Link-wise: Boundary located on lattice links (e.g., bounce-back).
◦
Wet-node: Boundary located on lattice nodes.
•
Key characteristics comparison:

3.1.2 Initial Conditions

•
Two approaches for consistent initial conditions:
◦
Explicit Chapman-Enskog decomposition into equilibrium and non-equilibrium parts.
◦
Modified LB scheme run before actual simulation.

3.2 Fundamentals

3.2.1 Concepts in Continuum Fluid Dynamics

•
Time-dependent problems require initial conditions: u(x,t0)=u0(x)\mathbf{u}(\mathbf{x}, t_0) = \mathbf{u}_0(\mathbf{x}).
•
General form of boundary conditions:
b1∂φ∂n∣(xB,t)+b2φ(xB,t)=b3.b_1\frac{\partial\varphi}{\partial n}\bigg|_{(\mathbf{x}_B,t)} + b_2\varphi(\mathbf{x}_B, t) = b_3.
◦
φ\varphi: the solution of a generic PDE.
◦
xB\mathbf{x}_B: the boundary location.
◦
n\mathbf{n}: the (outward) boundary normal.
•
Boundary conditions classify into three types:
◦
Dirichlet: Fixed value of solution at boundary (b1=0b_1=0 & b2≠0b_2 \neq 0).
◦
Neumann: Fixed flux at boundary (b1≠0b_1 \neq 0 & b2=0b_2 = 0).
◦
Robin: Relation between value and flux at boundary (b1≠0b_1 \neq 0 & b2≠0b_2 \neq 0).
•
For hydrodynamics:
◦
Dirichlet condition for the fluid velocity:
u(xB,t)=UB(xB,t).\mathbf{u}(\mathbf{x}_B, t) = \mathbf{U}_B(\mathbf{x}_B, t).
▪
UB\mathbf{U}_B: the boundary velocity.
◦
Neumann condition for the fluid velocity:
n⋅σ(xB,t)=TB(xB,t).\mathbf{n} \cdot \mathbf{\sigma}(\mathbf{x}_B, t) = \mathbf{T}_B(\mathbf{x}_B, t).
▪
σ\mathbf{\sigma}: the sum of pressure and viscous contributions (total stress tensor) as defined in σαβ=σαβ′−pδαβ\sigma_{\alpha\beta} = \sigma'_{\alpha\beta} - p\delta_{\alpha\beta}.
▪
TB\mathbf{T}_B: the traction vector prescribed at boundary.

3.2.2 Initial Conditions in Discrete Numerical Methods

•
Initial condition still specified by u(x,t0)=u0(x)\mathbf{u}(\mathbf{x}, t_0) = \mathbf{u}_0(\mathbf{x}). ⇒ But a more general step called initialization is also necessary.

3.2.3 Boundary Conditions in Discrete Numerical Methods

•
The order of accuracy of discretized boundary conditions should never be inferior to that of the bulk solution. ⇒ The accuracy of the solution may degrade everywhere; single point with lower accuracy can degrade solution everywhere.

3.2.4 Boundary Conditions for LBM: Introductory Concepts

Which Lattice Sites Need Boundary Conditions?

•
Fluid nodes refer to sites where the LBE applies.
•
Solid nodes are sites completely covered by the solid object where LBE doesn't apply.
•
Boundary nodes link fluid and solid nodes; they require special dynamical rules.
•
Problem: Unknown incoming populations from solid must be specified.
◦
The role of LB boundary conditions is to specify appropriate values for the incoming populations, meaning those moving from the solid object into the fluid region.

Key Differences from Traditional Methods

•
LBM must prescribe mesoscopic populations fif_i rather than macroscopic variables (e.g. velocity or pressure).
•
The fundamental difficulty is that the system of mesoscopic variables has more degrees of freedom than the corresponding macroscopic system. ⇒ Non-uniqueness problem leads to variety of boundary schemes.

How Can Boundary Conditions Be Formulated?

•
Despite the large number of methods available, all LB oundary conditions for straight boundaries belong to one of two groups:
◦
Link-wise family: boundary on lattice links.
◦
Wet-node family: boundary on lattice nodes.

What Defermines the Numerical Accuracy of LB Boundary Conditions?

•
Exactness: Ability to exactly resolve flows of certain order.
•
Accuracy: How error scales with resolution.
•
Chapman-Enskog expansion:
fi=fieq+fi(1)+ϵ2fi(2)+O(ϵ3).f_i = f_i^\mathrm{eq} + f_i^{(1)} + \epsilon^2f_i^{(2)} + O(\epsilon^3).
•
For LB boundaries:
◦
Second-order accurate: fi=fieq+fi(1)+O(Δx2)f_i = f_i^{eq} + f_i^{(1)} + O(\Delta x^2).
◦
Third-order accurate: fi=fieq+fi(1)+ϵ2fi(2)+O(Δx3)f_i = f_i^\mathrm{eq} + f_i^{(1)} + \epsilon^2f_i^{(2)} + O(\Delta x^3).

3.3 Boundary Condition Methods

3.3.1 Periodic Boundary Conditions

•
Periodic boundary conditions apply only to situations where the flow solution is periodic.
◦
They state that the fluid leaving domain on one side will, instantaneously, re-enters at the opposite side.
◦
They conserve mass and momentum at all times.
•
For NSEs, periodic flow conditions along a single dimension are applied as:
ρ(x,t)=ρ(x+L,t),ρu(x,t)=ρu(x+L,t)\begin{aligned} \rho(\mathbf{x}, t) = \rho(\mathbf{x} + \mathbf{L}, t), \\ \rho\mathbf{u}(\mathbf{x}, t) = \rho\mathbf{u}(\mathbf{x} + \mathbf{L}, t) \end{aligned}
◦
L\mathbf{L}: the periodicity direction and length of the flow pattern.
•
For LBM implementation:
fi∗(x,t)=fi∗(x+L,t).f_i^*(\mathbf{x}, t) = f_i^*(\mathbf{x} + \mathbf{L}, t).
◦
For the 2D flow problem:
fi∗(x,t)=fi∗(x+L,t)⟹{f1∗(x0,y2,t)=f1∗(xN,y2,t)f5∗(x0,y2,t)=f5∗(xN,y2,t)f8∗(x0,y2,t)=f8∗(xN,y2,t),fiˉ∗(x+L,t)=fiˉ∗(x,t)⟹{f3∗(xN+1,y2,t)=f3∗(x1,y2,t)f6∗(xN+1,y2,t)=f6∗(x1,y2,t)f7∗(xN+1,y2,t)=f7∗(x1,y2,t).\begin{aligned} f_i^*(\mathbf{x}, t) = f_i^*(\mathbf{x} + \mathbf{L}, t) &\Longrightarrow \begin{cases} f_1^*(x_0, y_2, t) = f_1^*(x_N, y_2, t) \\ f_5^*(x_0, y_2, t) = f_5^*(x_N, y_2, t) \\ f_8^*(x_0, y_2, t) = f_8^*(x_N, y_2, t) \end{cases}, \\ f_{\bar{i}}^*(\mathbf{x} + \mathbf{L}, t) = f_{\bar{i}}^*(\mathbf{x}, t) &\Longrightarrow \begin{cases} f_3^*(x_{N+1}, y_2, t) = f_3^*(x_1, y_2, t) \\ f_6^*(x_{N+1}, y_2, t) = f_6^*(x_1, y_2, t) \\ f_7^*(x_{N+1}, y_2, t) = f_7^*(x_1, y_2, t) \end{cases}. \end{aligned}

3.3.2 Periodic Boundary Conditions with Pressure Variations

•
For incompressible flows, the generalized periodic boundary conditions:
p(x,t)=p(x+L,t)+Δp,u(x,t)=u(x+L,t).\begin{aligned} p(\mathbf{x}, t) = p(\mathbf{x} + \mathbf{L}, t) + \Delta p,\\\mathbf{u}(\mathbf{x}, t) = \mathbf{u}(\mathbf{x} + \mathbf{L}, t). \end{aligned}
•
A simple and robust procedure proposed by Kim and Pitsch (2D flow problem):
1.
Decompose the populations on the virtual nodes (x0x_0, xN+1x_{N+1}) into equilibrium fieqf_i^{\mathrm{eq}} and non-equilibrium fineqf_i^{\mathrm{neq}} parts.
2.
Equilibrium part:
fieq(x0,y,t)=fieq(ρin,uN),fieq(xN+1,y,t)=fieq(ρout,u1).\begin{aligned}f_i^{\mathrm{eq}}(x_0, y, t) &= f_i^{\mathrm{eq}}(\rho_{\mathrm{in}}, \mathbf{u}_N), \\f_i^{\mathrm{eq}}(x_{N+1}, y, t) &= f_i^{\mathrm{eq}}(\rho_{\mathrm{out}}, \mathbf{u}_1). \end{aligned}
•
uN\mathbf{u}_N: the velocity at nodes xNx_N
•
u1\mathbf{u}_1: the velocity at nodes x1x_1
•
The subscripts “in” and “out” denote the pressure values at the left and right physical boundaries, respectively.
•
Δp=pin−pout\Delta p=p_\mathrm{in} - p_\mathrm{out}: the intended pressure drop along the periodicity length L=NΔxL=N \Delta x.
3.
Non-equilibrium part:
fineq(x0,y,t)=fineq(xN,y,t),fineq(xN+1,y,t)=fineq(x1,y,t).\begin{aligned}f_i^{\mathrm{neq}}(x_0, y, t) &= f_i^{\mathrm{neq}}(x_N, y, t), \\f_i^{\mathrm{neq}}(x_{N+1}, y, t) &= f_i^{\mathrm{neq}}(x_1, y, t).\end{aligned}
•
Non-equilibrium populations are determined after collision, i.e. fineq=fi∗−fieq.f_i^{\mathrm{neq}} = f_i^* - f_i^{\mathrm{eq}}.
4.
Finally, we merge the equilibrium part and the non-equilibrium part and perform the streaming step:
fi∗(x0,y,t)=fieq(x0,y,t)+fineq(x0,y,t).f_i^*(x_0, y, t) = f_i^{\mathrm{eq}}(x_0, y, t) + f_i^{\mathrm{neq}}(x_0, y, t).

3.3.3 Solid Boundaries: Bounce-Back Approach

Principle of the Bounce-Back Method

•
Populations hitting rigid wall during propagation are reflected back to where they originally came from.
•
The bounce-back of particles hitting a wall implies no flux across the boundary. → The wall is impermeable to the fluid.
•
The fact that particles are bounced back implies no relative transverse motion between fluid and boundary. → The fluid does not slip on the wall (no-slip boundary).

Implementation Types of Bounce-Back Method

•
Fullway bounce-back (a):
◦
Particles travel complete link path from the boundary to the solid node. ⇒ The particle velocity is inverted in the next collision step.
◦
Requires solid nodes for storage.
◦
2Δt2Δt delay for time-dependent problems.
•
Halfway bounce-back (b):
◦
Particles travel only half of the link distance. ⇒ The particle velocity is inverted during the streaming step.
◦
More accurate for unsteady flows.
•
Both the fullway and halfway approaches assume that the boundary is located approximately midway between solid and boundary nodes, not on the solid nodes themselves.

Resting Walls

•
Physical Mechanism:
◦
Populations leaving the boundary node xb\mathbf{x}_\mathrm{b} at time tt meet the wall surface at t+Δt/2t + \Delta t/2.
◦
Populations are reflected back with a reversed velocity: ciˉ→−ci\mathbf{c}_{\bar i} \rightarrow -\mathbf{c}_i.
◦
Reflected populations return to the boundary node xb\mathbf{x}_\mathrm{b} at t+Δtt + \Delta t.
•
Boundary Condition:
fiˉ(xb,t+Δt)=fi∗(xb,t).f_{\bar{i}}(\mathbf{x}_\mathrm{b}, t + \Delta t) = f_i^*(\mathbf{x}_\mathrm{b}, t).
◦
iˉ\bar{i}: opposite direction to ii.
•
Implementation Example (Bottom Wall):
f2(xb,t+Δt)=f4∗(xb,t),f5(xb,t+Δt)=f7∗(xb,t),f6(xb,t+Δt)=f8∗(xb,t).\begin{aligned} f_2(\mathbf{x}_\mathrm{b}, t + \Delta t) = f_4^*(\mathbf{x}_\mathrm{b}, t), \\ f_5(\mathbf{x}_\mathrm{b}, t + \Delta t) = f_7^*(\mathbf{x}_\mathrm{b}, t), \\ f_6(\mathbf{x}_\mathrm{b}, t + \Delta t) = f_8^*(\mathbf{x}_\mathrm{b}, t). \end{aligned}
•
Key Features:
◦
Simple specular reflection.
◦
No momentum transfer from wall.
◦
Maintains no-slip condition (u=0\mathbf{u} = 0 at wall).

Moving Walls

•
Physical Mechanism:
◦
Moving wall imparts momentum to bouncing particles.
◦
Momentum exchange must satisfy Galilean invariance.
◦
Understood via reference frame transformation approach: (i) transform to the rest frame of the wall, (ii) perform bounce-back there, and (iii) transform back to the initial frame.
•
Boundary Condition (modified bounce-back with momentum correction):
fiˉ(xb,t+Δt)=fi∗(xb,t)−2wiρwci⋅uwcs2.f_{\bar{i}}(\mathbf{x}_\mathrm{b}, t + \Delta t) = f_i^*(\mathbf{x}_\mathrm{b}, t) - 2w_i\rho_\mathrm{w}\frac{\mathbf{c}_i \cdot \mathbf{u}_\mathrm{w}}{c_s^2}.
◦
uw\mathbf{u}_w: prescribed wall velocity.
◦
ρw\rho_w: fluid density at wall location xw=xb+12ciΔt \mathbf{x}_\mathrm{w} = \mathbf{x}_\mathrm{b} + \frac{1}{2}\mathbf{c}_i\Delta t.
◦
Correction term vanishes when uw=0\mathbf{u}_\mathrm{w} = \mathbf{0} (reduces to resting wall).
•
Key Challenges:
◦
Density estimation: ρw\rho_w unknown at the wall (for the general compressible flow). → Approximated by ρ(xb)\rho(\mathbf{x}_\mathrm{b})or ⟨ρ⟩\langle\rho\rangle.
◦
Mass conservation: exact conservation difficult for inclined boundaries.
◦
Accuracy trade-off: mass-conserving strategies may reduce solution accuracy.

Advantages

•
Simple implementation.
•
Stable even near τ→Δt/2\tau \rightarrow \Delta t/2.
•
Exact mass conservation at resting boundaries.
•
It can be implemented straightforwardly for any number of spatial dimensions.

Disadvantages

•
The bounce-back rule can only approximates arbitrary surfaces through "staircase" shapes.
•
The exact location of the no-slip boundary is viscosity-dependent when the bounce-back scheme is used together with the BGK collision model. ⇒ This defect may be solved by replacing the BGK model by a more complex collision model such as TRT or MRT.

Numerical Evaluation of the Accuracy of the Bounce-Back Method

•
Couette flow: first-order analysis.
◦
For Couette flow we observed that: (i) bounce-back places the no-slip wall exactly midway between boundary and solid nodes, (ii) solutions are viscosity-independent. ⇒ Bounce-back, with the wall shifted by Δx/2\Delta x/2 from boundary nodes, correctly prescribes the boundary value of the velocity and its first-order derivative.
MATLAB Code of Couette flow
•
Poiseuille flow: second-order analysis.
◦
Unlike the numerical Couette flow solution, the Poiseuille flow simulation is affected by the choice of τ\tau. → Numerical Boundary Slip
◦
For Poiseuille flow we observed that: the bounce-back rule performs better with the walls located midway between boundary and solid nodes, just as for Couette flow. ⇒ However, the parabolic solution is generally not exact, due to an error depending on τ2\tau^2 and therefore on viscosity. ⇒ The slip error comes from the O(Δx2)O(\Delta x^2) term in the bounce-back closure relation.
◦
To avoid the unphysical viscosity-dependent slip, the safest procedure is to replace BGK by the TRT or MRT collision operators.
MATLAB Code of Poiseuille flow

3.3.4 Solid Boundaries: Wet-Node Approach

Finding the Density on Boundaries

•
Use continuity condition as constraint.
◦
For top wall:
ρw=cc+uw,y[f0+f1+f3+2(f2+f5+f6)].\rho_\mathrm{w} = \frac{c}{c + u_{\mathrm{w},y}}[f_0 + f_1 + f_3 + 2(f_2 + f_5 + f_6)].
◦
For bottom wall:
ρw=cc−uw,y[f0+f1+f3+2(f4+f7+f8)].\rho_\mathrm{w} = \frac{c}{c - u_{\mathrm{w},y}}[f_0 + f_1 + f_3 + 2(f_4 + f_7 + f_8)].

Equilibrium Scheme

•
Simplest approach:
fi(xb,t)=fieq(ρw,uw).f_i(\mathbf{x}_\mathrm{b}, t) = f_i^\mathrm{eq}(\rho_\mathrm{w}, \mathbf{u}_\mathrm{w}).
◦
Replaces all populations fif_i at boundary, instead of only the unknown ones.
•
Advantages:
◦
(i) Excellent stability, (ii) Easy to apply in 2D/3D.
•
Disadvantages:
◦
(i) Only accurate when τ=Δt\tau = \Delta t, (ii) First-order accuracy otherwise.
•
The linear incompressible equilibrium populations on boundary nodes:
fieq(xb,t)=wi(pwcs2+ρ0ci⋅uwcs2).f_i^{\mathrm{eq}}(\mathbf{x}_\mathrm{b}, t) = w_i \left( \frac{p_{\mathrm{w}}}{\mathbf{c}_s^2} + \rho_0 \frac{\mathbf{c}_i \cdot \mathbf{u}_{\mathrm{w}}}{\mathbf{c}_s^2} \right).
◦
For top wall:
pwcs2=ρ0uw,yc+[f0+f1+f3+2(f2+f5+f6)].\frac{p_{\mathrm{w}}}{\mathbf{c}_s^2} = \frac{\rho_0 u_{\mathrm{w},y}}{c} + [f_0 + f_1 + f_3 + 2(f_2 + f_5 + f_6)].
◦
For bottom wall:
pwcs2=ρ0uw,yc+[f0+f1+f3+2(f4+f7+f8)].\frac{p_{\mathrm{w}}}{\mathbf{c}_s^2} = \frac{\rho_0 u_{\mathrm{w},y}}{c} + [f_0 + f_1 + f_3 + 2(f_4 + f_7 + f_8)].

Non-Equilibrium Extrapolation Method

•
Includes non-equilibrium part from neighboring fluid node:
fi(xb,t)=fieq(ρw,uw)+(fi(xf,t)−fieq(ρf,uf)).f_i(\mathbf{x}_\mathrm{b}, t) = f_i^{\mathrm{eq}}(\rho_{\mathrm{w}}, \mathbf{u}_{\mathrm{w}}) + \left( f_i(\mathbf{x}_f, t) - f_i^{\mathrm{eq}}(\rho_f, \mathbf{u}_f) \right).
◦
Replaces all populations fif_i at boundary, instead of only the unknown ones.
◦
Second-order accurate. → Cannot capture parabolic solutions exactly
•
Non-equilibrium approximation:
fineq(x,t)≃−wiτρcs2Qiαβ∂βuα.f_i^{\mathrm{neq}}(\mathbf{x}, t) \simeq -w_i \frac{\tau \rho}{c_s^2} Q_{i\alpha\beta} \partial_\beta u_\alpha.
◦
where Qiαβ​=ciα​ciβ​−cs2​δαβ​.Q_{iαβ}​=c_{iα}​c_{iβ} ​− c_{s}^2 ​δ_{αβ}​.
•
Couette flow: NEEM provides exact solution regardless of τ\tau value.
•
Poiseuille flow: Method introduces errors by implicitly enforcing ∇u∣(x,y1)=∇u∣(x,y2)∇\mathbf{u}|_{(x,y_1)} = ∇\mathbf{u}|_{(x,y_2)}, which cannot be satisfied in parabolic solutions.

Non-Equilibrium Bounce-Back Method (Zou-He Method)

•
The non-equilibrium bounce-back (NEBB) method enforces bounce-back for the non-equilibrium part:
fiˉneq(xb,t)=fineq(xb,t)(ciˉ=−ci).f_{\bar{i}}^{\mathrm{neq}}(\mathbf{x}_\mathrm{b}, t) = f_i^{\mathrm{neq}}(\mathbf{x}_\mathrm{b}, t) \quad (\mathbf{c}_{\bar{i}} = -\mathbf{c}_i).
◦
There will be no way to guarantee that the tangential velocity obtains the intended value. ⇒ We used two macroscopic conditions: (i) density (or pressure) and (ii) normal velocity.
◦
The NEBB is modified with transverse momentum correction NtN_t:
fiˉneq(xb,t)=fineq(xb,t)−t⋅ci∣ci∣Nt.f_{\bar{i}}^{\mathrm{neq}}(\mathbf{x}_\mathrm{b}, t) = f_i^{\mathrm{neq}}(\mathbf{x}_\mathrm{b}, t) - \frac{\mathbf{t} \cdot \mathbf{c}_i}{|\mathbf{c}_i|} N_t.
▪
t\mathbf{t}: a tangent unit vector along the wall.
◦
Third-order accurate, viscosity-independent → Exact for parabolic solutions.
◦
The NEBB method modifies only the missing populations.

3.3.5 Open Boundaries

Open boundaries consist of inlets or outlets where the flow either enters or leaves the computational domain and where we should impose, for example, velocity or density profiles.
This is often non-trivial, and it can cause both physical and numerical difficulties.

Velocity Boundary Conditions

•
Bounce-back approach:
fiˉ(xb,t+Δt)=fi∗(xb,t)−2wiρwci⋅uwcs2.f_{\bar{i}}(\mathbf{x}_\mathrm{b}, t + \Delta t) = f_i^*(\mathbf{x}_\mathrm{b}, t) - 2w_i\rho_\mathrm{w}\frac{\mathbf{c}_i \cdot \mathbf{u}_\mathrm{w}}{c_s^2}.
◦
Same as the moving wall bounce-back with prescribed uw\mathbf{u}_\mathrm{w}. ⇒ The bounce-back algorithm still places the inlet/outlet boundary, with its imposed velocity uw\mathbf{u}_\mathrm{w}, midway between the lattice nodes.
•
Wet-node approach:
◦
The wet-node formulation of velocity boundary conditions applies similarly to open boundaries and walls. ⇒ The only difference is that the velocity represents the mass flux entering or leaving the domain, rather than the movement of a solid wall.

Pressure Boundary Conditions

•
Anti-bounce-back approach:
fiˉ(xb,t+Δt)=−fi∗(xb,t)+2wiρw[1+(ci⋅uw)22cs4−uw22cs2].f_{\bar{i}}(\mathbf{x}_\mathrm{b}, t + \Delta t) = -f_i^*(\mathbf{x}_\mathrm{b}, t) + 2w_i \rho_{\mathrm{w}} \left[ 1 + \frac{(\mathbf{c}_i \cdot \mathbf{u}_{\mathrm{w}})^2}{2\mathbf{c}_s^4} - \frac{\mathbf{u}_{\mathrm{w}}^2}{2\mathbf{c}_s^2} \right].
◦
The boundary is located Δx/2\Delta x / 2 outside the boundary node.
◦
It requires specification of uw\mathbf{u}_\mathrm{w} which is not generally known. → A possible way to estimate uw\mathbf{u}_\mathrm{w} is by extrapolation.
◦
Formally third-order accurate.
•
Wet-node approach:
◦
In the wet-node approach, the algorithm for imposing pressure and velocity boundary conditions is the same. ⇒ The only difference lies in the strategy to find the unknown macroscopic properties at the boundary.
◦
The impermeable wall condition is solved for boundary velocity:
uw,y=−c+cρw[f0+f1+f3+2(f2+f5+f6)].u_{\mathrm{w},y} = -c + \frac{c}{\rho_{\mathrm{w}}} \left[ f_0 + f_1 + f_3 + 2(f_2 + f_5 + f_6) \right].

3.3.6 Corners

In 2D domains, we have to deal with geometrical features where two straight surfaces intersect, called corners.
In 3D domains, we need to consider both edges and corners, i.e. places where two and three surfaces cross, respectively.

Corners and the Bounce-Back Rule: (a)

•
The bounce-back rule works in the same way for straight walls and corners, and it is the same for concave and convex cases ⇒ The unknown populations leaving corners are determined through the full reflection of the known incoming ones.
•
Concave corner: f1f_1, f2f_2, and f5f_5 unknown
•
Convex corner: only f5f_5 unknown
•
The advantages and disadvantages of the bounce-back rule for corners are similar to those at planar surfaces. ⇒ The ease of implementation, strict conservation of mass, and good stability characteristics (even close to τ=Δt/2\tau = \Delta t /2). ⇒ The lower accuracy, also featuring viscosity-dependent errors (if using the BGK model).

Corners and the NEBB Method: (b)

•
The main point is that we need to ensure that both known and unknown populations produce the desired macroscopic properties at corners, just as they do at straight boundaries. ⇒ Finding the corner density (pressure) generally requires extrapolation from neighbouring nodes. ⇒ All macroscopic quantities prescribed at corners (ρw\rho_\mathrm{w}, uw,xu_{\mathrm{w}, x}, uw,yu_{\mathrm{w}, y}) are known.
•
Concave corners (six unknowns: f0f_0, f1f_1, f2f_2, f5f_5, f6f_6, f8f_8):
◦
Apply NEBB rule:
f1=f3+23cρwuw,x,f2=f4+23cρwuw,y,f5=f7+16cρw(uw,x+uw,y).\begin{aligned}f_1 &= f_3 + \frac{2}{3c} \rho_{\mathrm{w}} u_{\mathrm{w},x}, \\f_2 &= f_4 + \frac{2}{3c} \rho_{\mathrm{w}} u_{\mathrm{w},y}, \\f_5 &= f_7 + \frac{1}{6c} \rho_{\mathrm{w}} (u_{\mathrm{w},x} + u_{\mathrm{w},y}).\end{aligned}
◦
Use conservation laws:
f6=112cρw(uw,y−uw,x),f8=112cρw(uw,x−uw,y),f0=ρw−∑i=18fi.\begin{aligned}f_6 &= \frac{1}{12c} \rho_{\mathrm{w}} (u_{\mathrm{w},y} - u_{\mathrm{w},x}), \\f_8 &= \frac{1}{12c} \rho_{\mathrm{w}} (u_{\mathrm{w},x} - u_{\mathrm{w},y}), \\f_0 &= \rho_{\mathrm{w}} - \sum_{i=1}^{8} f_i.\end{aligned}
•
Convex corners (one unknown: f5f_5):
◦
Over-specified with full NEBB approach. ⇒ Often use simple (node) bounce-back rule:
f5=f7.f_5 = f_7.

3.3.7 Symmetry and Free-Slip Boundaries

Mirror Symmetry

•
For a symmetry plane lying exactly midway on the links between two rows of nodes, the populations on one side of the boundary must exactly mirror those on the other.
◦
For a symmetry plane lying at y: fi(x,y−Δx/2,t)=fj(x,y+Δx/2,t)f_i(x, y - \Delta x/2, t) = f_j(x, y + \Delta x/2, t) where cj,n=−ci,nc_{j,n} = -c_{i,n} (normal velocity reversed).
•
Implementation:
◦
Treat symmetry boundary as bounce-back-like mechanism where populations leaving boundary node xb\mathbf{x}_b at time tt meet symmetry surface at time t+Δt2t + \frac{\Delta t}{2}.
◦
Populations undergo specular reflection with normal velocity component reversed from the incoming velocity ci\mathbf{c}_i: cj,n=−ci,nc_{j,n} = -c_{i,n} while maintaining same resulting velocity cj\mathbf{c}_j.
◦
Reflected populations arrive at time t+Δtt + \Delta t at node xb\mathbf{x}_b or neighboring boundary nodes.
◦
Standard streaming step replaced by:
fj(xb+cj,tΔt,t+Δt)=fi∗(xb,t).f_j(\mathbf{x}_\mathrm{b} + \mathbf{c}_{j,t} \Delta t, t + \Delta t) = f_i^*(\mathbf{x}_\mathrm{b}, t).
▪
cj,t=ci,t\mathbf{c}_{j,t}=\mathbf{c}_{i,t} represents tangential velocity component of populations, equating ci\mathbf{c}_i and cj\mathbf{c}_j with their normal velocity set to zero.

Free-Slip Boundary

•
A free-slip boundary condition enforces a zero normal fluid velocity un=0\mathbf{u}_n=0 but places no restrictions on the tangential fluid velocity ut\mathbf{u}_t.
•
It can be implemented in exactly the same way as the symmetry boundary condition.

3.4 Further Topics on Boundary Conditions

3.4.1 The Chapman-Enskog Analysis for Boundary Conditions

•
LB boundary conditions rely on specific closure rules that differ from bulk collision and propagation rules.
•
Chapman-Enskog analysis can determine macroscopic behavior of LB boundary conditions.

Bounce-Back Method

•
Macroscopic Dirichlet velocity condition applies at lattice links rather than grid nodes. ⇒ The understanding of what happens at links must be examined in terms of continuation behaviour.
•
For the third-order boundary accuracy, use the second-order Taylor expansion:
ux∣yw=ux∣yb+Δx2∂yux∣yb+12(Δx2)2∂y∂yux∣yb+O(Δx3).u_x|_{y_\mathrm{w}} = u_x|_{y_\mathrm{b}} + \frac{\Delta x}{2}\partial_y u_x|_{y_\mathrm{b}} + \frac{1}{2}\left(\frac{\Delta x}{2}\right)^2\partial_y\partial_y u_x|_{y_\mathrm{b}} + O(\Delta x^3).
◦
For a link-wise approach, to impose a wall velocity with formal third-order accuracy, its closure relation must satisfy above equation with yw=yb+Δx/2.y_\mathrm{w} = y_\mathrm{b} + \Delta x/2.
•
If the solution is time-independent, the second-order Chapman-Enskog analysis starts from linearized equilibrium:
fieq=wi(ρ+ρ0ciγuγcs2),fi(1)=−τwiciα∂α(1)(ρ+ρ0ciγuγcs2),fi(2)=τ(τ−Δt2)wiciαciβ∂α(1)∂β(1)(ρ+ρ0ciγuγcs2).\begin{aligned}f_i^{\mathrm{eq}} &= w_i \left( \rho + \rho_0 \frac{c_{i\gamma} u_\gamma}{c_s^2} \right), \\f_i^{(1)} &= -\tau w_i c_{i\alpha} \partial_\alpha^{(1)} \left( \rho + \rho_0 \frac{c_{i\gamma} u_\gamma}{c_s^2} \right), \\f_i^{(2)} &= \tau \left( \tau - \frac{\Delta t}{2} \right) w_i c_{i\alpha} c_{i\beta} \partial_\alpha^{(1)} \partial_\beta^{(1)} \left( \rho + \rho_0 \frac{c_{i\gamma} u_\gamma}{c_s^2} \right).\end{aligned}
•
The second-order Chapman-Enskog expansion of the bounce-back formula, fiˉ(xb,t+Δt)=fi∗(xb,t)−2wiρwci⋅uwcs2f_{\bar{i}}(\mathbf{x}_\mathrm{b}, t + \Delta t) = f_i^*(\mathbf{x}_\mathrm{b}, t) - 2w_i\rho_\mathrm{w}\frac{\mathbf{c}_i \cdot \mathbf{u}_\mathrm{w}}{c_s^2}, proceeds with the decomposition of its post-streaming and post-collision populations as follows (with uw,γ=uγ∣wu_{\mathrm{w}, \gamma} = u_\gamma|_\mathrm{w}):
fiˉeq+εfiˉ(1)+ε2fiˉ(2)=fieq+(1−Δtτ)(εfi(1)+ε2fi(2))−2wiρ0ciγuw,γcs2.f_{\bar{i}}^{\mathrm{eq}} + \varepsilon f_{\bar{i}}^{(1)} + \varepsilon^2 f_{\bar{i}}^{(2)} = f_i^{\mathrm{eq}} + \left( 1 - \frac{\Delta t}{\tau} \right) \left( \varepsilon f_i^{(1)} + \varepsilon^2 f_i^{(2)} \right) - 2w_i \rho_0 \frac{c_{i\gamma} u_{\mathrm{w},\gamma}}{c_s^2}.
•
The steady-state macroscopic content of the bounce-back formula at node yby_\mathrm{b} is obtained by performing some algebraic manipulations (using p=cs2ρp = c_s^2 \rho and ∂α=ϵ∂α(1)\partial_\alpha = \epsilon \partial_\alpha^{(1)}):
•
Assuming that this problem applies for a unidirectional horizontal flow aligned with the lattice so that ciγuγ=cixuxc_{i \gamma} u_\gamma = c_{i x} u_x and ciα∂α=ciy∂xc_{i \alpha} \partial_\alpha = c_{i y} \partial_x. → Only diagonal lattice links play a role. ⇒ The pressure terms have to be re-expressed in the form of velocity corrections (assuming a unidirectional pressure-driven flow): ∂x(p/ρ0)=cs2(τ−Δt/2)∂y2ux\partial_x (p/\rho_0) = c_s^2 (\tau - \Delta t/2) \partial_y^2 u_x and ∂x2(p/ρ0)=0.\partial_x^2 (p/\rho_0) = 0.
•
The steady closure relation for the bounce-back rule on a top wall can be re-expressed as:
uw,x=ux∣yb+Δx2∂yux∣yb+2c23(τ−Δt2)2⏟=Δx2/8∂y2ux∣yb.u_{\mathrm{w},x} = u_x \big|_{y_{\mathrm{b}}} + \frac{\Delta x}{2} \partial_y u_x \bigg|_{y_{\mathrm{b}}} + \underbrace{\frac{2c^2}{3} \left( \tau - \frac{\Delta t}{2} \right)^2}_{=\Delta x^2/8} \partial_y^2 u_x \bigg|_{y_{\mathrm{b}}} .
◦
Viscosity-independent boundary location requires (τ−Δt/2)2=(3/16)Δt2(\tau - \Delta t/2)^2 = (3/16)\Delta t^2 (remember that cs2=c2/3=Δx2/(3Δt2)c_s^2 = c^2/3 = \Delta x^2/(3 \Delta t^2)).
◦
Other τ\tau values lead to viscosity-dependent boundary shift or velocity slip.

Non-Equilibrium Bounce-Back Method

•
Wet boundary nodes operate with same dynamics as bulk nodes.
•
The NEBB method correctly reproduces macroscopic equation at wall:
∂xpw−c23(τ−Δt2)ρ0∂y2uw,x=0.\partial_x p_\mathrm{w} - \frac{c^2}{3}\left(\tau - \frac{\Delta t}{2}\right)\rho_0\partial_y^2 u_{\mathrm{w},x} = 0.
•
The NEBB method exactly describes the pressure-driven (Poiseuille) channel flow on the boundary, thereby agreeing with the LBE in the bulk. ⇒ The NEBB is a third-order accurate boundary scheme for straight walls coinciding with lattice nodes.

3.4.2 Mass Conservation at Solid Boundaries

•
Mass conservation is crucial for compressible, multi-phase, or multi-component flows.
•
The LB algorithm ensures that mass is conserved locally since collision leaves ∑ifi=ρ\sum_i f_i = \rho invariant. ⇒ However, there is no a priori guarantee that this property holds at boundary nodes, which necessarily must behave differently than bulk fluid nodes.

Link-Wise Approach: Bounce-Back

•
Bounce-back places wall Δx/2\Delta x/2 outside boundary node. → For straight walls, bounce-back is always mass conserving.
•
Mass is not conserved when:
◦
Normal mass flux is prescribed (e.g., normal velocity in Dirichlet correction).
◦
Boundary motions create or destroy fluid nodes to solid nodes and vice-versa.
◦
Boundary shapes don't align with lattice nodes. → Numerical errors due to the wall discretization contaminate the local mass balance.

Wet-Node Approach: Non-Equilibrium Bounce-Back

•
Boundary node considered in fluid domain, infinitesimally close to solid. → Reconstruction doesn't automatically conserve local mass.
•
For a top wall, the net mass Δm\Delta m on boundary node is:
Δm=min−mout=(f2+f5+f6)−(f4∗+f7∗+f8∗).\Delta m = m_\mathrm{in} - m_\mathrm{out} = (f_2 + f_5 + f_6) - (f_4^* + f_7^* + f_8^*).
•
Mass conservation up to O(ϵ3)O(\epsilon^3) for the NEBB:
Δm=ρwuw,y−Δt2∂y(ρwuw,y)+Δt6∂yρw.\Delta m = \rho_\mathrm{w} u_{\mathrm{w},y} - \frac{\Delta t}{2}\partial_y(\rho_\mathrm{w} u_{\mathrm{w},y}) + \frac{\Delta t}{6}\partial_y\rho_\mathrm{w}.
•
Exact mass conservation methods exist but may degrade overall solution accuracy

3.4.3 Momentum Exchange at Solid Boundaries

•
The total force f\mathbf{f} acting on a boundary area AA can be written as the surface integral of the traction vector t=σw⋅n^\mathbf{t} = \mathbf{\sigma_\mathrm{w}} \cdot \hat{\mathbf{n}}:
fα=∫dAα=∫dA σw,αβn^β.f_\alpha = \int \mathrm{d}A_\alpha = \int \mathrm{d}A \, \sigma_{\mathrm{w},\alpha\beta} \hat{n}_\beta.
◦
σw\mathbf{\sigma}_\mathrm{w}: the stress tensor at the surface.
◦
n^\hat{\mathbf{n}}: the unit normal vector pointing from the surface into the fluid.
•
LBM allows direct access to the stress tensor without complicated computations. ⇒ In the bulk, the viscous shear stress follows directly from the second-order velocity moment of fineqf_i^{\mathrm{neq}}.

Momentum Exchange in the Bounce-Back Method

•
Momentum Exchange Algorithm (MEA) by Ladd:
◦
Identify those populations that cross the boundary (both from the fluid into the solid and the other way around).
◦
Sum up all corresponding momentum contributions to obtain the momentum exchange at the wall.
◦
Note that the MEA provides only the traction vector t\mathbf{t} rather than the full stress tensor σw\mathbf{\sigma_\mathrm{w}}.
•
For each boundary link:
◦
Momentum exchange: Δp(xiw)=piin−piout\Delta p(x_i^\mathrm{w}) = p_i^{in} - p_i^{out}.
•
Total momentum exchange is the sum of all those contributions by all bounced back populations:
1.
Identify all links between boundary and solid nodes with locations xiw\mathbf{x}_i^\mathrm{w}.
2.
At each time step, evaluate the incoming and bounced back populations fiinf_i^\mathrm{in} and fiˉoutf_{\bar{i}}^\mathrm{out} at each of the identified links.
3.
The total momentum exchange during one streaming step is the sum over all identified links and, due to ciˉ=ci\mathbf{c}_{\bar{i}} = \mathbf{c}_i, can be written as:
ΔP=Δx3∑xiwΔp(xiw)=Δx3∑xiw(fiin+fiout)ci.\Delta \mathbf{P} = \Delta x^3 \sum_{\mathbf{x}_i^{\mathrm{w}}} \Delta \mathbf{p}(\mathbf{x}_i^{\mathrm{w}}) = \Delta x^3 \sum_{\mathbf{x}_i^{\mathrm{w}}} \left( f_i^{\mathrm{in}} + f_i^{\mathrm{out}} \right) \mathbf{c}_i.
◦
The prefactor Δx3\Delta x^3in 3D (or equivalently Δx2\Delta x^2 in 2D) ensures that the result is a momentum rather than a momentum density.
•
If we now assume that the momentum is exchanged smoothly during one time step Δt\Delta t, we can easily find the force acting on the boundary:
f=ΔPΔt.\mathbf{f} = \frac{\Delta \mathbf{P}}{\Delta t}.
◦
Similarly, the angular momentum exchange ΔL\Delta \mathbf{L} and therefore the total torque T=ΔL/Δt\mathbf{T} = \Delta \mathbf{L}/\Delta t acting on the wall can be computed. ⇒ Replace Δp(xiw)\Delta \mathbf{p}(\mathbf{x}_i^{\mathrm{w}}) by (xiw−xref)×Δp(xiw)(\mathbf{x}_i^\mathrm{w} - \mathbf{x}_\mathrm{ref}) \times \Delta \mathbf{p}(\mathbf{x}_i^{\mathrm{w}}) where xref\mathbf{x}_\mathrm{ref} is a fixed reference point, e.g. the origin of the coordinate system.

Accuracy of the Momentum Exchange Method

•
The MEA provides an accurate measure of wall shear stress. → It can be verified via a Chapman-Enskog analysis.
•
The net tangential momentum transferred from fluid to solid across the wall, which is located midway between the solid node xs\mathbf{x}_\mathrm{s} and the boundary node xb\mathbf{x}_\mathrm{b}:
Δpxw=cVΔtA[(f5−f6)∣xb−(f8−f7)∣xs]=c2[(f5−f6)−(f8∗−f7∗)]∣xb.\Delta p_x^\mathrm{w} = \frac{cV}{\Delta tA}[(f_5 - f_6)\big|_{x_\mathrm{b}} - (f_8 - f_7)\big|_{x_\mathrm{s}}] = c^2[(f_5 - f_6) - (f_8^* - f_7^*)]\big|_{x_\mathrm{b}}.
◦
In 3D: V=Δx3V = \Delta x^3 and A=Δx2A = \Delta x^2. | In 2D: V=Δx2V = \Delta x^2 and A=ΔxA = \Delta x.
•
First-order Taylor expansion of viscous stress σxy=−νρ0∂yux\sigma_{xy} = - \nu \rho_0 \partial_y u_x:
Δpxw=−cs2(τ−Δt2)ρ0∂yux∣xb+Δx2cs2(τ−Δt2)ρ0∂y2ux∣xb.\Delta p_x^{\mathrm{w}} = -c_s^2 \left( \tau - \frac{\Delta t}{2} \right) \rho_0 \partial_y u_x \bigg|_{\mathbf{x}_{\mathrm{b}}} + \frac{\Delta x}{2} c_s^2 \left( \tau - \frac{\Delta t}{2} \right) \rho_0 \partial_y^2 u_x \bigg|_{\mathbf{x}_{\mathrm{b}}} .
◦
The MEA is second-order accurate and exact for parabolic solutions.

3.4.4 Boundary Conditions in 3D

•
Geometry differences: 2D has edges and corners; 3D has planes, edges, and corners.
•
Number of unknown populations for different lattices and boundary configurations.
•
Link-wise boundary methods specify the missing populations based on simple reflection rules, e.g. bounce-back. → The link-wise methods naturally extend to 3D.
•
Wet-node methods face challenges:
◦
Increased number of unknowns.
◦
Buried links in concave edges and corners.
◦
Over-constrained problems at convex nodes.
•
Wet-node boundary conditions may yield superior accuracy over bounce-back on flat surfaces.
◦
3D implementation challenges: Difficult due to need for distinguishing populations by wall orientation and cumbersome treatment of edge and corner nodes.
•
Link-wise boundary schemes do not distinguish between populations and are easy to implement (particularly bounce-back).
◦
General 3D recommendation: Link-wise boundary methods offer better compromise between accuracy and implementation ease compared to wet-node techniques.

3.5 Initial Conditions

3.5.1 Steady and Unsteady Situations

•
LB simulations fall into one of the following categories:
1.
Steady flows: The initial state does not affect the final outcome.
2.
Unsteady flows after long times: Statistical long-time behavior is usually independent of the initial state.
3.
Time-periodic flows: Fully converged time-periodic flows are independent of the details of the initialization.
4.
Initialization-sensitive flows: Entire fate of a simulation depends on the initial state (e.g., turbulence, decaying Taylor-Green vortex).

3.5.2 Initial Conditions in LB Simulations

Role of Populations and Macroscopic Fields

•
The mos common situation is that the initial solenoidal (divergence-free) velocity field u0(x)\mathbf{u}_0(\mathbf{x}) is given, but neither the pressure p0(x)p_0(\mathbf{x}) nor the strain rate S0(x)\mathbf{S}_0(\mathbf{x}) are known.
◦
One could be tempted to initialize the populations at equilibrium:
fi(x,t=0)=fieq(ρ,u0(x)).f_i(\mathbf{x}, t = 0) = f_i^{\mathrm{eq}} (\rho, \mathbf{u}_0(\mathbf{x})).
▪
ρ\rho: some initial constant density.
◦
This will lead to an inconsistent state because the Poisson equation for the initial pressure, Δp0=−ρ∂βuα∂αuβ∣t=0\Delta p_0 = -\rho \partial_\beta u_\alpha \partial_\alpha u_\beta |_{t=0}, is not satisfied.
◦
If we assume for a moment that the initial pressure p0(x)p_0(\mathbf{x}) is known:
fi(x,t=0)=fieq(ρ0(x),u0(x)).f_i(\mathbf{x}, t = 0) = f_i^{\mathrm{eq}} (\rho_0(\mathbf{x}), \mathbf{u}_0(\mathbf{x})).
▪
where:
ρ0(x)=ρˉ+p0(x)−pˉcs2.\rho_0(\mathbf{x}) = \bar{\rho} + \frac{p_0(\mathbf{x}) - \bar{p}}{c_s^2}.
▪
The choice of ρˉ\bar{\rho} is arbitrary as it is a mere scaling factor in all euqations.
•
The non-equilibrium populations is needed:
◦
For given velocity gradients and low fluid velocities, the BGK non-equilibrium can be approximated as:
fineq≃−wiτρcs2Qiαβ∂αuβ.f_i^{\mathrm{neq}} \simeq -w_i \frac{\tau \rho}{c_s^2} Q_{i\alpha\beta} \partial_\alpha u_\beta.
▪
where:
Qiαβ=ciαciβ−cs2δαβ.Q_{i\alpha\beta} = c_{i\alpha} c_{i\beta} - c_s^2 \delta_{\alpha\beta}.
•
General initialization:
fi(x,t=0)=fieq(ρ0(x),u0(x))+fineq(ρ0(x),S0(x)).f_i(\mathbf{x}, t = 0) = f_i^{\mathrm{eq}} \big(\rho_0(\mathbf{x}), \mathbf{u}_0(\mathbf{x}) \big) + f_i^{\mathrm{neq}} \big(\rho_0(\mathbf{x}), \mathbf{S}_0(\mathbf{x}) \big).

Chicken or Egg? Order of Collision and Propagation

•
The equilibrium distribution, which is used for collision, is calculated using the known velocity and pressure fields.
◦
Initializing a simulation is basically the inverse of the velocity moments computation.
◦
Populations after initialization assume state compatible with post-propagation, pre-collision. ⇒ Initialization has to be followed by collision rather than propagation.

Consistent Initialization via Modified LB Scheme

•
Mei et al. algorithm:
1.
Initialize the populations fif_i, e.g. with the initial velocity u0(x)\mathbf{u}_0(\mathbf{x}) according to fi(x,t=0)=fieq(ρ,u0(x))f_i(\mathbf{x}, t = 0) = f_i^{\mathrm{eq}} (\rho, \mathbf{u}_0(\mathbf{x})).
2.
Compute the local density: ρ(x)=∑ifi(x)\rho(\mathbf{x}) = \sum_i f_i(\mathbf{x}).
3.
Perform collision using the modified incompressible equilibrium distributions fieq(ρ(x),u0(x))f_i^\mathrm{eq}\big( \rho(\mathbf{x}), \mathbf{u}_0(\mathbf{x}) \big). ⇒ Take the updated density field from step 2 but keep the initial velocity u0(x)\mathbf{u}_0(\mathbf{x}) rather than recomputing it.
4.
Propagate
5.
Go back to step 2 and iterate until the populations fif_i have converged to a user-defined degree. ⇒ It is important to end this algorithm with propagation rather than collision; otherwise the non-equilibrium would be incorrect.
◦
The above algorithm effectively solves:
▪
Advection-diffusion equation for the density ρ\rho as a passive scalar field.
▪
Poisson equation for the pressure pp.
◦
A disadvantage of this initialization approach is that an initial force field is not taken into account.
◦
For the BGK model with ω=1/τ=1/Δt\omega = 1/\tau = 1/\Delta t, ensures velocity field matches input.
◦
Produces consistent initial populations including non-equilibrium part.