Where Half the Force Goes in LBM — Four Forcing Schemes and 1−1/(2τ)

In 1993 Shan and Chen pulled water and oil apart on a lattice. The prescription was modest. All they did was add $\tau\mathbf{F}/\rho$ to the velocity fed to the equilibrium during collision. In 1998 He attached the force directly to the equilibrium distribution. In 2002 Guo argued that the earlier prescriptions leave an error in the stress, and put $(1-\frac{1}{2\tau})$ in front of the force term. In 2004 Kupershtokh did the same job with nothing but the difference of two equilibrium distributions.

Four answers to one question. So which one should you use? Today I put that question straight to a D2Q9 channel. The answer was not what I expected. In steady flow, every scheme gave the same error out to the fourth decimal place. The real trap is not where you pick a scheme. The derivation locates it. The numbers size it.

The moment dial changes scheme, $\tau$, $|u|$, and the prefactor together, and the M0, M1, and M2 bars move with them.

Spin the dial and M0 stays at zero. Raise $|u|$ and only plain's M2 misses the target. Kill the prefactor and M1 and $\Delta j$ break in the same step.

Cutting the force term out of the continuous equation#

The Boltzmann equation with an external force gains one term.

$f$ is the distribution function, $\boldsymbol{\xi}$ the particle velocity, $\mathbf{F}$ the body force per unit volume, and $\lambda$ the relaxation time.

The trouble is the third term. The velocity-space derivative $\nabla_{\boldsymbol{\xi}} f$ does not exist on a lattice. All we have is nine discrete velocities. So we swap $f$ for $f^{\rm eq}$ and handle the derivative analytically.

$c_s$ is the lattice sound speed ($c_s^2=1/3$ in D2Q9) and $\mathbf{u}$ is the fluid velocity.

Is it legal to swap $f$ for $f^{\rm eq}$? At low Mach number, yes. The non-equilibrium part $f^{(1)}$ is first order in the Knudsen number (mean free path over characteristic length), so once it is multiplied by the force it drops to a negligible order. What the swap does instead is leave a residual error in the stress. That is exactly where the schemes part ways.

Project this expression onto the discrete velocities $\mathbf{c}_i$ and truncate the Hermite series at second order. Out comes the force term on the lattice.

$w_i$ is the lattice weight. The reason for cutting at second order is plain. The moments needed to recover Navier–Stokes go up to second order.

Three moments are worth pinning down.

Zeroth order is mass. A force creates no mass, so it has to vanish. First order is momentum, second order is stress.

The half left behind by the trapezoidal rule#

Now integrate the discrete-velocity Boltzmann equation along a characteristic over $\delta t$. Integrating the collision and force terms with the trapezoidal rule buys second-order accuracy. The price is that $t+\delta t$ slips onto the right-hand side.

The superscripts $n$ and $n+1$ mean $t$ and $t+\delta t$. As written this is implicit. It would take a linear system solved every step.

LBM erases that implicitness with a change of variables.

Rewrite in terms of $\bar f_i$ and set $\tau = \lambda/\delta t + 1/2$, and you get the familiar one-liner.

The transformation left two traces. One is the $(1-\frac{1}{2\tau})$ in front of the force. The other is easier to miss. What we store in the array is $\bar f_i$, not $f_i$. So the momentum is off by that same half.

These two are one body. They fell out of the same transformation together. The half step the trapezoidal rule leaves in the viscosity comes from that same root. Keep one and drop the other and the books do not balance. How badly they fail to balance gets measured below.

A comparison table for four schemes#

Call $\mathbf{j}=\sum_i \bar f_i\mathbf{c}_i$ the raw momentum of the stored distributions, and $\delta\mathbf{u}=\delta t\,\mathbf{F}/\rho$ the velocity increment the force delivers in one step.

Read the five rows down the columns and the common ground shows. Before the gain multiplies it, every source term's first moment is exactly $\mathbf{F}$. No surprise — they were designed that way. What differs is the second moment, which is the stress.

But the target the change of variables demands is not $\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u}$. It is $(1-\frac{1}{2\tau})(\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u})$. Only the bottom row carries that coefficient explicitly. The rest go along carrying a stress deviation of size $\frac{1}{2\tau}$.

The deviation scales as force times velocity, order $u F$. At low Mach, $u$ is small, so the term itself is small. The $\mathbf{F}\mathbf{F}$ term that Shan–Chen and EDM leave behind is smaller still, order $F^2$. That is why it barely matters which one you use when the driving is weak gravity. The differences surface where the force is strong (when $\tau\mathbf{F}/\rho$ rivals $\mathbf{u}$), or where the force changes sharply across an interface, as in multiphase flow.

The M1 target is not $\mathbf{F}$ but $(1-\frac{1}{2\tau})\mathbf{F}$, because the bar measures the source after the gain. The missing $\mathbf{F}/(2\tau)$ comes back from the half-shifted equilibrium, and $\Delta j$ lands on $\mathbf{F}$. Only M2 splits by scheme. Plain is the wrong shape, so no gain fixes it. Push $\tau$ toward 2 and the Shan–Chen deviation grows.

Steady Poiseuille cannot tell them apart#

Let's go to the most common validation problem. A channel pushed by a body force between an upper and a lower wall. The steady solution is a known parabola.

$H$ is the channel height and $\nu$ the kinematic viscosity. With half-way bounce-back the wall sits half a cell outside the lattice nodes, so $H$ equals the number of fluid nodes and node $j$ sits at $y=j+0.5$.

Swapping the scheme leaves the parabola's amplitude alone. Maximum relative error after 40,000 steps, 33 fluid nodes, $F_x = 10^{-5}$.

There is nothing to compare across the columns. The values are identical. The error that remains comes from the discretization of bounce-back, not from the scheme.

The reason lies in where the second moment enters. This flow is steady and unidirectional. The velocity has a single component $u_x(y)$ and the force has only an $x$ component. The term where the schemes part ways is $\mathbf{u}\mathbf{F}$, the $xx$ component. But what actually governs the momentum balance of this flow is the $xy$ shear stress. The $xx$ deviation the force creates has no route into the balance.

That $\partial_t\mathbf{u}=0$ works in the same direction. There is no room for the force-term error to accumulate in time. To see the difference between schemes you have to go to unsteady flow, or to a spatially varying force. With this test, at least, the five cannot be told apart.

This is not bad news. It is something to know when choosing a validation problem. It means a passing Poiseuille does not let you say your forcing implementation is right.

Count the half twice and the force grows#

So what does actually go wrong? The two things called one body earlier — the half-step velocity correction and the coefficient in front of the force.

Count the momentum actually injected in one step. Add $\delta t\mathbf{F}/2\rho$ to the equilibrium velocity and collision carries $1/\tau$ of it across each step. On top of that comes the share the force term delivers directly. Write $s\in\{0,1\}$ for whether the half correction is on, and $g$ for the coefficient in front of the force, and

This has to equal $\delta t F$. The condition is a single one, $g = 1 - s/(2\tau)$. The two switches are not independent.

Mismatch them and the amplitude misses by a predictable ratio. Here are the amplitude ratios measured with the Guo force term held fixed and only the two switches changed.

Measurement and prediction agree to the third decimal place. At $\tau=0.6$, turning on only the half correction makes the force 83% larger. Turning on only the coefficient shrinks it to a sixth.

This is not a subtle accuracy issue. It is exactly what happens when you copy Guo's $F_i$ straight out of a paper and leave the velocity as the old code's $\mathbf{j}/\rho$. Start touching $\tau$ because the viscosity looks off and you go in deeper. Every change to $\tau$ moves the error ratio along with it.

There is a way to spot it from the symptoms. If doubling the grid does not shrink the error ratio, it is not a discretization error. If the ratio converges toward 1.5 as $\tau$ is brought near 1 and grows as if diverging as $\tau$ drops toward 0.5, that is the fingerprint of $1+\frac{1}{2\tau}$. And if raising the force leaves the ratio unchanged, it is a bookkeeping problem, not a nonlinear one. When all three line up, do not suspect the scheme — open up the velocity definition first.

At $\tau=0.6$, $2\tau=1.2$, so $1/(2\tau)=5/6\approx 0.8333$. Both switches on means $s=1$ and $g=1-1/(2\tau)=1/6$. Substituting into $(1/\tau)\cdot(s/2)+g$ gives $(1/0.6)\cdot 0.5 + 1/6 = 5/6 + 1/6 = 1$. The both column, 0.9991, is that 1 with the lattice shaved off.

Shift only freezes $g$ at 1, so the slot becomes $5/6+1=1.8333$. The table prints 1.8316. The gap of 0.0017 is about 0.09% of the prediction. Gain only sets $s=0$ and $g=1/6$, so the injection is 0.1667. Measured 0.1666 meets it at the fourth decimal.

At $\tau=1$, $1/(2\tau)=1/2$. Shift only predicts 1.5 and gain only 0.5. The table says 1.5002 and 0.5005. Both on gives 1.0003, 0.03% off a target of 1. At $\tau=1.8$, $1/(2\tau)=1/3.6=0.2778$. Predictions 1.2778 and 0.7222 sit next to 1.2867 and 0.7280. $(1.2867-1.2778)/1.2778$ is about 0.007, or 0.7%. The scheme error at that $\tau$ is 0.7358%. Same size.

The relative-error columns read 0.0875% at $\tau=0.6$, 0.0306% at $\tau=1.0$, and 0.7358% at $\tau=1.8$. With $\nu=c_s^2(\tau-1/2)$ and $c_s^2=1/3$, the viscosities are $1/30$, $1/6$, and $13/30$. The $(\tau-1/2)$ ratios run 1 to 5 to 13. The peak is $u(H/2)=F_x H^2/(8\nu)$, so at $F_x=10^{-5}$ and $H=33$ the $\tau=1.8$ peak is $1/13$ of the $\tau=0.6$ peak. The percent grows there. The five cells still share one number.

Python — mismatching the two switches#

A D2Q9 channel solver with the two switches pulled out as arguments. The flow is homogeneous in $x$, so only one column is kept. The five schemes in the table above only need source_terms and the velocity fed to the equilibrium swapped out, so here we fix on the Guo force term alone.

import numpy as np
 
EX = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
EY = 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])
OPP = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6])
CS2 = 1.0 / 3.0
 
 
def lattice_equilibrium(rho, ux, uy):
    feq = np.empty((9, rho.size))
    usq = ux**2 + uy**2
    for i in range(9):
        eu = EX[i] * ux + EY[i] * uy
        feq[i] = W[i] * rho * (1 + eu/CS2 + eu**2/(2*CS2**2) - usq/(2*CS2))
    return feq
 
 
def source_terms(rho, ux, uy, fx):
    """Guo forcing term, before the (1 - 1/2tau) gain."""
    src = np.empty((9, rho.size))
    for i in range(9):
        eu = EX[i]*ux + EY[i]*uy
        src[i] = W[i] * ((EX[i] - ux)/CS2 + eu*EX[i]/CS2**2) * fx
    return src
 
 
def run_forced_channel(tau, half_shift, prefactor, ny=33, fx=1.0e-5, steps=40000):
    rho = np.ones(ny)
    f = lattice_equilibrium(rho, np.zeros(ny), np.zeros(ny))
    gain = (1.0 - 1.0/(2*tau)) if prefactor else 1.0
    for _ in range(steps):
        rho = f.sum(axis=0)
        jx = (f*EX[:, None]).sum(axis=0)
        jy = (f*EY[:, None]).sum(axis=0)
        ux = (jx + (0.5*fx if half_shift else 0.0)) / rho   # switch 1
        uy = jy / rho
        fpost = f - (f - lattice_equilibrium(rho, ux, uy))/tau \
                + gain*source_terms(rho, ux, uy, fx)        # switch 2
        for i in range(9):                                   # streaming in y
            f[i] = fpost[i] if EY[i] == 0 else np.roll(fpost[i], EY[i])
        for i in range(9):                                   # half-way bounce-back
            if EY[i] > 0:
                f[i, 0] = fpost[OPP[i], 0]
            elif EY[i] < 0:
                f[i, -1] = fpost[OPP[i], -1]
    jx = (f*EX[:, None]).sum(axis=0)
    return (jx + 0.5*fx) / f.sum(axis=0)      # physical velocity, always shifted
 
 
NY, FX = 33, 1.0e-5
y = np.arange(NY) + 0.5
for tau in (0.6, 1.0, 1.8):
    exact = FX / (2*CS2*(tau - 0.5)) * y * (NY - y)
    ok = run_forced_channel(tau, True, True).max() / exact.max()
    m1 = run_forced_channel(tau, True, False).max() / exact.max()
    m2 = run_forced_channel(tau, False, True).max() / exact.max()
    print(f"tau={tau:4.2f}  both={ok:.4f}  shift_only={m1:.4f} (pred {1+1/(2*tau):.4f})"
          f"  gain_only={m2:.4f} (pred {1-1/(2*tau):.4f})")

The output looks like this.

tau=0.60  both=0.9991  shift_only=1.8316 (pred 1.8333)  gain_only=0.1666 (pred 0.1667)
tau=1.00  both=1.0003  shift_only=1.5002 (pred 1.5000)  gain_only=0.5005 (pred 0.5000)
tau=1.80  both=1.0074  shift_only=1.2867 (pred 1.2778)  gain_only=0.7280 (pred 0.7222)

Match the channel lab's scheme buttons, $\tau$, half-shift, and prefactor to the both, shift_only, and gain_only columns printed above.

All four scheme buttons leave the teal curve on the dashed parabola. Kill one switch and, at $\tau=0.6$, the amplitude jumps toward 1.83 or 0.17.

At $\tau=1.8$ the measurement misses the prediction by 0.7%, and under the same conditions the scheme's own discretization error is 0.74%. Same size. Which means the mismatch comes from the lattice, not from the bookkeeping.

Leave only the half-shift on and the force grows to 1.83#

At $\tau=0.6$, shift_only reads 1.8316. The prediction $1+\frac{1}{2\tau}=1.8333$ agrees to the third decimal. Copy Guo's source and keep the velocity as $\mathbf{j}/\rho$, and the run lands in exactly that cell: shift on, prefactor off. Doubling the grid does not shrink the ratio. Bring $\tau$ to 1 and the ratio sits near 1.5. Drop $\tau$ toward 0.5 and the ratio grows. That is the fingerprint of $1+\frac{1}{2\tau}$.

Every scheme's first moment is $\mathbf{F}$. The split is the second moment, and the target carries $(1-\frac{1}{2\tau})$. A steady channel has no path for that $xx$ error into the shear balance, so the five columns share 0.0875%, 0.0306%, and 0.7358%. Check M2 when the flow is unsteady or the force jumps across an interface. Filling the third moment on top of the second is a different bill, told in the third-moment limit of the force term.