把参考系推到0.2，粘性少了5.97% — LBM平衡分布缺失的三阶矩

LBM的速度上限不是由稳定性决定的，而是由平衡分布能匹配的矩阶数先定下来的。

静止的水和流动的水，粘性应该一样吗

把同一种流体量两次。一次在静止的盒子里，一次让盒子以匀速平移。两次得到的粘性系数应该完全相同。伽利略不变性，也就是在两个匀速相对运动的坐标系中物理定律形式相同这一性质，要求如此。

在D2Q9格子玻尔兹曼方法(LBM)上照做这个实验，第二个数值会小5.97%。加密网格没有用。改变松弛时间 $\tau$ 也几乎不改变这个比例。本文把这5.97%一路追到平衡分布的一行矩条件上。先说结论：这不是稳定性问题，而是代数问题。格子速度集合根本造不出那个缺失的项。

把格子上的那一行做Taylor展开#

LBM求解的方程只有一行。

f_i(\mathbf{x} + \mathbf{c}_i \Delta t,\; t + \Delta t) - f_i(\mathbf{x}, t) = -\frac{\Delta t}{\tau}\left[ f_i - f_i^{eq} \right]

其中 $f_i$ 是沿格子速度 $\mathbf{c}_i$ 运动的分布函数， $f_i^{eq}$ 是当地的平衡分布， $\tau$ 是松弛时间。

把左边按 $\Delta t$ 展开，除一阶项外还留下二阶项。

\Delta t \left( \partial_t + \mathbf{c}_i \cdot \nabla \right) f_i + \frac{\Delta t^2}{2} \left( \partial_t + \mathbf{c}_i \cdot \nabla \right)^2 f_i + O(\Delta t^3) = -\frac{\Delta t}{\tau} \left( f_i - f_i^{eq} \right)

正是这个二阶项，在Chapman–Enskog展开(把分布函数按克努森数展成幂级数的多尺度方法)中产生了著名的 $\tau - 1/2$ 。同一展开的下一步会给出一个条件：要让得到的连续介质方程成为Navier–Stokes方程，平衡分布必须精确满足四个矩。

\sum_i f_i^{eq} = \rho, \qquad \sum_i f_i^{eq} c_{i\alpha} = \rho u_\alpha, \qquad \sum_i f_i^{eq} c_{i\alpha} c_{i\beta} = \rho u_\alpha u_\beta + \rho c_s^2 \delta_{\alpha\beta}

\sum_i f_i^{eq} c_{i\alpha} c_{i\beta} c_{i\gamma} = \rho c_s^2 \left( u_\alpha \delta_{\beta\gamma} + u_\beta \delta_{\gamma\alpha} + u_\gamma \delta_{\alpha\beta} \right) + \rho u_\alpha u_\beta u_\gamma

前三个决定质量、动量和压力。第四个，也就是三阶矩，用来替换粘性应力的时间导数。这里一旦错位，连续方程和Euler层面都完好无损，唯独粘性项被污染。

把速度滑块从0往上推，黄色Maxwell曲线会离开蓝色格子直线。

velocity u : 0.180 (Ma = 0.31)

Σ f c³ = 0.180000  |  Maxwell = 0.185832  |  gap = 0.005832  |  gap / u³ = 1.0000

间距是否跟着 $u$ 的奇数次幂走，用 sweep u 就能看出来。打开 speeds ±2，柱子变成五根，间距闭合为零。

$c_x^3 = c_x$ — 格子造不出来的那一项#

D2Q9的 $x$ 分量只有 $-1, 0, 1$ 三个值。这三个数的立方都等于它们自己。于是对任意分布 $f_i$ ，

\sum_i f_i c_{ix}^3 = \sum_i f_i c_{ix} = \rho u_x

恒成立。平衡分布怎么设计都一样。哪怕往九个格子里填随机数，等式依然成立。三阶矩早已被一阶矩锁死。

而Maxwell分布要求的值是 $\rho(u_x^3 + 3 c_s^2 u_x)$ ，在D2Q9上 $c_s^2 = 1/3$ ，所以 $3c_s^2 = 1$ ，目标值折叠成 $\rho(u_x + u_x^3)$ 。两者之差恰好是 $\rho u_x^3$ 这一项。不是"大致是这个量级"，而是不多不少就这一项。

缺失项在动量方程里留下的样子

缺掉的三阶矩会直接加到粘性应力上。把Chapman–Enskog展开推到底，应力里会多出一项。

\sigma_{\alpha\beta}^{\text{err}} = -\left( \tau - \tfrac{1}{2} \right) \Delta t \, \partial_\gamma \left( \rho\, u_\alpha u_\beta u_\gamma \right)

设平均流 $U$ 沿 $x$ 方向流动，其上叠加一个小扰动 $u'$ 。只保留对 $u'$ 的一阶部分，该项变为 $-(\tau - 1/2)\Delta t\, U^2 \partial_x^2 (\rho u'_\alpha)$ ，形式与粘性项完全一致。也就是说， $x$ 方向的粘性系数被整体改写。

\nu_{xx}^{\text{eff}} = \left( \tau - \tfrac{1}{2} \right) \Delta t \left( c_s^2 - U^2 \right)

相对误差是 $-U^2 / c_s^2$ 。注意 $\tau$ 在分子和分母中同时消去：无论怎样调节松弛，这个百分比都不变。而误差正比于速度的平方，用马赫数写就是 $-\mathrm{Ma}^2$ 。

把算术摊开。 $c_s^2 = 1/3$ 、 $U = 0.2$ 时 $U^2 = 0.04$ ， $U^2/c_s^2 = 0.12$ 。只看 $x$ 方向，粘性差12%。Taylor–Green满足 $k_x = k_y$ ，衰减率取一半，预测是6%。表里 $U_0 = 0.20$ 的 $\nu = 0.0940296$ ，对照理论值 $0.100000$ ，相对误差 $-5.9704\%$ 。6.00和5.97相差0.03个百分点。

晃动 $\tau$ 也停在同一处。 $\tau = 0.6$ 的理论粘性是 $c_s^2 \times 0.1 = 0.0333333$ ， $\tau = 1.0$ 是 $c_s^2 \times 0.5 = 0.1666667$ ，五倍。可 $U_0 = 0.15$ 的相对误差只从 $-3.4076\%$ 移到 $-3.3099\%$ 。把 $U^2 = 0.0225$ 代入 $-U^2/(2c_s^2)$ ，得到 $-0.03375$ ，也就是 $-3.375\%$ 。三行都贴着这个数。 $\tau - 1/2$ 如何进入粘性公式，写在把离散LBE展开到二阶的那一篇。这里只剩下那个因子在分子分母里被整项消掉。

亏损的大小已经印在审计列里。 $u = 0.20$ 时格子和是 $0.20000000$ ，Maxwell一侧是 $0.20800000$ ，间距 $0.00800000$ 。因为 $0.2^3 = 0.008$ ，比值是1.00000。降到 $u = 0.05$ ，间距变成 $0.00012500$ ，比值仍是1。 $u = 0.30$ 的间距 $0.02700000$ 等于 $0.3^3$ 。塞进随机数，两个和的差仍是 $0.000\mathrm{e}{+}00$ 。把帧速放慢和平衡装不下 $\rho u^3$ ，是同一项的两张脸。

用Python在两个参考系里测同一个涡#

把Taylor–Green涡放在64×64的周期格子上，再叠加一个匀速 $U_0$ 。用振幅衰减的对数斜率反算粘性系数。只改变 $U_0$ ，其余全部固定。

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):
    # 展开到二阶的标准平衡分布
    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):
    # 给同一个涡叠加匀速U0，再从衰减率把粘性测回来
    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()   # 扣掉平均流之后的脉动分量
            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碰撞
        for i in range(9):                            # 迁移
            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%

把 frame velocity U 调高，右侧的涡一边平移，衰减一边变慢。

frame velocity U : 0.18 (Ma = 0.31)relaxation τ : 0.80

ν = 0.10000  |  νeff = 0.09514  |  error -4.86 %  |  step 0  |  amplitude ratio 1.000

同一种流体，右边却更慢地淡掉。晃动 $\tau$ 滑块，下方 error 几乎不动。那就是表里的 $-3.41\%$ 和 $-3.31\%$ 。

误差挂在 $U^2$ 上，而不是 $\tau$ 上#

[1]最后一列全部是1.00000。这说明"缺失项恰好是 $\rho u^3$ "这一断言精确到小数点后五位。

[2]中静止参考系的误差是 $+0.0175\%$ ，这是测量方法本身的噪声基准线。把 $U_0$ 提高，误差依次走到 $-0.36\%$ 、 $-1.48\%$ 、 $-3.35\%$ 、 $-5.97\%$ 。最后一列是误差除以 $U_0^2$ ， $U_0$ 越小，它越贴近 $-1.5$ 。

这个 $-1.5$ 正是上一节的预测。相对误差 $-U^2/c_s^2$ 只作用在 $x$ 分量上，而Taylor–Green模式是 $k_x = k_y$ 的对角模式，衰减率取两个方向的平均。所以是一半，即 $-U^2/(2c_s^2) = -1.5 U^2$ 。在 $U_0 = 0.2$ 处，预测是 $-6.00\%$ ，实测是 $-5.97\%$ 。

[3]更值得注意。把 $\tau$ 从0.6提到1.0，粘性变成五倍，相对误差却只从 $-3.41\%$ 移动到 $-3.31\%$ 。也就是说，靠加大粘性并不能把这个误差埋掉。

三种修法，以及各自的账单

降低速度。 误差按 $U^2$ 走，格子速度减半，误差就降到四分之一。代价是复现同样物理时间所需的步数增加。LBM里"马赫数保持在0.1以下"的经验法则通常被当作稳定性规则介绍，但真正先起作用的是这个缺失项。

加补偿项。 用差分计算 $\partial_\gamma(\rho u_\alpha u_\beta u_\gamma)$ ，以相反符号注入碰撞步。额外成本是一次梯度计算，同时格子玻尔兹曼完全局部碰撞这一优点会被削弱一点。即便如此，对旋转坐标系这类平均流很大的问题，它仍是最便宜的选项。

扩大速度集合。 包含 $\pm 2$ 的多速度格子满足 $c^3 \neq c$ ，因此可以精确匹配三阶矩。第一个viz里的 speeds ±2 按钮做的就是这个求解。代价是模板变宽、内存增加、边界处理更复杂，而且部分速度的权重可能变负，稳定性需要重新评估。

这个缺失项在真实计算中露面的地方

凡是平均流很大的场合都是候选。旋转机械的旋转坐标系、滑移网格、叠加在高速均匀来流上的湍流、贴在运动物体上的参考系。一份在静止流基准算例上表现完美的代码，偏偏在这些问题上把有效雷诺数算偏，就是这个症状。

按层级重取 $\tau$ 也抹不掉这项亏损，因为相对误差与 $\tau$ 无关。没有当作物性输入的量，会被格子结构自己定死。

诊断很便宜。给同一个算例叠一个平均流再跑一次，比较衰减率或阻力系数。差距若随 $U^2$ 增大，就是这一项。不必先去怀疑边界处理。

把 $u/c_s$ 压低，还是把速度扩到 $\pm 2$ #

D2Q9在三阶丢掉 $\rho u^3$ 。 $U_0 = 0.2$ 的账单是 $-5.97\%$ 。加密网格，或把 $\tau$ 从0.6抬到1.0，这个比例几乎不动。

能缩小它的只有三条路。压低 $u/c_s$ ，误差按 $U^2$ 下降。把亏损项 $\partial_\gamma(\rho u_\alpha u_\beta u_\gamma)$ 加回碰撞。或者扩大速度集合。三个速度为什么锁成 $c^3 = c$ ，接在厄密根挑选格子速度的那一篇。不改这个选择，参考系推到0.2，粘性还会少5.97%。