Raise the Velocity, Lose 27% of Your Diffusion — The Extra Flux in LBM Convection-Diffusion

In lattice Boltzmann, one relaxation time sets the diffusivity. $D = c_s^2(\tau - 1/2)\delta t$. Short line. No velocity in it. Does the same $D$ survive once the flow is on? Under a uniform advection speed $u$ it drops to $D(1 - u^2/c_s^2)$. At $u = 0.3$ that is 27% gone. The Chapman–Enskog expansion is where the extra term leaks out, and one source term is what puts it back.

Phase-field multiphase work is where the miss hurts. Cahn–Hilliard or Allen–Cahn on a lattice Boltzmann kernel ties interface thickness to the mobility. A mobility 27% off gives a thickness 27% off, and the surface-tension coefficient follows it.

The diffusivity is set, yet the interface keeps thinning#

Start by looking at it. Below is the one-dimensional convection-diffusion equation

solved with a D1Q3 lattice Boltzmann scheme. The initial condition is a single Gaussian; the exact solution is also a Gaussian, widening as $\sigma^2(t) = \sigma_0^2 + 2Dt$. Push u from 0 to 0.30 and the hump narrows against the dashed curve; tau does not move that ratio.

At u = 0.30 the solid curve (computed) stands sharper than the dashed one (exact). The $\sigma^2$ trace below misses the dashed slope $2D$. Park tau at 0.6 or at 1.5 and the deficit ratio stays in the same slot.

D1Q3 has only three moments to honor#

Unlike a Navier–Stokes lattice Boltzmann scheme, a distribution $g_i$ built for scalar transport has fewer moments to satisfy. The equilibrium Guo used in his 2009 nonlinear convection-diffusion model is this:

Here $w_i$ are the lattice weights, $\mathbf{c}_i$ the lattice velocities, and $c_s^2 = 1/3$ the squared lattice sound speed. It satisfies three moment constraints.

The third one is where this parts ways with a flow equilibrium: there is no $\phi\mathbf{u}\mathbf{u}$ term. The second-order term in velocity was never put in. A scalar equation has no stress tensor, so an isotropic second moment is enough, and D2Q5 can replace D2Q9. The half step inside $\tau - 1/2$ in front of the diffusivity is the same leftover the second-order expansion of the discrete LBE leaves in the viscosity.

That looks sufficient, and at $u = 0$ it is exact. The catch is that the truncation is not free.

The fourth term Chapman-Enskog leaves behind#

Write the BGK lattice Boltzmann equation with a source $S_i$ attached:

Expand $g_i = g_i^{(0)} + \varepsilon g_i^{(1)} + \varepsilon^2 g_i^{(2)}$ and split the time derivative as $\partial_t = \varepsilon\partial_{t_1} + \varepsilon^2\partial_{t_2}$. The zeroth moment at order $\varepsilon$ hands back the advective part of the target equation unchanged.

The first moment of that same order-$\varepsilon$ equation is where everything is decided, because it fixes the flux carried by $g^{(1)}$.

The second term in the bracket, $c_s^2\nabla\phi$, is the diffusive flux we asked for. The first one is riding along with it. In a flow scheme the $\rho\mathbf{u}\mathbf{u}$ term in the equilibrium second moment cancels most of it; the scalar model has no such term. So it survives.

Collecting the zeroth moment at order $\varepsilon^2$ gives the final form.

$D = c_s^2(\tau - 1/2)\delta t$ appears as expected. The second term on the right is a flux nobody ordered. It vanishes in a genuine steady state, where $u$ is constant in time and $\phi$ has stopped changing. But as long as $\phi$ is still moving — as long as the computation is running — $\partial_t(\phi\mathbf{u})$ is not zero.

Under uniform $u$ that term is negative diffusion#

Take the simplest case, $u$ constant in space and time. Then $\partial_t(\phi u) = u\,\partial_t\phi$, and to leading order $\partial_t\phi \simeq -u\,\partial_x\phi$. Substituting,

A diffusion term with the sign reversed. Adding it to the physical one,

The deficit is one square. Lattice sound speed squared is one third, so a speed of 0.3 squares to 0.09 and triples to 0.27. At relaxation time 1 the dialed diffusivity is one third times one half, 0.166667. What remains is 0.166667 times 0.73, which is 0.121667. The gap, 0.045, is the dialed value times 0.27. The $u = 0.30$ row of the printout is that subtraction. At low speed the square shrinks. A speed of 0.05 squares to 0.0025, triples to 0.0075, and the deficit is 0.75%. The hump sits next to the exact curve. Easy to miss by eye.

The ratio refuses to track relaxation time for the same reason. The factor of relaxation time minus one half multiplies both the physical diffusion and the extra term. Divide the two and the factor drops out. The ratio does not depend on $\tau$. Raising $\tau$ from 0.6 to 1.5 takes $D$ from 0.033333 to 0.333333, a factor of ten. At speed 0.25, three times 0.0625 is 0.1875, and one minus that is 0.8125, fixed on all four rows. Only the absolute gap grows, from about 0.006 to about 0.062. There is no room to bury the error under a larger diffusivity.

When speed reaches the lattice sound speed, the reciprocal square root of one third, near 0.577, the parenthesis hits zero. Past that point $D_{\text{eff}}$ is negative. The hump sharpens itself instead of spreading, and the distribution blows up.

$1 - 1/(2\tau)$ shows up again#

Now fix the coefficient of $S_i$. For the spurious term and the source term to cancel in the equation above,

Only one simple form satisfies this while also keeping $\sum_i S_i = 0$.

There is $1 - 1/(2\tau)$ again. This factor sits on the same line as the third-moment limit a forcing term has to respect. On a discrete-time lattice a source acts twice — once through $g^{(1)}$, once through the second-order term of the Taylor expansion — and the factor of $1/2$ is the residue of that double counting.

What happens if you drop the coefficient and set $\sum_i \mathbf{c}_i S_i = \partial_t(\phi\mathbf{u})$? At $\tau = 1$ you apply exactly twice the correction, and a 27% deficit becomes a 27% excess. The magnitude of the error is unchanged, so a log-log convergence plot will not catch it.

Computing $\partial_t(\phi\mathbf{u})$ in code means storing $\phi u$ from the previous step and taking a backward difference. One extra array is the entire cost.

Measuring $D_{\text{eff}}$ in 60 lines of Python#

Rather than argue, measure. Advect a Gaussian, fit the growth slope of its second moment $\sigma^2$ by least squares, and that slope is $2D_{\text{eff}}$.

import numpy as np
 
CS2 = 1.0 / 3.0
C = np.array([0, 1, -1])
W = np.array([2 / 3, 1 / 6, 1 / 6])
 
 
def d1q3_equilibrium(phi, u):
    """g_i^eq = w_i phi (1 + c_i u / cs^2) — an equilibrium fixing only three moments"""
    return np.stack([W[i] * phi * (1.0 + C[i] * u / CS2) for i in range(3)])
 
 
def gaussian_moments(x, phi):
    m0 = phi.sum()
    mean = (x * phi).sum() / m0
    return mean, (((x - mean) ** 2) * phi).sum() / m0
 
 
def run_cde_lbm(L, steps, tau, u, sigma0, x0, corrected):
    x = np.arange(L, dtype=float)
    phi = np.exp(-((x - x0) ** 2) / (2 * sigma0**2))
    g = d1q3_equilibrium(phi, u)
    phi_old = phi.copy()
    hist = []
    for n in range(steps + 1):
        if n % 100 == 0:
            hist.append((n, gaussian_moments(x, phi)[1]))
        src = np.zeros_like(g)
        if corrected and n > 0:
            # S_i = w_i (1 - 1/(2 tau)) c_i d_t(phi u) / cs^2,  dt = 1
            dt_phiu = (1.0 - 1.0 / (2 * tau)) * u * (phi - phi_old)
            for i in range(3):
                src[i] = W[i] * C[i] * dt_phiu / CS2
        geq = d1q3_equilibrium(phi, u)
        g = g - (g - geq) / tau + src          # collision
        for i in range(3):
            g[i] = np.roll(g[i], C[i])         # streaming
        phi_old = phi
        phi = g.sum(axis=0)
    return np.array(hist)
 
 
def fit_diffusivity(hist):
    """read D_eff off the slope of sigma^2 = sigma0^2 + 2 D_eff t"""
    return np.polyfit(hist[:, 0], hist[:, 1], 1)[0] / 2.0
 
 
L, STEPS, SIG0, X0 = 800, 1600, 10.0, 80.0
 
tau = 1.0
D = CS2 * (tau - 0.5)
print(f"tau = {tau},  D = cs^2 (tau-1/2) = {D:.6f},  cs^2 = {CS2:.6f}")
print(f"{'u':>6} {'u^2/cs^2':>9} | {'D_eff (no src)':>14} {'ratio':>7} {'1-u^2/cs^2':>11} |"
      f" {'D_eff (src)':>12} {'ratio':>7}")
for u in [0.05, 0.10, 0.20, 0.30]:
    d_raw = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, False))
    d_fix = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, True))
    print(f"{u:>6.2f} {u * u / CS2:>9.4f} | {d_raw:>14.6f} {d_raw / D:>7.4f} {1 - u * u / CS2:>11.4f} |"
          f" {d_fix:>12.6f} {d_fix / D:>7.4f}")
 
u = 0.25
print(f"\nu = {u} fixed, tau sweep   (theory: ratio = 1 - u^2/cs^2 = {1 - u * u / CS2:.4f}, tau-independent)")
print(f"{'tau':>6} {'D':>10} | {'D_eff (no src)':>14} {'ratio':>7} | {'D_eff (src)':>12} {'ratio':>7}")
for tau in [0.6, 0.8, 1.0, 1.5]:
    D = CS2 * (tau - 0.5)
    d_raw = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, False))
    d_fix = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, True))
    print(f"{tau:>6.1f} {D:>10.6f} | {d_raw:>14.6f} {d_raw / D:>7.4f} | {d_fix:>12.6f} {d_fix / D:>7.4f}")

The output:

tau = 1.0,  D = cs^2 (tau-1/2) = 0.166667,  cs^2 = 0.333333
     u  u^2/cs^2 | D_eff (no src)   ratio  1-u^2/cs^2 |  D_eff (src)   ratio
  0.05    0.0075 |       0.165417  0.9925      0.9925 |     0.166666  1.0000
  0.10    0.0300 |       0.161667  0.9700      0.9700 |     0.166666  1.0000
  0.20    0.1200 |       0.146667  0.8800      0.8800 |     0.166663  1.0000
  0.30    0.2700 |       0.121667  0.7300      0.7300 |     0.166658  0.9999
 
u = 0.25 fixed, tau sweep   (theory: ratio = 1 - u^2/cs^2 = 0.8125, tau-independent)
   tau          D | D_eff (no src)   ratio |  D_eff (src)   ratio
   0.6   0.033333 |       0.027096  0.8129 |     0.033345  1.0004
   0.8   0.100000 |       0.081258  0.8126 |     0.100006  1.0001
   1.0   0.166667 |       0.135417  0.8125 |     0.166661  1.0000
   1.5   0.333333 |       0.270794  0.8124 |     0.333275  0.9998

Raise u and only the rose ghost flux grows as $u^2$. Drag alpha from 0 to 1 and the source term cancels that width.

At alpha = 0 the rose lobe mirrors the blue physical flux upside down. At alpha = 1 green lands on rose and the amber sum returns to the blue curve. That is the column where the source is on and ratio snaps to 1.0000.

In the first table the ratio column matches 1-u^2/cs^2 to four decimal places. The prediction is not an estimate; it is the exact leading order. Switch the source term on and all four velocities return to 1.0000.

The second table bites harder. Sweeping $\tau$ from 0.6 to 1.5 changes $D$ by a factor of ten, and the deficit ratio moves only from 0.8129 to 0.8124. Trying to bury the error under a larger diffusivity does not work: raise $D$ tenfold and the amount that disappears rises tenfold too.

The lattice velocity budget was already spent#

$u^2/c_s^2$ is the squared lattice Mach number. At $u = 0.1$, $0.01 \times 3 = 0.03$, and 3% of the diffusion is gone. At $u = 0.2$, $0.04 \times 3 = 0.12$, which is 12%. At or below 0.05 the loss is 0.75% and it hides in the grid error. The habit of pinning $u$ near 0.1 because of compressibility has not seen the scalar bill yet. If an interface thickness or a Sherwood number is going out as a number, the source term goes on the moment $u$ passes 0.1. Below that, skipping one array still leaves the printed 0.9925 standing.

Phase-field work cannot trust $u^2/c_s^2$ alone. The extra term is the whole $\partial_t(\phi\mathbf{u})$, and the square holds only when $u$ is uniform and steady. Near an interface $u$ jumps and $\partial_t u$ is not zero, so the $\phi\,\partial_t u$ piece survives too. $\partial_t(\phi\mathbf{u})$ goes in as a backward difference of $\phi u$ from the previous step. The $\tau$ inside $1 - 1/(2\tau)$ is the relaxation time of the first moment. How viscosity drops when the velocity frame is shifted is the same kind of bill as the missing third moment of the equilibrium. If the mobility arithmetic checks out and the interface still thins, $\tau$ was not the miss. The flow took the diffusion budget first.