Search

[LBM Study] 2. Analysis of the Lattice Boltzmann Equation

Category
Study
Kewords
LBM
Kruger
3 more properties
Table of Contents

2.1 Chapman-Enskog Analysis

2.1.1 The Perturbation Expansion

•
Key Concept: The distribution function fif_i is expanded around the equilibrium distribution fieqf_i^\mathrm{eq} with Knudsen number Kn\mathrm{Kn} as the expansion parameter:
fi=fieq+ϵfi(1)+ϵ2fi(2)+⋯ .f_i = f_i^\mathrm{eq} + \epsilon f_i^{(1)} + \epsilon^2 f_i^{(2)} + \cdots.
◦
where ϵn\epsilon^n is a label indicating the order of Knn\mathrm{Kn}^n.
Important Assumption: In perturbation analysis, we assume that only the two lowest orders in Kn\mathrm{Kn} are sufficient to find the NSE.
•
LBE with BGK collision operator:
fi(x+ciΔt,t+Δt)−fi(x,t)=−Δtτ[fi(x,t)−fieq(x,t)].f_i(\mathbf{x} + \mathbf{c}_i\Delta t, t + \Delta t) - f_i(\mathbf{x}, t) = -\frac{\Delta t}{\tau}[f_i(\mathbf{x}, t) - f_i^\mathrm{eq}(\mathbf{x}, t)].
•
Mass and Momentum Conservation (Solvability Conditions):
∑ifineq=0,∑icifineq=0.\sum_i f_i^\mathrm{neq} = 0, \quad \sum_i \mathbf{c}_i f_i^\mathrm{neq} = \mathbf{0}.
•
These can be assumed to hold individually at each order:
∑ifi(n)=0 and ∑icifi(n)=0 for all n≥1.\sum_i f_i^{(n)} = 0 \text{ and } \sum_i \mathbf{c}_i f_i^{(n)} = \mathbf{0} \text{ for all } n \geq 1.

2.1.2 Taylor Expansion, Perturbation, and Separation

•
Taylor expanding the LBE gives:
Δt(∂t+ciα∂α)fi+Δt22(∂t+ciα∂α)2fi+O(Δt3)=−Δtτfineq.\Delta t(\partial_t + c_{i\alpha}\partial_\alpha)f_i + \frac{\Delta t^2}{2}(\partial_t + c_{i\alpha}\partial_\alpha)^2f_i + O(\Delta t^3) = -\frac{\Delta t}{\tau}f_i^\mathrm{neq}.
•
After neglecting third and higher order terms and eliminating second-order derivative terms:
Δt(∂t+ciα∂α)fi=−Δtτfineq+Δt(∂t+ciα∂α)Δt2τfineq.\Delta t(\partial_t + c_{i\alpha}\partial_\alpha)f_i = -\frac{\Delta t}{\tau}f_i^\mathrm{neq} + \Delta t(\partial_t + c_{i\alpha}\partial_\alpha)\frac{\Delta t}{2\tau}f_i^\mathrm{neq}.
•
Expansion of time and space derivatives:
Δt∂tfi=Δt(ϵ∂t(1)fi+ϵ2∂t(2)fi+⋯),Δtciα∂αfi=Δt(ϵciα∂α(1)fi).\Delta t\partial_t f_i = \Delta t\bigl(\epsilon\partial_t^{(1)}f_i + \epsilon^2\partial_t^{(2)}f_i + \cdots \bigl), \quad\quad \Delta t c_{i\alpha} \partial_\alpha f_i = \Delta t\bigl(\epsilon c_{i\alpha}\partial_\alpha^{(1)}f_i \bigl).
•
Separating by orders of Kn\mathrm{Kn}:
O(ϵ):(∂t(1)+ciα∂α(1))fieq=−1τfi(1),O(ϵ2):∂t(2)fieq+(∂t(1)+ciα∂α(1))(1−Δt2τ)fi(1)=−1τfi(2).\begin{aligned} O(\epsilon): \quad(\partial_t^{(1)} + c_{i\alpha}\partial_\alpha^{(1)})f_i^\mathrm{eq} = -\frac{1}{\tau}f_i^{(1)}, \\\\ O(\epsilon^2): \quad\partial_t^{(2)}f_i^\mathrm{eq} + (\partial_t^{(1)} + c_{i\alpha}\partial_\alpha^{(1)})(1 - \frac{\Delta t}{2\tau})f_i^{(1)} = -\frac{1}{\tau}f_i^{(2)}. \end{aligned}

2.1.3 Moments and Recombination

•
Equilibrium moments:
Παβeq=∑iciαciβfieq=ρuαuβ+ρcs2δαβ,Παβγeq=∑iciαciβciγfieq=ρcs2(uαδβγ+uβδαγ+uγδαβ),Παβ(1)=∑iciαciβfi(1).\begin{aligned} \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},\\\\ \Pi_{\alpha\beta\gamma}^{\mathrm{eq}} &= \sum_i c_{i\alpha} c_{i\beta} c_{i\gamma} f_i^{\mathrm{eq}} = \rho c_s^2 (u_\alpha \delta_{\beta\gamma} + u_\beta \delta_{\alpha\gamma} + u_\gamma \delta_{\alpha\beta}), \\\\ \Pi_{\alpha\beta}^{(1)} &= \sum_i c_{i\alpha} c_{i\beta} f_i^{(1)}. \end{aligned}
•
Taking the 0th to 2nd moments of O(ϵ)O(\epsilon) Kn\mathrm{Kn} equation yields the O(ϵ)O(\epsilon) moment equations:
∂t(1)ρ+∂γ(1)(ρuγ)=0,∂t(1)(ρuα)+∂β(1)Παβeq=0,∂t(1)Παβeq+∂γ(1)Παβγeq=−1τΠαβ(1).\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} = 0, \\\\ \partial_t^{(1)}\Pi_{\alpha\beta}^\mathrm{eq} + \partial_\gamma^{(1)}\Pi_{\alpha\beta\gamma}^\mathrm{eq} = -\frac{1}{\tau}\Pi_{\alpha\beta}^{(1)}. \end{aligned}
•
Taking the 0th and 1st moments of O(ϵ2)O(\epsilon^2) Kn\mathrm{Kn} equation yields the O(ϵ2)O(\epsilon^2) moment 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}
•
Assembling the mass and momentum equations from their O(ϵ)O(\epsilon) and O(ϵ2)O(\epsilon^2) component equations, we find
(ϵ∂t(1)+ϵ2∂t(2))ρ+ϵ∂β(1)(ρuγ)=0,(ϵ∂t(1)+ϵ2∂t(2))(ρuα)+ϵ∂β(1)Παβeq=−ϵ2∂β(1)(1−Δt2τ)Παβ(1).\begin{aligned}\left( \epsilon \partial_t^{(1)} + \epsilon^2 \partial_t^{(2)} \right) \rho + \epsilon \partial_\beta^{(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^2 \partial_\beta^{(1)} \left( 1 - \frac{\Delta t}{2\tau} \right) \Pi_{\alpha\beta}^{(1)}.\end{aligned}
◦
Reversing the derivative expansions, these equations become the continuity equation and a momentum conservation equation.
◦
With an as-of-yet unknown viscous stress tensor:
σαβ=−(1−Δt2τ)Παβ(1).\sigma_{\alpha\beta} = -\left( 1 - \frac{\Delta t}{2\tau} \right) \Pi_{\alpha\beta}^{(1)}.
Through Chapman-Enskog analysis, we can finally find Perturbation Moment:
Παβ(1)=−ρcs2τ(∂β(1)uα+∂α(1)uβ)+τ∂γ(1)(ρuαuβuγ).\Pi_{\alpha\beta}^{(1)} = -\rho c_s^2 \tau \left( \partial_\beta^{(1)} u_\alpha + \partial_\alpha^{(1)} u_\beta \right) + \tau \partial_\gamma^{(1)} (\rho u_\alpha u_\beta u_\gamma).
The second term is an error term arising from the lack of a correct O(u3)O(u^3) term in the equilibrium distribution fieqf_i^\mathrm{eq}.

2.1.4 Macroscopic Equations

•
Final Macroscopic Equations can be obtained by 1) inserting Perturbation Moment Παβ(1)\Pi_{\alpha\beta}^{(1)} into mass and momentum equations from their O(ϵ)O(\epsilon) and O(ϵ2)O(\epsilon^2) component equations, 2) reversing the derivative expansion:
∂tρ+∂γ(ρuγ)=0,∂t(ρuα)+∂β(ρuαuβ)=−∂αp+∂β[η(∂βuα+∂αuβ)]\begin{aligned}\partial_t \rho + \partial_\gamma (\rho u_\gamma) &= 0, \\\partial_t (\rho u_\alpha) + \partial_\beta (\rho u_\alpha u_\beta) &= -\partial_\alpha p + \partial_\beta \left[ \eta (\partial_\beta u_\alpha + \partial_\alpha u_\beta) \right]\end{aligned}
◦
Pressure: p=ρcs2p = \rho c_s^2.
◦
Kinematic viscosity: ν=cs2(τ−Δt2)\nu = c_s^2(\tau - \frac{\Delta t}{2}).
◦
Dynamic viscosity: η=ρν=ρcs2(τ−Δt2)\eta = \rho \nu = \rho c_s^2(\tau - \frac{\Delta t}{2}).
◦
Bulk viscosity: ηB=23η\eta_B = \frac{2}{3}\eta.
•
Bulk viscosity differences:
◦
In monatomic kinetic theory, the bulk viscosity ηB\eta_{B} is normally zero.
◦
However, in LBE, the bulk viscosity to be of the order of the shear viscosity η\eta.
◦
Difference caused by isothermal equation of state usage, fundamentally incompatible with monatomic assumption.
Stability Condition (for positive viscosity):
•
Necessary condition for stability: τ/Δt≥1/2\tau/\Delta t \geq 1/2
•
The same condition was found in the discrete LBGK equation’s behavior; τ/Δt<1/2\tau / \Delta t < 1/2 would lead to a divergent under-relaxation.

2.2 Discussion of the Chapman-Enskog Analysis

2.2.1 Dependence of Velocity Moments

•
Third-Order Moment Problem: In standard velocity sets (D2Q9, D3Q19, etc.), third-order moments are not independent:
Πxxxeq=∑icix3fieq=(ΔxΔt)2∑icixfieq=(ΔxΔt)2Πxeq.\Pi_{xxx}^{eq} = \sum_i c_{ix}^3 f_i^{eq} = \left(\frac{\Delta x}{\Delta t}\right)^2 \sum_i c_{ix} f_i^{eq} = \left(\frac{\Delta x}{\Delta t}\right)^2 \Pi_x^{eq}.
◦
This results in O(u3)O(u^3) errors in the macroscopic momentum equation (stress tensor).

2.2.2 Time Scale Interpretation

•
Common Misinterpretation:
∂t=ϵ∂t(1)+ϵ2∂t(2)+⋯.∂_t=ϵ∂_t^{(1)}+ϵ^2∂_t^{(2)}+⋯.
◦
Time derivative expansion is often misinterpreted as a decomposition into different time scales. → Viewed as "clocks ticking at different speeds".
◦
This interpretation can lead to false conclusions.
•
Example: Steady Poiseuille Flow
◦
Steady Poiseuille flow: ∇p=∇⋅σ\nabla p = \nabla \cdot \sigma (pressure gradient balanced by viscous stress).
◦
Correct physics: ∂t(uα)=0\partial_t(u_\alpha) = 0 (velocity doesn't change with time).
◦
False conclusion from time scale interpretation:
▪
∂t(1)(uα)=0\partial_t^{(1)}(u_\alpha) = 0 and ∂t(2)(uα)=0\partial_t^{(2)}(u_\alpha) = 0 independently. → This would lead to: ∇p=0\nabla p = 0 and ∇⋅σ=0\nabla \cdot \sigma = 0 separately.
▪
Result: No flow at all! (contradicts physical reality).
◦
Correct Understanding:
▪
∂t(n)∂_t^{(n)} are NOT actual time derivatives. → They are perturbation expansion terms at different orders in ϵ\epsilon. → Only their sum equals the physical time derivative!
▪
Correct relation: ∂t(uα)=∂t(1)(uα)+∂t(2)(uα)=0.\partial_t(u_\alpha) = \partial_t^{(1)}(u_\alpha) + \partial_t^{(2)}(u_\alpha) = 0.
▪
This gives: ∇p=∇⋅σ\nabla p = \nabla \cdot \sigma (physically correct balance).

2.3 Alternative Equilibrium Models

2.3.1 Linear Fluid Flow

•
Discrete Equilibrium Distribution:
fieq=wiρ(1+ci⋅ucs2+(ci⋅u)22cs4−u⋅u2cs2)=wiρ(1+ciαuαcs2+uαuβ(ciαciβ−cs2δαβ)2cs4).f_i^{\mathrm{eq}} = w_i ρ \left( 1 + \frac{\mathbf{c}_i \cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i \cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u} \cdot \mathbf{u}}{2c_s^2} \right) = w_i ρ \left( 1 + \frac{c_{iα} u_α}{c_s^2} + \frac{u_α u_β (c_{iα} c_{iβ} - c_s^2 δ_{αβ})}{2c_s^4} \right).
•
Linearized Equilibrium Distribution (used in acoustics):
fieq=wi(ρ+ρ0ciαuαcs2).f_i^\mathrm{eq} = w_i\left(\rho + \rho_0\frac{c_{i\alpha}u_\alpha}{c_s^2}\right).
◦
ρ0\rho_0: the rest state density.
•
This produces linearized macroscopic equations (Chapman-Enskog analysis results):
∂tρ+ρ0∂αuα=0,ρ0∂tuα=−∂αp+η∂β(∂βuα+∂αuβ).\begin{aligned} \partial_t \rho + \rho_0 \partial_\alpha u_\alpha = 0, \\\\ \rho_0 \partial_t u_\alpha = -\partial_\alpha p + \eta \partial_\beta \left( \partial_\beta u_\alpha + \partial_\alpha u_\beta \right). \end{aligned}
◦
Pressure: p=ρcs2p = \rho c_s^2.
◦
Dynamic shear viscosity: η=ρ0cs2(τ−Δt/2)\eta = \rho_0 c_s^2 (\tau - \Delta t/2).

2.3.2 Incompressible Flow

•
Incompressible Equilibrium Distribution:
◦
Using the fact that ρ′/ρ0=O(Ma2)\rho'/\rho_0 = O(Ma^2) in steady flow:
fieq=wiρ+wiρ0(ciαuαcs2+uαuβ(ciαciβ−cs2δαβ)2cs4).f_i^\mathrm{eq} = w_i\rho + w_i\rho_0\left(\frac{c_{i\alpha}u_\alpha}{c_s^2} + \frac{u_\alpha u_\beta(c_{i\alpha}c_{i\beta} - c_s^2\delta_{\alpha\beta})}{2c_s^4}\right).
•
This recovers the exact incompressible Navier-Stokes equations in steady state.

2.4 Stability

2.4.1 Stability Analysis

•
Courant number insufficient for LBM: ⇒ Unlike standard CFD where C=∣u∣ΔtΔx≤1C = \frac{|u|\Delta t}{\Delta x} \leq 1 often determines stability, LBM has additional degrees of freedom (relaxation times) that make Courant number alone inadequate.
•
LBM stability analysis requires inversion of q×qq \times q matrix (where qq is number of populations) due to one equation per velocity direction ci\mathbf{c}_i, making it more complex than standard CFD.
•
Stability map exists: ∣u∣max⁡(τ)|u|_{\max}(\tau) defines maximum achievable velocity magnitude for given relaxation time τ\tau before instability sets in.
◦
High Reynolds number dilemma: for Re=∣u∣NΔxνRe = \frac{|u|N\Delta x}{\nu}, achieving high ReRe requires compromise between maximum velocity and minimum viscosity, both constrained by stability limits.

2.4.2 BGK Stability

•
Stability Conditions for BGK Collision Operator:
◦
Sufficient stability condition is the non-negativity of all equilibrium populations. → τ/Δt>1/2τ/Δt > 1/2 & fieq≥0f_i^\mathrm{eq} \ge 0 for all ii
◦
Since the equilibrium populations are functions of the velocity u\mathbf{u}, this can be expressed as a sufficient stability condition for the velocity u\mathbf{u}.
▪
For D2Q9D2Q9, D3Q15D3Q15, D3Q19D3Q19, D3Q27D3Q27:
∣umax∣<13ΔxΔt≈0.577ΔxΔt.|\mathbf{u}_\mathrm{max}| < \sqrt{\frac{1}{3}}\frac{\Delta x}{\Delta t} \approx 0.577\frac{\Delta x}{\Delta t}.
◦
Optimal stability condition is the non-negativity of the rest equilibrium population. → τ/Δt>1τ/Δt > 1 & f0eq>0f_0^\mathrm{eq} \gt 0
◦
Velocity magnitude condition:
∣u∣<23ΔxΔt.|\mathbf{u}| < \sqrt{\frac{2}{3}}\frac{\Delta x}{\Delta t}.
Stability in Practical Simulations:
•
In real simulations with boundaries, as τ/Δtτ/Δt → 12\frac{1}{2}, ∣umax∣|u_\mathrm{max}| should approaches zero:
∣umax∣(τ)=8(τΔt−12)ΔxΔt for τΔt<0.55.|\mathbf{u}_\mathrm{max}|(\tau) = 8\left(\frac{\tau}{\Delta t} - \frac{1}{2}\right)\frac{\Delta x}{\Delta t} \quad\quad \text{ for } \quad \frac{\tau}{\Delta t} < 0.55.
•
As a guideline to find stable parameters for small viscosities:
◦
Start with the sufficient stability condition for all τ/Δt>12.\tau / \Delta t > \frac{1}{2}.
◦
If simulations are unstable, then perform a few simulations with different values of τ\tau and u\mathbf{u} to find an empirical relation ∣umax∣(τ)|u_\mathrm{max}|(\tau).

2.4.3 Stability for Advanced Collision Operators

•
TRT (Two-Relaxation-Time) Model ⇒ In the TRT framework, there is a certain combination of τ+\tau^+ and τ−\tau^- that governs the stability and accuracy of simulations:
Λ=(τ+Δt−12)(τ−Δt−12).\Lambda = \left(\frac{\tau^+}{\Delta t} - \frac{1}{2}\right)\left(\frac{\tau^-}{\Delta t} - \frac{1}{2}\right).
◦
A recommended choice is Λ=1/4\Lambda=1/4. → This corresponds to τ/Δt=1\tau / \Delta t = 1 in the BGK case, and allows the same optimal stability from which the velocity condition in ∣u∣<2/3ΔxΔt|\mathbf{u}| < \sqrt{2/3}\Delta x \Delta t.
◦
For any valude of τ+\tau^+, one can always select the free parameter τ−\tau^- such that Λ=1/4\Lambda=1/4.
◦
The advantage of the TRT model is therefore that the stability condiiton and the kinematic viscosity ν=cs2(τ−Δt2)\nu = c_s^2(\tau - \frac{\Delta t}{2}) are decoupled.

2.4.4 Stability Improvement Guidelines

•
Stability Improvement Process:
1.
Start with BGK collision operator with any τ\tau and Reynolds number Re\text{Re}:
Re=∣u∣NΔxcs2(τ−Δt2).\text{Re} = \frac{|u|N\Delta x}{c_s^2\left(\tau - \frac{\Delta t}{2}\right)}.
•
NN is the number of lattice nodes along a characteristic length scale l=NΔxl=N \Delta x.
2.
Set τ=Δt\tau = \Delta t and adjust ∣u∣|\mathbf{u}|to matchthe Reynolds number.
3.
Check if ∣u∣≥cs|\mathbf{u}| ≥ c_s:
•
If No → Check if simulation is stable.
•
If Yes → Maintain Re\text{Re}, reduce τ\tau to find new ∣u∣＜cs|\mathbf{u}| ＜ c_s. Check if ∣u∣<∣umax∣(τ)|\mathbf{u}| < |\mathbf{u}_\mathrm{max}|(\tau).
4.
If unstable → Use TRT/MRT with Λ = 1/4.

2.5 Accuracy

2.5.1 Formal Order of Accuracy

•
Truncation error analysis through Taylor expansion is used:
∂u/∂t−ν(∂2u/∂y2)=(Δt/2)(∂2u/∂t2)+ν(Δy2/4!)(∂4u/∂y4)+O(Δt2)+O(Δy4).∂u/∂t - ν(∂²u/∂y²) = (Δt/2)(∂²u/∂t²) + ν(Δy²/4!)(∂⁴u/∂y⁴) + O(Δt²) + O(Δy⁴).
◦
This demonstrates the scheme is first-order in time and second-order in space.

2.5.2 Accuracy Measurement

•
L₂ Error Norm:
ϵq(t):=∑x[qn(x,t)−qa(x,t)]2∑xqa2(x,t).\epsilon_q(t) := \sqrt{\frac{\sum_\mathbf{x}[q_n(\mathbf{x},t) - q_a(\mathbf{x},t)]^2}{\sum_\mathbf{x} q_a^2(\mathbf{x},t)}}.

2.5.3 Numerical Errors

Error Type
Description
Mitigation
Round-off Error
Due to finite precision
Use double precision
Iterative Error
Incomplete convergence
Convergence criterion: ε < 10 ⁻ ⁷
Discretization Error
Main error from discretizing PDEs
Can be controlled through τ selection

2.5.4 Modeling Errors

•
O(u3)O(u^3) error:
◦
Due to insufficient isotropy of standard lattices.
◦
Negligible when Ma2≪1\text{Ma}^2 ≪ 1.
•
Compressibility error:
◦
O(Ma2)O(\text{Ma}^2) magnitude.
◦
Difference from incompressible NSE.

2.5.5 Lattice Boltzmann Accuracy

•
Optimization Conditions:
Condition
BGK τ/Δt
TRT Λ
Application
Cancel O(ε ³)
≈ 0.789
1/12
Advection-dominated
Cancel O(ε ⁴)
≈ 0.908
1/6
Diffusion-dominated
Optimal stability
1.0
1/4
General cases

2.5.6 Accuracy Improvement Guidelines

•
Primary Objective:
◦
Optimize accuracy while respecting non-dimensional groups of the problem.
◦
Match non-dimensional numbers (Reynolds number, Péclet number) between physical problem and LB scheme.
•
Accuracy Improvement Process:
1.
Initial Setup and Parameter Selection
a.
Determine relevant non-dimensional numbers (Reynolds number).
b.
Select grid number NN based on computational resources.
c.
Set final parameter τ\tau as main control variable.
2.
Collision Operator Assessment and Selection
a.
Start with BGK collision operator.
b.
Evaluate τ\tau-sensitivity for velocity condition:
u/(ΔxΔt)≪1.\mathbf{u} / \big( \frac{Δx}{Δt} \big) ≪ 1.
c.
Keep BGK if solutions almost τ\tau-independent or Change to TRT/MRT if significantly τ\tau-dependent.
3.
Accuracy Optimization and Final Adjustment
a.
Assess current accuracy satisfaction.
b.
If unsatisfactory → Increase grid number NN while keeping Reynolds number constant.
c.
Proceed with simulations using optimized parameters.

2.6 Summary

Key Points

•
Chapman-Enskog analysis: Establishes connection between LBE and macroscopic NSE.
•
Macroscopic equations: LBE solves weakly compressible NSE, converging to incompressible NSE as Ma2→0Ma² → 0.
•
Stability: τ>Δt/2τ > Δt/2 is essential condition, with more stringent constraints in practice.

Practical Recommendations

1.
High Reynolds number simulations: Use TRT/MRT recommended.
2.
Accuracy optimization: Choose τ\tau or Λ\Lambda based on problem characteristics.
3.
Stability issues: Set τ=Δtτ = Δt then adjust velocity.
4.
Minimize compressibility error: Maintain Ma<0.1Ma < 0.1.