Surviving the Low-Viscosity Regime Where BGK Blows Up — MRT Lattice Boltzmann Collision

Same lattice, same equilibrium. One model spins. The other shatters into a checkerboard. BGK pushes all nine distributions with one relaxation time. MRT (Multiple Relaxation Time) maps them to moments and pushes each at its own rate.

Shove the $s_\nu$ bar to 1.99 and the viscosity readout falls toward 0. The $s_e$ and $s_q$ bars stay where they are.

The conserved bars $\rho$, $j_x$, $j_y$ sit on zero. BGK ties every bar to one number. That tie is where the checkerboard shows up at low viscosity.

The Weakness of Relaxing Everything in One Beat#

The lattice Boltzmann method (LBM) repeats two steps: distribution functions $f_i$ stream across the lattice, and at each point they relax toward equilibrium (collision). The most common BGK collision is written as

where $f_i$ is the distribution along velocity $\mathbf{c}_i$, $\tau$ is the relaxation time, and $f_i^{\text{eq}}$ is the local equilibrium. Viscosity is fixed by $\tau$ alone.

The trouble starts here. To lower viscosity (to raise the Reynolds number), you must push $\tau$ close to $1/2$. But as $\tau \to 1/2$, the high-order modes that have nothing to do with viscosity barely relax and roam the lattice. When those modes grow near the grid-resolution limit, a checkerboard oscillation appears and the simulation diverges. In BGK, the viscosity knob and the stability knob are tied to the same handle.

Looking at Collision in Moment Space#

Do not relax the nine distributions as they stand. Map them into physically meaningful moments first, then relax. Moments are linear combinations of the distributions.

Here $\mathbf{f} = (f_0,\dots,f_8)^\top$ is the distribution vector, $M$ is the $9\times 9$ transformation matrix, and $\mathbf{m}$ is the moment vector. Collision relaxes toward the equilibrium moments $\mathbf{m}^{\text{eq}}$ in moment space, then maps back to velocity space.

$S = \mathrm{diag}(s_0,\dots,s_8)$ is the diagonal matrix of per-moment relaxation rates. If you set every $s_i$ to the same value $1/\tau$, this collapses exactly back to BGK. MRT does nothing more than solve those diagonal entries moment by moment.

The Nine D2Q9 Moments and the Transformation Matrix M#

On the D2Q9 lattice (two dimensions, nine velocities) the moments are these nine: density $\rho$, energy $e$, energy square $\varepsilon$, momentum $j_x, j_y$, energy flux $q_x, q_y$, and the stress-like $p_{xx}, p_{xy}$.

Among them, $\rho, j_x, j_y$ are conserved quantities unchanged by collision, so their relaxation rates are fixed at $s=0$. The other six are the tunable handles. The standard Lallemand–Luo transformation matrix is built so its rows are orthogonal, so the inverse is simply the transpose divided by the row norms.

The equilibrium moments close in terms of density and momentum alone.

Here $j_x = \rho u_x$, $j_y = \rho u_y$ are the momentum components. The equilibria of the stress moments $p_{xx}, p_{xy}$ couple to velocity gradients, so their relaxation rates set the viscosity.

A Different Rate per Moment — Degrees of Freedom Unrelated to Viscosity#

The shear viscosity is set by the relaxation rate of the stress moments, $s_\nu$.

The other rates $s_e$ (bulk), $s_\varepsilon$, and $s_q$ (energy flux) never enter the viscosity. They are free. Park $s_\nu$ against $2$ and the ghost modes can still sit in the stable band $1.5\sim1.8$. The two handles BGK had lashed together come apart here.

The listing's default viscosity is $\nu=0.004$. Then $s_\nu = 1/(3\nu + 1/2)$, so $3\times 0.004 + 0.5 = 0.512$ and $s_\nu = 1/0.512 = 1.953125$. Shove the bar to 1.99 and $\nu = (1/3)(1/1.99 - 1/2)$. Subtract 0.5 from $1/1.99 \approx 0.5025$ and divide by 3: $\nu \approx 8.4\times 10^{-4}$. Near zero, not zero. At that viscosity BGK also drives $s_e$, $s_\varepsilon$, and $s_q$ to 1.99. The MRT array in the listing pins those three slots at $1.64$, $1.54$, and $1.7$. The vector is [0, 1.64, 1.54, 0, 1.7, 0, 1.7, s_nu, s_nu], so the only zeros are the three conserved slots.

Equilibrium moments close on density and momentum only after the distribution has already been cut at low order. That cut is Maxwell–Boltzmann inside nine velocities.

Collision in Moment Space, Streaming in Velocity Space#

The algorithm flows like this.

1. Compute the macroscopic $\rho, u_x, u_y$ from $f_i$
2. Transform to moments with $\mathbf{m} = M\mathbf{f}$
3. Compute equilibrium moments $\mathbf{m}^{\text{eq}}$ (from density and momentum)
4. Relax per moment: $m_a \leftarrow m_a - s_a(m_a - m_a^{\text{eq}})$
5. Restore to velocity space with $\mathbf{f}^{*} = M^{-1}\mathbf{m}$
6. Stream: move $f_i^{*}$ to the neighbor along $\mathbf{c}_i$
7. Apply boundary conditions and return to step 1

Only collision happens in moment space; streaming proceeds in velocity space as usual. The extra cost is just two $9\times 9$ matrix products per cell.

Python — Measuring Viscous Decay with a Taylor–Green Vortex#

For verification we use a Taylor–Green vortex. This flow has an exact viscous-decay solution, so you can check the implementation against numbers. The kinetic energy should decay as $\exp(-4\nu k^2 t)$, with $k$ the vortex wavenumber.

import numpy as np
 
# D2Q9 lattice
cx = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
cy = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
 
# Lallemand-Luo transformation matrix (rows: rho,e,eps,jx,qx,jy,qy,pxx,pxy)
M = np.array([
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [-4,-1,-1,-1,-1, 2, 2, 2, 2],
    [ 4,-2,-2,-2,-2, 1, 1, 1, 1],
    [ 0, 1, 0,-1, 0, 1,-1,-1, 1],
    [ 0,-2, 0, 2, 0, 1,-1,-1, 1],
    [ 0, 0, 1, 0,-1, 1, 1,-1,-1],
    [ 0, 0,-2, 0, 2, 1, 1,-1,-1],
    [ 0, 1,-1, 1,-1, 0, 0, 0, 0],
    [ 0, 0, 0, 0, 0, 1,-1, 1,-1],
], dtype=float)
Minv = np.linalg.inv(M)
 
def d2q9_equilibrium(rho, ux, uy):
    eq = np.zeros((9,) + rho.shape)
    for i in range(9):
        cu = cx[i] * ux + cy[i] * uy
        eq[i] = w[i] * rho * (1 + 3*cu + 4.5*cu**2 - 1.5*(ux**2 + uy**2))
    return eq
 
def relaxation_rates(nu):
    s_nu = 1.0 / (3*nu + 0.5)          # the rate that sets viscosity
    # conserved (rho,jx,jy)=0, high-order modes fixed in the stable band
    return np.array([0, 1.64, 1.54, 0, 1.7, 0, 1.7, s_nu, s_nu])
 
def taylor_green_init(N, U0):
    x = np.arange(N)
    X, Y = np.meshgrid(x, x, indexing='ij')
    k = 2*np.pi / N
    ux = -U0 * np.cos(k*X) * np.sin(k*Y)
    uy =  U0 * np.sin(k*X) * np.cos(k*Y)
    rho = 1 - 0.75 * U0**2 * (np.cos(2*k*X) + np.cos(2*k*Y))
    return d2q9_equilibrium(rho, ux, uy), k
 
def macros(f):
    rho = f.sum(axis=0)
    ux = (cx[:, None, None] * f).sum(axis=0) / rho
    uy = (cy[:, None, None] * f).sum(axis=0) / rho
    return rho, ux, uy
 
def mrt_collide(f, s):
    rho, ux, uy = macros(f)
    jx, jy = rho*ux, rho*uy
    m = np.tensordot(M, f, axes=([1], [0]))        # into moment space
    meq = np.zeros_like(m)
    meq[0] = rho
    meq[1] = -2*rho + 3*(jx**2 + jy**2)
    meq[2] = rho - 3*(jx**2 + jy**2)
    meq[3], meq[4] = jx, -jx
    meq[5], meq[6] = jy, -jy
    meq[7] = jx**2 - jy**2
    meq[8] = jx*jy
    m -= s[:, None, None] * (m - meq)               # per-moment relaxation
    return np.tensordot(Minv, m, axes=([1], [0]))    # back to velocity space
 
def stream(f):
    for i in range(9):
        f[i] = np.roll(f[i], (cx[i], cy[i]), axis=(0, 1))
    return f
 
def run_mrt_lbm(N=64, nu=0.004, U0=0.04, steps=2000):
    f, k = taylor_green_init(N, U0)
    s = relaxation_rates(nu)
    e0 = None
    for t in range(steps):
        f = mrt_collide(f, s)
        f = stream(f)
        if t % 400 == 0:
            _, ux, uy = macros(f)
            E = np.mean(ux**2 + uy**2)
            e0 = E if e0 is None else e0
            print(f"t={t:5d}  KE/KE0={E/e0:.4f}  analytic={np.exp(-4*nu*k**2*t):.4f}")
 
run_mrt_lbm()

The analytic value printed beside the kinetic energy is $\exp(-4\nu k^2 t)$. Defaults are $N=64$, $\nu=0.004$, $k=2\pi/N$, so $4\nu k^2 = 0.016\,(\pi/32)^2 \approx 1.54\times 10^{-4}$. At $t=400$ the exponent is about $-0.0617$ and the factor is about 0.940. At $t=1600$ it is about $-0.247$ and the factor is about 0.781. range(2000) means the last printed $t$ is 1600. Matching to two decimals is that factor against the numerical kinetic energy. BGK on the same lattice and the same $\nu$ diverges first, at a lower viscosity.

BGK and MRT Side by Side at the Same Re#

The viscosity slider slows the decay. The BGK/MRT button decides whether the lattice survives at that same $\nu$.

Drop the viscosity to around $0.002$ and leave BGK on. A rough pattern creeps across the lattice and "diverged" appears. Switch to MRT at the same viscosity and reset. The vortex decays smoothly. At $\nu=0.002$, $s_\nu = 1/(3\times 0.002+0.5) = 1/0.506 \approx 1.976$. BGK pushes the ghost modes to that same number. The listing's $1.64$, $1.54$, and $1.7$ stay beside it.

Leave the ghost rates when s_ν goes to 2#

MRT only separates the viscosity handle $s_\nu$ from the stability handles $s_e$ and $s_q$. Copying $s_\nu=1.953$ at $\nu=0.004$ onto every mode is BGK. Keeping $1.64$, $1.54$, and $1.7$ while $s_\nu$ alone climbs is MRT. The extra cost is two $9\times 9$ products per cell.

Without that split, a low-viscosity Taylor–Green dies as a checkerboard. Another way to stop the same runaway is an H-theorem inside the collision. The decay check is still the one line $\exp(-4\nu k^2 t)$.