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

Two collision models relax toward the same equilibrium on the same lattice, yet one spins smoothly while the other shatters into a checkerboard. The difference is a single choice. BGK relaxes all nine distribution functions with one relaxation time, while MRT (Multiple Relaxation Time) maps them into moments and relaxes each at its own rate. This post traces why that choice decides stability, using the D2Q9 lattice Boltzmann method. We move distributions into moments with a transformation matrix $M$, tighten the degrees of freedom that have nothing to do with viscosity, and watch with a Taylor–Green vortex how MRT holds the low-viscosity regime where BGK blows up.

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#

The key idea is simple. Instead of relaxing the nine distribution functions directly, first map them into physically meaningful moments, 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 point is that the remaining rates $s_e$ (bulk), $s_\varepsilon$, and $s_q$ (energy flux) have no effect on viscosity at all. They are free parameters. Even when you push $s_\nu$ near $2$ for low viscosity, the high-order modes can be kept separately in the stable band of $1.5\sim1.8$. The two handles BGK had tied together come apart.

Try adjusting the relaxation rates directly in the bars below.

Push $s_\nu$ up to $1.99$ and the viscosity $\nu$ drops nearly to $0$, yet you can still keep $s_e, s_q$ in the stable band. Notice too that the conserved-moment bars ($\rho, j_x, j_y$) stay pinned at zero.

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 output shows the numerical decay matching the analytic curve $\exp(-4\nu k^2 t)$ to two decimal places. Running 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#

In the simulation below, lower the viscosity slider and switch the collision model with the BGK/MRT button.

Drop the viscosity to around $0.002$ and leave it on BGK, and a rough pattern creeps across the lattice as "diverged" appears. Switch to MRT at the same viscosity and reset, and the vortex decays smoothly. Even though both models target the same Reynolds number, their stability limits differ at a glance.

The Next Time You Touch Lattice Boltzmann#

MRT is not magic. It merely physically separates the handle that sets viscosity ($s_\nu$) from the handles that set stability ($s_e, s_q$). That separation widens the operating range of lattice Boltzmann at low viscosity and high Reynolds number.

Three things to keep. First, BGK blows up at low viscosity because the high-order modes unrelated to viscosity slow down along with it. Second, MRT moves only collision into moment space to tighten those modes separately, at a cost of two matrix products per cell. Third, you can always verify the implementation against the $\exp(-4\nu k^2 t)$ decay of a Taylor–Green vortex.