Pushing the Frame to 0.2 Cut the Viscosity by 5.97% — The Missing Third Moment in LBM Equilibria

Should still water and moving water have the same viscosity?#

Measure the same fluid twice. Once in a box at rest, once while the box slides along at constant speed. The viscosity must come out identical both times. Galilean invariance — the rule that physics reads the same in two frames moving at constant relative velocity — demands it.

Run that experiment on a D2Q9 lattice Boltzmann solver and the second number is 5.97% smaller. Refining the grid does not help. Changing the relaxation time $\tau$ barely moves the ratio. This post traces that 5.97% down to a single moment of the equilibrium distribution. The short answer: this is not a stability problem, it is an algebra problem. The lattice velocity set simply cannot build the term that is missing.

Taylor-expanding the one line the lattice solves#

LBM advances one equation.

Here $f_i$ is the population riding the lattice velocity $\mathbf{c}_i$, $f_i^{eq}$ is the local equilibrium, and $\tau$ is the relaxation time.

Expand the left side in $\Delta t$ and a second-order term survives.

That second-order term is what produces the famous $\tau - 1/2$ in the Chapman–Enskog expansion, the multiscale technique that splits the distribution into a power series in the Knudsen number. The next step of the same expansion hands you a requirement: to recover Navier–Stokes, the equilibrium must match four moments exactly.

The first three fix mass, momentum, and pressure. The fourth, the third-order moment, is what replaces the time derivative of the viscous stress. Get it wrong and continuity stays clean, the Euler level stays clean, and only the viscous term is contaminated.

Push the velocity slider up from zero and the amber Maxwellian curve splits off the blue lattice line.

sweep u shows the gap tracking an odd power of $u$. Turn on speeds ±2 and the bars grow to five while the gap closes to zero.

$c_x^3 = c_x$, so the term cannot exist#

The $x$ components of D2Q9 take only three values: $-1$, $0$, $1$. Every one of them equals its own cube. So for any set of populations $f_i$,

holds identically. It does not matter how the equilibrium is designed. Fill the nine slots with random numbers and the identity still holds. The third moment is already nailed to the first.

What a Maxwellian asks for is $\rho(u_x^3 + 3 c_s^2 u_x)$, and on D2Q9 $c_s^2 = 1/3$, so $3c_s^2 = 1$ and the target collapses to $\rho(u_x + u_x^3)$. The difference between the two is exactly $\rho u_x^3$. Not approximately that size — that term, with no remainder.

What the defect leaves behind in the momentum equation#

The missing moment lands directly on the viscous stress. Carry the Chapman–Enskog expansion through and one extra term appears.

Now let a mean flow $U$ run along $x$ with a small disturbance $u'$ on top. Keep only the part linear in $u'$ and the term becomes $-(\tau - 1/2)\Delta t\, U^2 \partial_x^2 (\rho u'_\alpha)$, which has exactly the shape of a viscous term. The $x$-direction viscosity is rewritten wholesale.

The relative error is $-U^2 / c_s^2$. Notice that $\tau$ cancels between numerator and denominator: no amount of tuning the relaxation changes the percentage. And the error scales with the square of the velocity, which in Mach number is simply $-\mathrm{Ma}^2$.

Put the arithmetic on the page. With $c_s^2 = 1/3$ and $U = 0.2$, $U^2 = 0.04$ and $U^2/c_s^2 = 0.12$, so the $x$ viscosity alone is off by 12%. Taylor–Green has $k_x = k_y$, and half of that error is 6%. The row $U_0 = 0.20$ measures $\nu = 0.0940296$ against $0.100000$, which is $-5.9704\%$. 6.00 and 5.97 differ by 0.03 percentage points.

At $\tau = 0.6$ theory gives $c_s^2 \times 0.1 = 0.0333333$. At $\tau = 1.0$ it gives $c_s^2 \times 0.5 = 0.1666667$, five times larger. The relative error at $U_0 = 0.15$ only moves from $-3.4076\%$ to $-3.3099\%$. Substituting $U^2 = 0.0225$ into $-U^2/(2c_s^2)$ yields $-0.03375$, or $-3.375\%$. Those three rows hug that single number. The algebra that produces $\tau - 1/2$ is in the second-order expansion of the discrete LBE. In this error the factor cancels.

The audit column already holds the defect. At $u = 0.20$ the lattice sum is $0.20000000$, the Maxwellian value is $0.20800000$, and the gap is $0.00800000$. Because $0.2^3 = 0.008$, the ratio is 1.00000. At $u = 0.05$ the gap is $0.00012500$ and the ratio is still 1. At $u = 0.30$ the gap $0.02700000$ equals $0.3^3$. Random populations leave a difference of $0.000\mathrm{e}{+}00$. Same term, two faces.

Measuring one vortex from two frames in Python#

Put a Taylor–Green vortex on a 64×64 periodic lattice and add a uniform velocity $U_0$ to it. Fit the log slope of the decaying amplitude and the viscosity falls out. Only $U_0$ changes; everything else is held fixed.

import numpy as np
 
CX  = np.array([0, 1, 0, -1,  0, 1, -1, -1,  1], dtype=float)
CY  = np.array([0, 0, 1,  0, -1, 1,  1, -1, -1], dtype=float)
W   = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
CS2 = 1.0 / 3.0
 
def f_equilibrium(rho, ux, uy):
    # standard equilibrium, expanded to second order
    cu = CX[:, None, None] * ux + CY[:, None, None] * uy
    usq = ux * ux + uy * uy
    return W[:, None, None] * rho * (1 + cu / CS2 + cu * cu / (2 * CS2**2) - usq / (2 * CS2))
 
print("[1] third moment audit:  sum f_i c_ix^3   vs   Maxwell rho*(u^3 + 3*cs2*u)")
print("   u        lattice      Maxwell        gap        gap/u^3")
for u in (0.05, 0.10, 0.20, 0.30):
    feq = f_equilibrium(np.ones((1, 1)), np.full((1, 1), u), np.zeros((1, 1)))
    lat = float((feq * (CX**3)[:, None, None]).sum())
    exact = u**3 + 3 * CS2 * u
    print(f"  {u:4.2f}  {lat:11.8f}  {exact:11.8f}  {exact-lat:11.8f}   {(exact-lat)/u**3:8.5f}")
 
rng = np.random.default_rng(7)
frand = rng.random((9, 1, 1))
d = float((frand * (CX**3)[:, None, None]).sum() - (frand * CX[:, None, None]).sum())
print(f"  any random f:  sum f c_x^3 - sum f c_x = {d:.3e}   (c_x^3 = c_x, so always 0)")
 
def taylor_green_run(U0, tau=0.8, N=64, steps=900, amp=0.04, warmup=300, sample=25):
    # same vortex, ridden by a uniform U0; recover viscosity from the decay rate
    k = 2 * np.pi / N
    x = np.arange(N)[:, None] * np.ones(N)[None, :]
    y = np.ones(N)[:, None] * np.arange(N)[None, :]
    ux = U0 - amp * np.cos(k * x) * np.sin(k * y)
    uy = amp * np.sin(k * x) * np.cos(k * y)
    rho = np.ones((N, N))
    f = f_equilibrium(rho, ux, uy)
    ts, logs = [], []
    for n in range(steps + 1):
        rho = f.sum(axis=0)
        ux = (f * CX[:, None, None]).sum(axis=0) / rho
        uy = (f * CY[:, None, None]).sum(axis=0) / rho
        if n % sample == 0:
            up, vp = ux - ux.mean(), uy - uy.mean()   # fluctuation left after removing the mean flow
            ts.append(n)
            logs.append(0.5 * np.log(2 * np.mean(up * up + vp * vp)))
        f += (f_equilibrium(rho, ux, uy) - f) / tau   # BGK collision
        for i in range(9):                            # streaming
            f[i] = np.roll(np.roll(f[i], int(CX[i]), axis=0), int(CY[i]), axis=1)
    ts, logs = np.array(ts, dtype=float), np.array(logs)
    m = ts >= warmup
    slope = np.polyfit(ts[m], logs[m], 1)[0]
    return -slope / (2 * k * k)
 
print("\n[2] same vortex measured in a uniformly moving frame  (tau=0.8, nu_theory=0.100000)")
nu_th = CS2 * (0.8 - 0.5)
print("   U0      nu_eff      rel.err     rel.err/U0^2")
for U0 in (0.00, 0.05, 0.10, 0.15, 0.20):
    nu = taylor_green_run(U0)
    e = nu / nu_th - 1
    tail = f"{e/U0**2:10.4f}" if U0 > 0 else "         -"
    print(f"  {U0:4.2f}  {nu:10.7f}  {e:+9.4%}  {tail}")
 
print("\n[3] does the error depend on tau?  (U0=0.15 fixed)")
print("   tau     nu_theory    nu_eff      rel.err")
for tau in (0.6, 0.8, 1.0):
    nu = taylor_green_run(0.15, tau=tau)
    th = CS2 * (tau - 0.5)
    print(f"  {tau:4.2f}  {th:10.7f}  {nu:10.7f}  {nu/th-1:+9.4%}")

[1] third moment audit:  sum f_i c_ix^3   vs   Maxwell rho*(u^3 + 3*cs2*u)
   u        lattice      Maxwell        gap        gap/u^3
  0.05   0.05000000   0.05012500   0.00012500    1.00000
  0.10   0.10000000   0.10100000   0.00100000    1.00000
  0.20   0.20000000   0.20800000   0.00800000    1.00000
  0.30   0.30000000   0.32700000   0.02700000    1.00000
  any random f:  sum f c_x^3 - sum f c_x = 0.000e+00   (c_x^3 = c_x, so always 0)
 
[2] same vortex measured in a uniformly moving frame  (tau=0.8, nu_theory=0.100000)
   U0      nu_eff      rel.err     rel.err/U0^2
  0.00   0.1000175   +0.0175%           -
  0.05   0.0996433   -0.3567%     -1.4268
  0.10   0.0985207   -1.4793%     -1.4793
  0.15   0.0966495   -3.3505%     -1.4891
  0.20   0.0940296   -5.9704%     -1.4926
 
[3] does the error depend on tau?  (U0=0.15 fixed)
   tau     nu_theory    nu_eff      rel.err
  0.60   0.0333333   0.0321975   -3.4076%
  0.80   0.1000000   0.0966495   -3.3505%
  1.00   0.1666667   0.1611502   -3.3099%

Raise frame velocity U and the right-hand vortex drifts while its decay slows.

Same fluid, and the right-hand vortex fades more slowly. Wiggle the $\tau$ slider and the error readout barely moves. That is the $-3.41\%$ to $-3.31\%$ band in the table.

The error rides on $U^2$, not on $\tau$#

The last column of [1] is 1.00000 all the way down. The claim that the defect is exactly $\rho u^3$ holds to five decimals.

In [2] the resting frame reports $+0.0175\%$, which is the noise floor of the measurement itself. Raise $U_0$ and the error walks to $-0.36\%$, $-1.48\%$, $-3.35\%$, $-5.97\%$. The last column divides that error by $U_0^2$, and as $U_0$ shrinks it settles on $-1.5$.

That $-1.5$ is the prediction from the previous section. The relative error $-U^2/c_s^2$ hits the $x$ component only, and a Taylor–Green mode has $k_x = k_y$, so its decay rate averages the two directions. Half of $-U^2/c_s^2$ is $-U^2/(2c_s^2) = -1.5\,U^2$. At $U_0 = 0.2$ the prediction is $-6.00\%$ against a measured $-5.97\%$.

Result [3] matters more. Raising $\tau$ from 0.6 to 1.0 multiplies the viscosity by five, and the relative error moves only from $-3.41\%$ to $-3.31\%$. You cannot bury this error under more viscosity.

Three ways to fix it, and what each one costs#

Slow down. The error goes as $U^2$, so halving the lattice velocity cuts it to a quarter. The bill is step count: reproducing the same physical time now takes more iterations. The LBM folklore rule of "keep the Mach number below 0.1" is usually sold as a stability rule, but this defect is what actually binds first.

Add a correction term. Compute $\partial_\gamma(\rho u_\alpha u_\beta u_\gamma)$ by finite differences and inject it into the collision with the opposite sign. The cost is one gradient evaluation plus a partial loss of the strictly local collision that makes LBM attractive. For rotating frames and other large-mean-flow problems, it is still the cheapest option on the table.

Widen the velocity set. A multi-speed lattice that includes $\pm 2$ has $c^3 \neq c$, so the third moment can be matched exactly. The speeds ±2 button in the first visualization is doing precisely that solve. You pay with a wider stencil, more memory, and messier boundary treatment — and some weights can go negative, which puts stability back on the table.

Where the defect actually shows up#

Anywhere the mean flow is large. Rotating frames in turbomachinery, sliding meshes, turbulence riding on a fast uniform stream, a frame attached to a moving body. A solver that looks flawless on static benchmarks and then reports the wrong effective Reynolds number on exactly these problems is showing this symptom.

Resetting $\tau$ on every grid level leaves this defect in place, because the relative error does not depend on $\tau$. A quantity nobody typed in as a material property is fixed by the lattice itself.

Diagnosis is cheap. Rerun the same case with a mean flow added and compare the decay rate or the drag coefficient. If the gap grows like $U^2$, this is the term. The boundary treatment is the wrong first suspect.

Lower $u/c_s$, or widen the speeds to $\pm 2$#

D2Q9 drops $\rho u^3$ at third order. At $U_0 = 0.2$ the bill is $-5.97\%$. Refining the grid, or lifting $\tau$ from 0.6 to 1.0, barely moves that ratio.

Three cuts actually shrink it. Lower $u/c_s$ and the error falls with $U^2$. Put $\partial_\gamma(\rho u_\alpha u_\beta u_\gamma)$ back into the collision. Or widen the velocity set. Why three speeds lock into $c^3 = c$ continues in the Hermite roots that pick the lattice velocities. Leave that choice alone and a frame at 0.2 cuts the viscosity by 5.97% again.