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.

Check that third moment yourself in the simulation below.

Drag the velocity slider and the amber curve — what a Maxwellian demands — pulls away from the blue line, which is what D2Q9 actually delivers. Hit sweep u and the gap tracks 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$.

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%

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.

Raise frame velocity U and the right-hand vortex drifts across the panel while fading more slowly than the left one. Same fluid, and yet the right one refuses to die. Wiggle the $\tau$ slider and watch the error readout below stay put.

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.

In the non-equilibrium rescaling of refined LBM grids the fix was to reset $\tau$ per level, and this defect survives that operation untouched, because its relative size does not depend on $\tau$. It also has the same shape as the two constants buried in thermal LBM: a quantity nobody entered as a material property, decided instead by the structure of the lattice.

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 your term. No need to go suspecting the boundary treatment, as in the parity test in STL voxelization.

Next time the viscosity depends on the frame#

Stability is not the only thing capping velocity in lattice Boltzmann. The moment order the equilibrium can reach caps it first. D2Q9 drops $\rho u^3$ at third order, and the invoice arrives as a viscosity error of $-\mathrm{Ma}^2$.

When numbers look off, refining the grid and retuning $\tau$ are the reflexes. Neither one touches this error. It shrinks only when $u/c_s$ comes down, when the missing term is put back by hand, or when the velocity set gets wider.