Filling the Third Moment Honestly Made the Error 2.8x Worse — What a Forcing Term Actually Owes

How many moments does a forcing term owe?#

When we put gravity or surface tension into an LBM code, we match three lines. Zeroth moment zero, first moment $\mathbf{F}$, second moment $\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u}$. What about the third? The force term of the continuous Boltzmann equation has a definite third moment too. Matching that value on the lattice looks like a free accuracy gain.

The numbers went the other way. Filling the third moment with what kinetic theory demands pushed the error from 0.053% to 0.150%. The sliders $u_x$ and $u_y$ are the velocity. Raising them splits the moment rows that had matched. The checkbox closes only the yellow rows.

Pull the velocity back toward zero and the rows line up again. The red row stays open no matter how far the sliders go.

What the continuous Boltzmann equation demands#

For a body force $\mathbf{F}$, the force term of the Boltzmann equation is a derivative in velocity space.

Here $\xi$ is molecular velocity and $\rho$ is density. One integration by parts closes the third moment.

with $M_{\alpha\beta}=\int f^{eq}\xi_\alpha\xi_\beta\,d\boldsymbol{\xi} = \rho(u_\alpha u_\beta + c_s^2\delta_{\alpha\beta})$. So the third moment being asked for reads as follows.

A $c_s^2$ piece and a $uuF$ piece sit side by side. The first is velocity-independent, the second is quadratic in velocity.

What Guo forcing actually supplies#

The forcing term on the lattice is usually written in Guo's form. Strip the half-step factor and call the bracket $\hat{S}_i$.

$w_i$ are the lattice weights and $\mathbf{c}_i$ the discrete velocities. The actual $S_i$ carries a factor $(1-1/2\tau)$ on top. Why that factor is there is covered separately in the half step the discrete LBE leaves behind.

Compute the third moment of this $\hat{S}_i$ using the fourth-order isotropy of D2Q9 and this is what survives.

The $c_s^2$ piece is exactly right. Not one of the $uuF$ terms is there. The second bracket builds a fifth moment, which vanishes on D2Q9.

That is the split the velocity sliders already opened: two yellow rows move, and the red row stays locked.

$c_x^3 = c_x$: the two components the lattice locked#

The red row is the lattice itself. D2Q9 velocity components are only $-1, 0, 1$. That makes $c_{ix}^3 = c_{ix}$ in every direction. The two vectors are the same vector.

$\Phi_{xxx}$ is not a value you get to choose. Fixing the first moment to $F_x$ fixes it at the same time. The demand is $3F_x(u_x^2+c_s^2) = F_x + 3F_x u_x^2$, so $3F_x u_x^2$ is a permanent deficit.

$\Phi_{xxy}$ and $\Phi_{xyy}$, by contrast, are independent. Of nine degrees of freedom the 0th to 2nd moments take six, leaving three. So of the four components, two are free and two are locked. You can fix half of a tensor.

A term that carries no mass and no momentum still moves momentum#

You can build a correction vector $v_i$ that changes only $\Phi_{xyy}$, by an amount $d$. Its 0th, 1st and 2nd moments all vanish, so on the spot it touches neither density, nor momentum, nor stress. On D2Q9 that vector looks like this.

One streaming step changes the story. Each component flies to a different neighbour. Collect only the $x$-momentum and this is what is left.

A second central difference. It is exactly zero when $d$ is uniform, and grows with the square of the wavenumber when $d$ oscillates. A uniform waveform kills the cell-to-cell change. Shorten a sine and that change grows with the square of the wavenumber. The sum of the cells stays zero.

A channel with a pinned force pattern and a cross-flow#

Verification needs a problem with an exact solution. Pin a sinusoidal $x$-body-force in a periodic box, then blow a uniform cross-flow $V$ through it in $y$. Porous walls injecting on one side and sucking on the other.

$\nu$ is kinematic viscosity, $k$ the wavenumber, $a_0$ the force amplitude. A single mode gives a closed form.

The cross-flow does two things. It shaves the amplitude and it lags the phase by $\arctan\!\left(V/\nu k\right)$. At the same time $u_y = V$ stays alive, so $d = F_x u_y^2$ oscillates in space. That is the switch for the leak above.

Four combinations on the same lattice#

A 16-node periodic lattice, $\tau = 1$, $a_0 = 4\times10^{-3}$, four runs. Plain Guo, the force third moment filled, the equilibrium third moment filled, and both.

import numpy as np
 
CX  = np.array([0, 1, 0,-1, 0, 1,-1,-1, 1], float)
CY  = np.array([0, 0, 1, 0,-1, 1, 1,-1,-1], float)
WT  = 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 hermite3_vectors():
    """Zero 0th-2nd moments, unit third moment for xxy / xyy only."""
    rows = np.vstack([np.ones(9), CX, CY, CX*CX, CX*CY, CY*CY, CX*CX*CY, CX*CY*CY])
    pinv = np.linalg.pinv(rows)
    return pinv @ np.eye(8)[6], pinv @ np.eye(8)[7]
 
VXXY, VXYY = hermite3_vectors()
 
def feq_second_order(rho, ux, uy):
    cu = CX[:, None, None]*ux + CY[:, None, None]*uy
    return WT[:, None, None]*rho*(1 + cu/CS2 + cu*cu/(2*CS2**2)
                                  - (ux*ux + uy*uy)/(2*CS2))
 
def force_populations(ux, uy, fx, fy, tau):
    """Guo forcing: first moment F, second moment uF+Fu."""
    cF  = CX[:, None, None]*fx + CY[:, None, None]*fy
    gxx, gyy, gxy = ux*fx, uy*fy, 0.5*(ux*fy + uy*fx)
    quad = (gxx*(CX[:, None, None]**2 - CS2) + gyy*(CY[:, None, None]**2 - CS2)
            + 2*gxy*(CX*CY)[:, None, None])
    return WT[:, None, None]*(1 - 1/(2*tau))*(cF/CS2 + quad/CS2**2)
 
def analytic_wave(ny, nu, V, a0):
    """Periodic steady solution of nu u'' - V u' + a0 sin(ky) = 0."""
    k = 2*np.pi/ny
    y = np.arange(ny, dtype=float)
    c1 = a0/(nu*k*k + V*V/nu)
    return c1*np.sin(k*y) - (V*c1/(nu*k))*np.cos(k*y)
 
def run_cross_flow(ny, tau, V, a0, cubic_eq, cubic_force, nstep=6000, tol=1e-12):
    nx = 3
    k  = 2*np.pi/ny
    fx = np.tile(a0*np.sin(k*np.arange(ny)), (nx, 1))
    fy = np.zeros((nx, ny))
    rho = np.ones((nx, ny)); ux = np.zeros((nx, ny)); uy = np.full((nx, ny), V)
    f = feq_second_order(rho, ux, uy)
    half, prev = 1 - 1/(2*tau), ux.copy()
    for s in range(nstep):
        rho = f.sum(0)
        ux  = ((CX[:, None, None]*f).sum(0) + 0.5*fx)/rho
        uy  = ((CY[:, None, None]*f).sum(0) + 0.5*fy)/rho
        feq = feq_second_order(rho, ux, uy)
        if cubic_eq:                      # fill the cubic defect rho*u*u*u of the equilibrium
            feq = feq + VXXY[:, None, None]*(rho*ux*ux*uy) \
                      + VXYY[:, None, None]*(rho*ux*uy*uy)
        src = force_populations(ux, uy, fx, fy, tau)
        if cubic_force:                   # fill the cubic defect F*u*u of the force term
            src = src + half*(VXXY[:, None, None]*(2*fx*ux*uy + fy*ux*ux)
                              + VXYY[:, None, None]*(fx*uy*uy + 2*fy*ux*uy))
        f = f - (f - feq)/tau + src
        for i in range(9):
            f[i] = np.roll(f[i], (int(CX[i]), int(CY[i])), axis=(0, 1))
        if s % 1000 == 999:
            if np.max(np.abs(ux - prev)) < tol:
                break
            prev = ux.copy()
    return ux[0]
 
print("D2Q9 lock check:  max|cx^3 - cx| =", float(np.abs(CX**3 - CX).max()))
for tag, mom in (("xxy", CX*CX*CY), ("xyy", CX*CY*CY), ("xxx", CX**3), ("xy", CX*CY)):
    print(f"  correction vector v_xxy :  sum v_i c_{tag:3s} = {float(mom @ VXXY): .6f}")
 
NY, TAU = 16, 1.0
NU, A0  = CS2*(TAU - 0.5), 0.004
print(f"\nny={NY}  tau={TAU}  nu={NU:.5f}  k={2*np.pi/NY:.4f}  a0={A0}")
print("\n   V    |u|max   Guo       +force    +feq3     +both")
for V in (0.05, 0.10, 0.20):
    ref = analytic_wave(NY, NU, V, A0)
    errs = []
    for ce, cf in ((0, 0), (0, 1), (1, 0), (1, 1)):
        u = run_cross_flow(NY, TAU, V, A0, bool(ce), bool(cf))
        errs.append(100*np.linalg.norm(u - ref)/np.linalg.norm(ref))
    print(f"{V:5.2f}  {np.abs(ref).max():.5f}  " + "  ".join(f"{e:7.4f}%" for e in errs))

D2Q9 lock check:  max|cx^3 - cx| = 0.0
  correction vector v_xxy :  sum v_i c_xxy =  1.000000
  correction vector v_xxy :  sum v_i c_xyy = -0.000000
  correction vector v_xxy :  sum v_i c_xxx =  0.000000
  correction vector v_xxy :  sum v_i c_xy  = -0.000000
 
ny=16  tau=1.0  nu=0.16667  k=0.3927  a0=0.004
 
   V    |u|max   Guo       +force    +feq3     +both
 0.05  0.12258   0.5993%   0.5915%   0.0181%   0.0143%
 0.10  0.08374   1.6422%   1.6207%   0.0295%   0.0375%
 0.20  0.04826   3.7307%   3.6840%   0.0531%   0.1503%

The gap between +feq3 and +both in that table grows as the wavelength shrinks. The waveform toggle and the wavelength slider are the bars for that leak.

On uniform the green bars sit at zero. Shorten the wavelength on sine and the bars grow with the square of the wavenumber, while the sum of the cells stays zero.

The correction that turned 0.053% into 0.150%#

Read the table one row at a time and the ordering shows up. Numbers are at $V=0.20$.

Filling the third moment of the force buys 1.3% at best. Filling the third moment of the equilibrium cuts the error 70-fold. Put both in and it gets worse again. The three numbers tell one story.

The reason lives in the order of the Chapman-Enskog expansion. The second moment of the force enters the $O(\varepsilon)$ stress equation directly. The third moment of the force enters the $O(\varepsilon)$ third-moment equation, descends from there into the stress, and from the stress into momentum. It arrives in the momentum equation at $O(\varepsilon^3)$, that is, at Burnett order. For an isothermal LBM aiming at Navier-Stokes, that term is not something to match.

The third moment of the equilibrium, by contrast, enters the $O(\varepsilon)$ stress directly. That is why the cubic defect of the equilibrium dominates.

The path by which the correction manufactured error#

So why does filling the force third moment actively hurt? The leak formula above answers it. Since $d = F_x u_y^2 = a_0 V^2 \sin(ky)$, streaming manufactures a spurious body force of this size.

Put in $V=0.2$, $k=0.3927$, $\tau=1$ and it is $0.154\%$ of $a_0$. The largest gap measured between the two solutions in the code was $0.156\%$ of $|u|_{\max}$. The two numbers agree. The lattice faithfully delivered a term the momentum equation never asked for.

The $V^2 k^2$ dependence shows up cleanly too. Raising $V$ through $0.05 \to 0.10 \to 0.20$, the gap between the two solutions went $0.010\% \to 0.039\% \to 0.156\%$. A factor of four each time.

Where the third moment genuinely starts to matter#

None of this makes the derivation of $\Phi_{\alpha\beta\gamma}$ pointless. There are places that need it.

In thermal or compressible LBM that solves an energy equation, the third moment is the heat flux. There the third moment of the force enters the $O(\varepsilon^2)$ energy equation directly. No trip to Burnett order required. Multiphase LBM that builds interfaces from a non-ideal equation of state is similar. The force is large and spatially steep. Once $k$ approaches one, $V^2k^2$ is no longer small.

At $V=0.20$, plain Guo is 3.7307%. Filling only the force third moment gives 3.6840%. The difference is 0.0467 percentage points, about 1.3% of 3.7307. Equilibrium alone is 0.0531%. Divide 3.7307 by 0.0531 and the factor is about 70. Both together give 0.1503%. Divide 0.1503 by 0.0531 and the factor is about 2.83, the 2.8× in the table. At $V=0.05$, both together is 0.0143%, smaller than equilibrium alone at 0.0181%. At $V=0.10$ the order flips, 0.0375% against 0.0295%. At $V=0.20$ the gap is 0.0972 percentage points.

Plug $V=0.2$, $k=0.3927$, $\tau=1$ into the leak. Then $1-1/(2\tau)$ is $1/2$. Times the leading $1/2$ that is $1/4$. $V^2=0.04$, and $k^2$ is about $0.1542$. One quarter times 0.04 times 0.1542 is 0.001542, which is 0.154% of $a_0$. The code measured 0.156% of $|u|_{\max}$. The two sit within 0.002 percentage points. Doubling $V$ from 0.05 to 0.10 to 0.20 multiplies $V^2$ by four, and the gaps 0.010%, 0.039%, 0.156% do the same. $|u|_{\max}$ falls through 0.12258, 0.08374, 0.04826. The denominator of $C_1$ is $\nu k^2+V^2/\nu$ with $\nu=0.16667$, so a larger $V$ shaves the amplitude.

Where the missing half of the force goes is what the side-by-side of Guo, Shan–Chen, and EDM splits on the same factor $(1-1/2\tau)$. Count the target equation, and the order at which the term arrives, before adding the line.

Why you cannot half-fill a term the lattice does not carry#

The real lesson of this experiment sits ahead of the order counting. $\Phi_{\alpha\beta\gamma}$ is a tensor. D2Q9 grants $xxy$ and $xyy$ and locks $xxx$ and $yyy$. Fill half of a tensor and you have inserted a correction with broken isotropy. Errors that were cancelling can surface, and here they did.

When you find a term the lattice cannot carry, there are two options. Enlarge the velocity set until you own that moment (D2Q17, D2Q37), or leave it alone entirely. There is no middle. Before adding one line to a forcing term, count whether the rest of that line exists on your lattice.