At tau = 0.51 I read the viscosity 51x too high — the half step left behind by the discrete LBE

The continuum viscosity and the lattice viscosity differed by a factor of 51#

Look only at the continuous BGK equation and the relaxation time $\tau$ ties to the kinematic viscosity $\nu$ as $\nu = c_s^2 \tau$. Set $\tau = 0.51$ and $\nu$ should be 0.17. Feed the same $\tau$ to a D2Q9 code, measure the decay rate of a Taylor–Green vortex, and you get 0.00333. A factor of 51.

That is not a bug, and not a sign error. It is the relaxation time short by exactly $\Delta t/2$. The half step comes down two paths: a second-order Taylor expansion along a characteristic, and the midpoint rule on the collision integral. The two paths meet at the same place. Push the $\tau$ slider toward 0.51 and the orange dashes for $c_s^2\tau$ peel away from the lattice amplitude.

The blue curve is the vortex amplitude alive on the 32×32 lattice. The green dashes are $c_s^2(\tau - 1/2)$. At $\tau = 0.51$ green sits at 0.00333 and orange at 0.17.

The one line on the lattice is not a differential equation#

What we actually run is this single line.

Here $f_i$ is the distribution along the discrete velocity $\mathbf{c}_i$, $\Omega_i = f_i - f_i^{\rm eq}$ is the non-equilibrium part (how far the population sits from equilibrium), and $\Delta t$ is one step.

The left side is not a derivative. It is a difference taken across a finite gap $\Delta t$ that the code really jumps. Turning it into derivatives requires an expansion, and that expansion has a second-order term. The common move in CFD is to keep the first-order term and file the rest under "truncation error." Here the second-order term is not an error — it changes a physical coefficient.

Counting orders shows why it cannot be dropped. The second-order term scales as $\Delta t^2 D^2 f / 2$, while the collision term on the right carries $\Delta t / \tau$. In lattice units $\tau$ is an $O(1)$ number between 0.5 and 2, and $\Delta t = 1$. So the ratio of second to first order never shrinks. Freezing at $\Delta t = 1$ what should have gone to zero is the starting point of lattice gas methods, and this term is the price.

What survives when you expand to second order along a characteristic#

Write $D_i \equiv \partial_t + \mathbf{c}_i \cdot \nabla$ and the left side becomes an exponential of $D_i$. Keep terms to second order and divide by $\Delta t$.

The $D_i^2$ term is the one in question. Drop it and you get $\nu = c_s^2\tau$; keep it and you get something else.

Run a Chapman–Enskog expansion. With $f_i = f_i^{(0)} + \varepsilon f_i^{(1)} + \varepsilon^2 f_i^{(2)}$, $\partial_t = \varepsilon \partial_{t_1} + \varepsilon^2 \partial_{t_2}$, $\nabla = \varepsilon \nabla_1$, and $f_i^{(0)} = f_i^{\rm eq}$, the order-$\varepsilon$ balance reads

Reusing that relation inside the order-$\varepsilon^2$ balance is the whole trick. Since $D_{i1}^2 f_i^{(0)} = -\tfrac{1}{\tau} D_{i1} f_i^{(1)}$,

A bracket has appeared in front of $D_{i1} f_i^{(1)}$. The viscous stress comes from the second moment of $f_i^{(1)}$ multiplied by $\tau$, so that $\tau$ becomes $\tau\left(1 - \frac{\Delta t}{2\tau}\right) = \tau - \frac{\Delta t}{2}$.

In lattice units with $\Delta t = 1$ that is $c_s^2(\tau - 1/2)$. With $c_s^2 = 1/3$, a $\tau$ of 0.51 gives $\nu = 0.00333$. The green dashes above are exactly this number.

Pushing the leftover term inside the distribution#

There is a second route to the same answer that needs no expansion. Integrate the continuous BGK equation along a characteristic, but use the trapezoidal rule. Call the relaxation time $\lambda$.

The right side holds the future time level, so this is implicit. As written, it will not run. So change variables.

Rewritten in $\bar f_i$, the future-time terms cancel and vanish. Substituting $\bar\Omega_i = \bar f_i - f_i^{\rm eq} = \left(1 + \frac{\Delta t}{2\lambda}\right)\Omega_i$ leaves this.

Fully explicit, and character for character the line we ship — but only when $\tau = \lambda + \Delta t/2$. So the $\tau$ in the code is not the physical relaxation time. The physical one is $\lambda = \tau - \Delta t/2$. This is not a rounding issue. It is the bill for a change of variables. The trapezoidal integral that leaves this half step continues in the $\Delta t/2$ stuck on the relaxation time.

Reading the same line as a midpoint rule#

The trapezoidal rule is the $\theta = 1/2$ case of a $\theta$-weighted rule. Sweeping $\theta$ shows why exactly one half. The one-step amplification is $A_\theta(z) = \frac{1 - (1-\theta)z}{1 + \theta z}$ against an exact $e^{-z}$, where $z = \Delta t/\lambda$.

Step $\theta$ off 1/2 and the one-step amplification error goes as $(\theta - \tfrac12)z^2$, riding slope 2. Only at $\theta = 1/2$ does it drop to slope 3. The half is not a tuning knob. It is the unique value that meshes with the $1/2$ in the second-order coefficient of $e^{-z}$.

I decayed a Taylor–Green vortex at six values of tau#

That is as far as argument goes. Now ask the lattice directly. A 2D Taylor–Green vortex is a ruler for $\nu$, because its amplitude falls as $e^{-2\nu k^2 t}$.

import numpy as np
 
C = np.array([[0,0],[1,0],[0,1],[-1,0],[0,-1],[1,1],[-1,1],[-1,-1],[1,-1]])
W = np.array([4/9] + [1/9]*4 + [1/36]*4)
CS2 = 1/3
 
def equilibrium_d2q9(rho, ux, uy):
    cu = C[:,0,None,None]*ux + C[:,1,None,None]*uy
    u2 = ux*ux + uy*uy
    return W[:,None,None]*rho*(1 + cu/CS2 + cu*cu/(2*CS2**2) - u2/(2*CS2))
 
def taylor_green_init(n, u0):
    k = 2*np.pi/n
    x, y = np.meshgrid(np.arange(n), np.arange(n), indexing='ij')
    ux = -u0*np.cos(k*x)*np.sin(k*y)
    uy =  u0*np.sin(k*x)*np.cos(k*y)
    return np.ones((n,n)), ux, uy, k
 
def advance_lattice(f, tau):
    rho = f.sum(axis=0)
    ux = (C[:,0,None,None]*f).sum(axis=0)/rho
    uy = (C[:,1,None,None]*f).sum(axis=0)/rho
    f = f - (f - equilibrium_d2q9(rho, ux, uy))/tau      # collide
    for i in range(9):                                    # advect
        f[i] = np.roll(f[i], (C[i,0], C[i,1]), axis=(0,1))
    return f, ux, uy
 
def measure_decay_rate(n, tau, u0=0.01, nstep=400, skip=100):
    rho, ux, uy, k = taylor_green_init(n, u0)
    f = equilibrium_d2q9(rho, ux, uy)
    amp = []
    for s in range(nstep):
        f, ux, uy = advance_lattice(f, tau)
        amp.append(np.sqrt((ux**2 + uy**2).mean()))
    t = np.arange(skip, nstep, dtype=float)
    a = np.log(np.array(amp[skip:]))                      # drop the initial transient
    slope = np.polyfit(t, a, 1)[0]
    return -slope/(2*k*k), k
 
def viscosity_scan(n=32):
    print(f"{'tau':>6} {'nu_measured':>12} {'cs2*(tau-1/2)':>14} {'cs2*tau':>10} "
          f"{'err_half[%]':>11} {'err_naive[%]':>12}")
    for tau in (0.51, 0.55, 0.6, 0.8, 1.0, 1.5):
        nu, k = measure_decay_rate(n, tau)
        half = CS2*(tau - 0.5)
        naive = CS2*tau
        print(f"{tau:6.2f} {nu:12.6f} {half:14.6f} {naive:10.6f} "
              f"{100*abs(nu-half)/half:11.2f} {100*abs(nu-naive)/naive:12.2f}")
 
viscosity_scan()

   tau  nu_measured  cs2*(tau-1/2)    cs2*tau err_half[%] err_naive[%]
  0.51     0.003333       0.003333   0.170000        0.01        98.04
  0.55     0.016667       0.016667   0.183333        0.00        90.91
  0.60     0.033337       0.033333   0.200000        0.01        83.33
  0.80     0.100115       0.100000   0.266667        0.12        62.46
  1.00     0.167205       0.166667   0.333333        0.32        49.84
  1.50     0.337688       0.333333   0.500000        1.31        32.46

Pull theta off 1/2 and the slope on the right sticks at 2. Only at 1/2 does it fall to 3.

The blue region is the area under $e^{-s/\lambda}$, and the pink dashes are the single height the $\theta$ rule substitutes. Only $\theta = 1/2$ meshes with the $1/2$ in the second-order coefficient of $e^{-z}$. The 0.01% half-step error in the table is the residual when that midpoint is right.

The $\tau - 1/2$ column is off by 0.01 to 1.31 percent, the level of high-order error the lattice discretization leaves behind. The $\tau$ column is off by 32 to 98 percent. That is not a magnitude you can call an error.

The initial amplitude sits at 0.01 to erase compressibility. At a Mach number around 0.02, the truncated third-order term of $f^{\rm eq}$ does not leak into the decay rate. Discarding the first 100 steps serves the same purpose. Initializing with the equilibrium distribution leaves $f^{(1)}$ at zero, so the non-equilibrium part needs a few dozen steps to settle. Fit the slope through that window and the viscosity comes out too small.

The slow growth of the $\tau - 1/2$ error with $\tau$ is worth reading too. At $\tau = 1.5$ it reaches 1.31 percent, because the $O(\Delta t^2)$ terms dropped in the expansion grow with $\tau$. The half-step correction is not the final word — it is the story to second order.

Read the table as $\tau/(\tau - 1/2)$ and the factor of 51 is one line. At $\tau = 0.51$, $0.51/0.01 = 51$, and $0.170000/0.003333$ is 51.0 as well. At $\tau = 0.55$, $0.55/0.05 = 11$ and $0.183333/0.016667$ is 11. Then $\tau = 0.60$ gives 6, $\tau = 1.0$ gives 2, and $\tau = 1.5$ gives 1.5. The naive errors 98.04%, 90.91%, 83.33%, 62.46%, 49.84%, and 32.46% are that ratio walking toward 1.

The half-step column is much quieter. At $\tau = 0.51$ both measured and predicted values are 0.003333, and the error is 0.01%. At $\tau = 0.80$ the measurement is 0.100115 against 0.100000, a difference of 0.000115, or 0.12%. At $\tau = 1.50$ it is 0.337688 against 0.333333, a difference of 0.004355, or 1.31%. The dropped $O(\Delta t^2)$ growing with $\tau$ is those three cells. Amplitude 0.01, the first 100 steps discarded, and a 32-node grid were chosen to peel compressibility off that residual. At Mach 0.02 the truncated third-order term does not leak into the decay rate. The naive error 90.91% at $\tau = 0.55$ equals $1 - 0.016667/0.183333$. That ratio is 1/11, so reading $\nu = c_s^2\tau$ writes an eleven-fold viscosity into the code. Even at $\tau = 1.5$ the ratio is 1.5, and 32.46% is not a small error.

The boundary and the force term bill the same half step again#

This half step does not show up only once.

It is why bounce-back puts the wall halfway between two nodes rather than on a node: the reflection happens mid-step. As $\tau$ moves away from 0.5, that effective wall position starts drifting with $\tau$, which is the well-known loss of second-order accuracy.

It is also why the Guo-type force term carries a factor $\left(1 - \frac{1}{2\tau}\right)$. The force is integrated along the same characteristic under the same trapezoidal rule, and since the $\bar f$ transform acts on $\Omega_i$ alone, a coefficient is left over. Drop that bracket and a fluid that should sit still drifts slowly, or a hydrostatic profile under gravity tilts a little.

It is why halving the lattice spacing does not mean doubling $\tau$. Holding $\nu$ fixed requires solving $\tau_{\rm fine} - 1/2 = 2(\tau_{\rm coarse} - 1/2)$. Where that single line sits inside a refinement loop is what splits the error.

When $\tau$ varies node by node, this half step is the floor. Since $\nu \to 0$ as $\tau \to 0.5$, any problem where the viscosity drops locally toward zero runs into a wall built into the lattice itself.

0.5 is not a lower bound, it is the origin#

Memorizing $\tau > 0.5$ as a stability condition misses the cause. The physical relaxation time is $\tau - 1/2$. At $\tau = 0.51$, $c_s^2\tau = 0.170000$ and $c_s^2(\tau - 1/2) = 0.003333$, so the ratio is 51. A negative relaxation time is not physics.

What belongs in the code is $\nu$, not $\tau$. Put $\tau = 1/2 + \nu/c_s^2$ at one input and leave it there. The factor $\left(1 - \frac{1}{2\tau}\right)$ on the force is the same half-step bill. The limit of keeping that bracket next to the third moment of the force continues in the line a forcing term has to respect. If the decay slows by a factor of 51, somewhere $\tau$ is being read as $\nu$.