LBM Flow Solver

Background & Mathematical Foundation

1. Why the Lattice Boltzmann Method for Hemodynamics?

Traditional CFD discretizes the Navier–Stokes equations directly on body-fitted meshes. That approach is excellent for well-defined engineering geometries but carries two burdens that matter enormously for blood flow: mesh generation becomes a bottleneck when the vessel geometry is drawn, segmented from a scan, or changes interactively, and parallelization requires domain decomposition that does not map cleanly onto GPU fragment shaders. The Lattice Boltzmann Method (LBM) was developed in the late 1980s and early 1990s (Higuera & Jiménez 1989; Qian, d'Humières, Lallemand 1992) and has since become a workhorse for vascular CFD because:

• Arbitrary geometry is trivial. Walls are simply flagged cells on a uniform Cartesian grid. The "bounce-back" boundary condition — at a wall cell, swap each distribution function with its opposite — enforces no-slip without any body-fitted mesh.
• Perfect parallelism. Each cell updates using only its immediate neighbours, with no global operators (no pressure Poisson solve, no Riemann problem). This maps directly onto a fragment shader: one texel = one fluid cell, one draw call = one timestep.
• Kinetic foundation. LBM evolves discrete particle distribution functions $f_i(\mathbf{x}, t)$ rather than the macroscopic $\rho, \mathbf{u}, p$ fields. Non-Newtonian rheology, multiphase flow, and thermal coupling all insert naturally at the kinetic level.
• Recovers Navier–Stokes. A Chapman–Enskog multi-scale expansion of the lattice Boltzmann equation shows that in the limit of small Knudsen number and low Mach number, it recovers the incompressible Navier–Stokes equations to second order in both space and time (Succi 2001).

2. The D2Q9 Lattice

In two dimensions we discretize velocity space with nine lattice vectors — one at rest plus eight directions connecting a cell to its Moore neighbourhood:

Each cell stores nine distribution functions $f_i(\mathbf{x}, t)$, $i = 0 \ldots 8$, representing the density of particles at position $\mathbf{x}$ moving with velocity $\mathbf{c}_i$. The associated lattice weights $w_i$ follow from requiring isotropy and the correct moments of the Gauss–Hermite quadrature:

Macroscopic density and momentum emerge as moments of the distribution:

3. The Lattice Boltzmann Equation: Stream + Collide

Each timestep splits into two substeps. In collision, distributions at each cell relax toward local equilibrium; in streaming, they translate to their neighbours along their respective velocity directions.

3a. BGK Collision

The Bhatnagar–Gross–Krook (BGK) single-relaxation-time collision operator is (Qian et al. 1992):

The discrete equilibrium $f_i^{eq}$ is the Maxwell–Boltzmann distribution expanded to second order in velocity:

The dimensionless relaxation time $\tau$ controls how quickly distributions relax toward equilibrium; it directly sets the kinematic viscosity.

3b. Streaming

Post-collision distributions translate one cell in the direction of their lattice velocity:

This simulator uses the pull formulation: each cell fetches $f_i$ from its upstream neighbour at $\mathbf{x} - \mathbf{c}_i$. Pull streaming maps to a single fragment-shader pass with nine texel fetches per cell, which is exactly what WebGL2 is good at.

4. Recovery of Incompressible Navier–Stokes

The lattice Boltzmann equation governing the continuous-time, continuous-space evolution of the distributions, before discretization onto a lattice, is:

We now derive the macroscopic equations this evolves in the low-Mach, small-gradient limit using the Chapman–Enskog expansion — a multi-scale technique originally developed for the continuous Boltzmann equation (Chapman 1916; Enskog 1917) and adapted to LBM by Frisch, Hasslacher & Pomeau (1986) and later rigorized by Qian et al. (1992).

4a. Multi-scale expansion

Introduce a formal smallness parameter $\varepsilon$ (proportional to the Knudsen number, $\text{Kn} = \ell / L$ = mean free path over macroscopic length) and expand both the distribution and derivatives at multiple scales:

Here $t_1$ is the advective (convective) timescale and $t_2$ is the diffusive (viscous) timescale — these must be kept separate because viscous effects are $O(\varepsilon^2)$ smaller than advection. Substituting into the Boltzmann equation and matching orders of $\varepsilon$ yields a hierarchy of equations.

4b. Order-by-order solution

$O(\varepsilon^0)$: $\;f_i^{(0)} = f_i^{eq}$. The zeroth-order distribution is the equilibrium itself — this is the definition of local thermodynamic equilibrium.

$O(\varepsilon^1)$:

This gives the first-order correction $f_i^{(1)} = -\tau\!\left(\partial_{t_1} + \mathbf{c}_i\cdot\nabla_1\right)\!f_i^{eq}$.

$O(\varepsilon^2)$:

4c. Moments yield continuity and momentum

Summing over all lattice velocities ($\sum_i \cdot$) gives the macroscopic continuity equation. Using the moment constraints $\sum_i f_i^{eq} = \rho$ and $\sum_i f_i^{(k)} = 0$ for $k \geq 1$ (required to preserve $\rho$):

Multiplying by $\mathbf{c}_i$ before summing ($\sum_i \mathbf{c}_i \cdot$) gives the momentum equation. The zeroth moment of $\mathbf{c}_i \mathbf{c}_i f_i^{eq}$ is the equilibrium momentum-flux tensor $\Pi^{(0)}_{\alpha\beta} = \rho u_\alpha u_\beta + p\,\delta_{\alpha\beta}$ (where $p = c_s^2\rho$), and the first-order correction produces a viscous stress:

Combining all moments through $O(\varepsilon^2)$ and dropping $O(\text{Ma}^3)$ terms:

with the identifications:

4d. The $-\tfrac{1}{2}$ correction and why it matters

A subtle but critical point: the discrete-time lattice update differs from the continuous-time Boltzmann equation by a second-order Taylor error. Expanding $f_i(\mathbf{x} + \mathbf{c}_i\Delta t, t + \Delta t)$ in $\Delta t$ and comparing with the continuous equation produces an effective relaxation time $\tau_{\text{eff}} = \tau - \tfrac{1}{2}$. This half-step shift is what produces the $(\tau - \tfrac{1}{2})$ in the viscosity formula — not $\tau$ alone. The constraint $\tau > \tfrac{1}{2}$ is thus a physical requirement: negative effective viscosity is unstable. At $\tau \to \tfrac{1}{2}$ the scheme approaches the inviscid (Euler) limit but loses stability. Practically we keep $\tau \in [0.55, 2.5]$ in this simulator.

4e. Weak compressibility

LBM recovers the compressible Navier–Stokes equations with an isothermal equation of state $p = c_s^2 \rho$. Incompressibility is only reached in the limit $\text{Ma} = |\mathbf{u}|/c_s \to 0$. The relative density variation scales as $\Delta\rho/\rho \sim \text{Ma}^2$, so at lattice Mach $0.1$ the compressibility error is about $1\%$ — which is why LBM is described as weakly compressible. This simulator keeps lattice Mach comfortably below $0.1$ (adaptive $U_{\text{lat}} \in [0.02, 0.04]$), so compressibility error is well under $1\%$ for all supported inlet velocities.

5. Boundary Conditions

5a. No-slip Walls: Full-way Bounce-back

At a wall cell, the post-streaming distributions are swapped pairwise with their opposites:

On the next timestep, fluid cells adjacent to the wall stream in these reversed distributions, producing zero velocity at the halfway point between the last fluid cell and the wall cell — this is the half-way bounce-back interpretation, which is second-order accurate at the wall interface (Ladd 1994).

5b. Velocity Inlet (Plug Flow)

On inlet cells we prescribe $\mathbf{u} = (U_{\text{in}}, 0)$, set $\rho = 1$, and overwrite all nine distributions with their equilibrium values $f_i = f_i^{eq}(\rho, \mathbf{u})$. The bounce-back walls downstream enforce the no-slip condition and the flow develops the characteristic Poiseuille profile over the entry length $L_e \approx 0.06\, \text{Re}\, D$. Students can observe this entry-flow development directly.

5c. Pressure Outlet / Inlet

When pressure is prescribed, density is fixed at the target value and velocity is extrapolated from the nearest interior fluid cell. In LBM, pressure and density are linked through the isothermal equation of state:

6. Non-Newtonian Blood: Carreau–Yasuda

Blood is a non-Newtonian shear-thinning fluid due to reversible erythrocyte aggregation at low shear. The Carreau–Yasuda model (Yasuda et al. 1981) captures this over the full physiological range:

Recommended parameters for human blood at 37 °C (Cho & Kensey 1991):

$\mu_0$

56 mPa·s (zero-shear viscosity, dominated by RBC aggregates)

$\mu_\infty$

3.45 mPa·s (infinite-shear viscosity, dispersed RBCs)

$\lambda$

3.313 s (relaxation time)

$n$

0.3568 (shear-thinning exponent)

$a$

2.0 (Yasuda transition exponent)

In LBM we need the local shear rate $\dot\gamma$ at each cell. The Chapman–Enskog expansion provides this directly from the non-equilibrium part of the distribution — the momentum-flux tensor $\Pi_{\alpha\beta}^{\text{neq}}$ — without any finite differences of $\mathbf{u}$:

The Chapman–Enskog shear rate is then fed back into Carreau–Yasuda to produce a local $\mu(\dot\gamma)$, which sets a local relaxation time $\tau_{\text{local}} = \tfrac{1}{2} + 3\nu_{\text{local}}$. This closes the loop: each cell collides with its own effective viscosity every step.

7. Lattice-to-Physical Unit Mapping

The simulator runs internally in lattice units (everything dimensionless, grid spacing $\Delta x = 1$, timestep $\Delta t = 1$). Physical SI quantities are recovered on the fly through a chosen length scale and lattice velocity:

$\rho_0$

1060 kg/m³ (blood density, fixed)

$L$

1 cm (characteristic vessel width = grid height $N_y$)

$\Delta x$

$L / N_y$ m per lattice cell

$U_{\text{lat}}$

0.05 (fixed, keeps lattice Mach $\ll$ 1)

$\Delta t$

$U_{\text{lat}}\cdot \Delta x / U_{\text{phys}}$ s per lattice step

$\nu_{\text{lat}}$

$\nu_{\text{phys}}\cdot \Delta t / \Delta x^2$

$\tau$

$\tfrac{1}{2} + 3\,\nu_{\text{lat}}$

$p_{\text{phys}}$

$\rho_0\, (\Delta x / \Delta t)^2 \cdot (\rho_{\text{lat}} - 1) / 3$ Pa

8. Reynolds Regimes in the Vasculature

Reynolds number is the dimensionless ratio of inertial to viscous forces — it determines whether flow is dominated by friction (low Re, Stokes) or momentum (high Re, inertial, turbulent):

Typical values across the circulation:

This simulator's default preset (sphere in channel, U=0.30 m/s, μ=3.5 mPa·s, L=1 cm) gives Re ≈ 909, squarely in the arterial regime. Drop μ to 1 mPa·s and raise U to 1.0 m/s to push toward vortex shedding around Re ≈ 3000.

9. Preset Scenarios

• Empty rectangular channel — baseline Poiseuille flow. Plug inlet develops into parabolic profile over the entry length. Good for sanity checks and for observing the centerline velocity reach $\tfrac{3}{2}U_{\text{in}}$ at full development.
• Sphere (disk) in channel — classic bluff-body wake. At low Re the flow is symmetric; above Re ≈ 45 (for a cylinder) the von Kármán vortex street emerges. This simulator reproduces the transition visibly. Inviscid theory predicts zero drag (d'Alembert's paradox); the finite numerical viscosity of LBM — plus the physical viscosity you dial in — gives a realistic drag and wake.
• Cube (square block) — sharp-edged bluff body. Separation is forced at the upstream corners, producing a wider wake and higher pressure drag than the cylinder at the same frontal area. Useful for illustrating the geometric basis of flow separation.
• Stenosis (50% symmetric) — a cosine-shaped symmetric narrowing representative of a concentric atherosclerotic lesion. By continuity, the throat velocity is roughly doubled; by Bernoulli, the throat pressure drops accordingly. Downstream, flow decelerates and can separate on the expansion wall — the post-stenotic recirculation zone that promotes platelet adhesion. Clinically, a ≥70% area stenosis is considered hemodynamically significant.
• Atherosclerotic plaque (asymmetric) — one-sided Gaussian-shaped lesion on the lower wall with a steep distal edge, representative of a shoulder plaque. The asymmetry produces a skewed jet, an elongated recirculation on the healthy-wall side, and sustained wall-shear gradients at the distal shoulder — geometries associated with plaque rupture risk.
• Saccular aneurysm — a dome carved into the upper wall. Main-flow fluid convects past the ostium while slow recirculating flow fills the dome. Dome-to-parent-vessel flow exchange and intra-aneurysmal shear are the quantities of interest for rupture-risk modelling (Cebral et al. 2005).
• Bifurcation (Y-split) — idealized carotid-type bifurcation. Flow divides into two daughters of equal diameter. The outer wall of each daughter experiences a stagnation-like stress concentration (the apex); the inner walls see accelerated flow. This is where early atherosclerotic lesions preferentially develop.
• Trifurcation — three-way division (less common anatomically; approximates the celiac/SMA region). Useful for comparing flow splitting behaviour with the bifurcation.

10. Numerical Limitations & Known Behavior

• Weakly compressible. LBM recovers weakly compressible Navier–Stokes; strict incompressibility is only reached in the limit $U_{\text{lat}} \to 0$. This simulator fixes $U_{\text{lat}} = 0.05$, giving a lattice Mach of 0.087 — compressibility error is below 1% for these flows.
• Half-way bounce-back is first-order in pressure. Standard bounce-back is second-order accurate for velocity but first-order for pressure near walls. For high-fidelity pressure fields along curved walls, higher-order BCs (Bouzidi, interpolated) are needed. For teaching, bounce-back is the right trade-off: simple, stable, robust against user-drawn geometries.
• Uniform grid. True adaptive mesh refinement in LBM requires multi-grid rescaling at interfaces (Filippova & Hänel 1998; Lagrava et al. 2012) — a substantial implementation. This simulator uses a uniform fine grid (up to 640×240 = 153,600 cells); at Ultra resolution, fine geometry features resolve cleanly.
• Plug inlet not fully developed. The inlet BC is plug flow, not parabolic Poiseuille. Flow develops over the entry length downstream of the inlet. For a channel of height $D$ and Reynolds Re, $L_e \approx 0.06\,\text{Re}\,D$, which at Re = 1000 and D = 1 cm is 60 cm — longer than the visible domain. Therefore flow in the domain is developing, not fully developed. Don't compare centerline velocity with the $1.5 U_\text{in}$ Poiseuille value unless Re is small.
• Non-Newtonian uses last-step τ. The Chapman–Enskog shear rate uses the previous step's relaxation time, not an iteratively-converged value. For steady state this introduces a small lag but converges correctly.
• 2D only. This solver is 2D. Real vascular flow is three-dimensional; key features (secondary flow in curved vessels, Dean vortices, fully developed 3D helical flow) are not representable. Use this for intuition, not quantitative clinical prediction.