LBM 里那一半的力去哪了 — 四种 forcing 格式与 1−1/(2τ)

1993 年，Shan 和 Chen 在格子上把水和油分了开来。处方很朴素。只是在碰撞所用的平衡速度上加了一项 $\tau\mathbf{F}/\rho$，仅此而已。1998 年，He 把力直接挂进了平衡分布。2002 年，Guo 指出前面这些处方在应力上留了误差，于是在力项前面添上 $(1-\frac{1}{2\tau})$。2004 年，Kupershtokh 只用两个平衡分布之差就做成了同一件事。

同一个问题有四个答案。那么到底该用哪一个。今天把这个问题直接摆到 D2Q9 通道上问了一遍。答案与预想不同。在定常流动里，无论用哪种格式，误差直到小数点后第四位都完全一样。真正的陷阱不在挑选格式这一步。推导指出位置，数字量出大小。

矩旋钮同时改格式、$\tau$、$|u|$ 和 prefactor，M0、M1、M2 三条柱的高度跟着变。

旋钮怎么转，M0 都停在 0。$|u|$ 推上去，只有 plain 的 M2 够不着目标。关掉 prefactor，M1 和 $\Delta j$ 一起走样。

从连续方程里截出力项

带外力的 Boltzmann 方程会多出一项。

$f$ 是分布函数，$\boldsymbol{\xi}$ 是粒子速度，$\mathbf{F}$ 是单位体积的体积力，$\lambda$ 是弛豫时间。

麻烦出在第三项。速度空间的导数 $\nabla_{\boldsymbol{\xi}} f$ 在格子上并不存在。我们手里只有九个方向的离散速度。于是把 $f$ 换成 $f^{\rm eq}$，把这个导数解析地处理掉。

$c_s$ 是格子声速（D2Q9 中 $c_s^2=1/3$），$\mathbf{u}$ 是流体速度。

把 $f$ 换成 $f^{\rm eq}$ 到底行不行。低马赫数下行。非平衡部分 $f^{(1)}$ 是 Knudsen 数（平均自由程/特征长度）的一阶量，与力相乘后就落到可以忽略的阶。代价是这次替换后来会在应力上留下残余误差。各种格式分道扬镳的地方，正是这里。

把这个表达式投影到离散速度 $\mathbf{c}_i$ 上，并把 Hermite 级数截到二阶。于是得到格子上的力项。

$w_i$ 是格子权重。偏偏截在二阶，理由很明确。恢复 Navier–Stokes 所需要的矩，到二阶为止。

先把三个矩确认下来。

零阶矩是质量。力不制造质量，所以必须为 0。一阶矩是动量，二阶矩是应力。

梯形积分留下的那一半

现在把离散速度 Boltzmann 方程沿特征线积分 $\delta t$。碰撞项和力项用梯形法则积分，可得二阶精度。代价是右端混进了 $t+\delta t$。

上标 $n$ 和 $n+1$ 分别表示 $t$ 与 $t+\delta t$。就这样写是隐式的。每一步都得解联立方程。

LBM 用一次变量变换抹掉了这个隐式性。

改用 $\bar f_i$ 重写，并令 $\tau = \lambda/\delta t + 1/2$，就得到熟悉的那一行。

这次变换留下了两处痕迹。一处是力前面挂上的 $(1-\frac{1}{2\tau})$。另一处不那么显眼。我们存进数组的是 $\bar f_i$，不是 $f_i$。于是动量也偏掉了一半。

这两处是一体的。它们从同一次变换里一起掉出来。梯形积分留在黏性里的半步也从同一条根上长出来。只顾上一处、漏掉另一处，账就对不上。差多少，下面来量。

四种格式对照表

记 $\mathbf{j}=\sum_i \bar f_i\mathbf{c}_i$ 为存储分布函数的生动量，$\delta\mathbf{u}=\delta t\,\mathbf{F}/\rho$ 为力在一步内给出的速度增量。

把这五行竖着读，共同点就出来了。在乘上系数之前，力项的一阶矩全都精确等于 $\mathbf{F}$。本来就是这么设计的，理所当然。分道扬镳的是二阶矩，也就是应力。

可是变量变换所要求的目标值不是 $\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u}$，而是 $(1-\frac{1}{2\tau})(\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u})$。明确扛起这个系数的，只有最后一行。其余几种都带着 $\frac{1}{2\tau}$ 那么大的应力偏差往前走。

偏差的量级是力与速度之积，即 $u F$ 阶。低马赫数下 $u$ 很小，这一项本身就小。Shan–Chen 与 EDM 留下的 $\mathbf{F}\mathbf{F}$ 项更小，是 $F^2$ 阶。所以在用弱重力驱动的问题里，用哪一种都差不多。差别显现的地方，是力很强的时候（$\tau\mathbf{F}/\rho$ 可与 $\mathbf{u}$ 相比），或者像多相流那样力在界面处急剧变化的位置。

M1 的目标不是 $\mathbf{F}$，而是 $(1-\frac{1}{2\tau})\mathbf{F}$。量的是乘上系数之后的力项。缺掉的 $\mathbf{F}/(2\tau)$ 由平移半格的平衡态还回来，$\Delta j$ 落到 $\mathbf{F}$。分道的只有 M2。plain 错在形状，系数救不回来。把 $\tau$ 推向 2，Shan–Chen 的偏差会长大。

定常 Poiseuille 分辨不出它们#

去看最常见的验证问题。用体积力推动上下壁面之间的通道。定常态的解已知是抛物线。

$H$ 是通道高度，$\nu$ 是运动黏度。用 half-way bounce-back 时，壁面落在格点外半格处，因此 $H$ 等于流体节点数，第 $j$ 个节点的坐标是 $y=j+0.5$。

换成别的格式，抛物线幅值也不动。33 个流体节点、$F_x = 10^{-5}$、40,000 步之后的最大相对误差如下。

列与列之间没什么可比的。数值相同。剩下的误差不是来自格式，而是来自 bounce-back 的离散化。

理由就在二阶矩进入方程的位置上。这个流动是定常且单向的。速度只有 $u_x(y)$ 一个分量，力也只有 $x$ 分量。各格式分道扬镳的那一项是 $\mathbf{u}\mathbf{F}$，也就是 $xx$ 分量。可真正支配这个流动动量平衡的是 $xy$ 剪应力。力造成的 $xx$ 偏差没有通道进入平衡式。

$\partial_t\mathbf{u}=0$ 也在一起起作用。力项的误差没有随时间累积显现的余地。要看出格式差异，就得走向非定常流动或空间上变化的力。至少用这个测试，分不出这五种。

这不是坏消息。这是挑选验证问题时该知道的信息。因为它意味着，Poiseuille 通过了并不能说明 forcing 的实现是对的。

把那一半数两遍，力就变大

那么实际出错的是什么。是前面说过属于一体的那两样，即速度的半格修正与力前面的系数。

来数一数一步之内实际注入的动量。在平衡速度上加了 $\delta t\mathbf{F}/2\rho$，碰撞每步会把其中的 $1/\tau$ 搬运过去。再加上力项直接给出的那一份。把是否开启半格修正记为 $s\in\{0,1\}$，力前面的系数记为 $g$，则

它必须等于 $\delta t F$。条件只有 $g = 1 - s/(2\tau)$ 这一个。两个开关并不独立。

一旦错开，幅值就按可预测的比例偏掉。固定 Guo 力项，只改这两个开关所测得的幅值比如下。

测量与预测在小数点后第三位上吻合。$\tau=0.6$ 时只开半格修正，力会大出 83%。只开系数，力会缩到六分之一。

这不是什么微妙的精度问题。从论文里原样抄来 Guo 的 $F_i$，速度却仍旧沿用旧代码里的 $\mathbf{j}/\rho$，发生的正是这件事。要是觉得黏性不对劲而开始动 $\tau$，就会陷得更深。因为每改一次 $\tau$，误差比例也跟着一起动。

有办法从症状上区分。把网格加倍而误差比例不下降，那就不是离散化误差。$\tau$ 越靠近 1，比例越收敛到 1.5 附近；$\tau$ 越降到 0.5，比例越像发散一样变大 — 这是 $1+\frac{1}{2\tau}$ 的指纹。反过来，把力加大而比例不变，那就是记账问题，不是非线性问题。三条同时出现，就别怀疑格式了，先去打开速度的定义看看。

$\tau=0.6$ 时 $2\tau=1.2$，$1/(2\tau)$ 等于 $5/6$，也就是 0.8333。两个都开意味着 $s=1$，$g=1-1/(2\tau)=1/6$。一步注入量除以 $\delta t F$，得到 $(1/\tau)\cdot(s/2)+g$。代入就是 $(1/0.6)\cdot 0.5 + 1/6 = 5/6 + 1/6 = 1$。和目标相同。表里 both 一列的 0.9991，只比这个 1 少了格子那一点。

只开修正时 $g$ 固定为 1。同一格变成 $5/6+1=1.8333$。表上是 1.8316。差 0.0017，大约是预测的 0.09%。只开系数时 $s=0$，$g=1/6$，注入量是 0.1667。测量值 0.1666 在第四位上碰上。

$\tau=1$ 时 $1/(2\tau)=1/2$。只开修正是 1.5，只开系数是 0.5。表上是 1.5002 和 0.5005。两个都开是 1.0003，离目标 1 还剩 0.03%。$\tau=1.8$ 时 $1/(2\tau)=1/3.6=0.2778$。预测 1.2778 和 0.7222，测量是 1.2867 和 0.7280。$(1.2867-1.2778)/1.2778$ 约 0.007，也就是 0.7%。同一 $\tau$ 下格式相对误差 0.7358%，大小叠在一起。

相对误差表也是这么读。$\tau=0.6$ 为 0.0875%，$\tau=1.0$ 为 0.0306%，$\tau=1.8$ 为 0.7358%。$\nu=c_s^2(\tau-1/2)$，而 $c_s^2=1/3$，所以黏度分别是 $1/30$、$1/6$、$13/30$。$(\tau-1/2)$ 之比是 0.1 比 0.5 比 1.3，也就是 1 比 5 比 13。峰值为 $u(H/2)=F_x H^2/(8\nu)$。同样的 $F_x=10^{-5}$、$H=33$ 下，高度与 $\nu$ 成反比。$\tau=1.8$ 的峰是 $\tau=0.6$ 的 $1/13$。百分比在那里变大，五格仍然是同一个数。那不是格式分开的痕迹。

Python — 故意把两个开关错开#

下面是把两个开关提成参数的 D2Q9 通道求解器。沿 $x$ 方向均匀，所以只保留一列。上表的五种格式只需换掉 source_terms 和放进平衡态的速度即可，所以这里固定用 Guo 力项一种。

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 力项，尚未乘上 (1 - 1/2tau) 增益。"""
    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   # 开关 1
        uy = jy / rho
        fpost = f - (f - lattice_equilibrium(rho, ux, uy))/tau \
                + gain*source_terms(rho, ux, uy, fx)        # 开关 2
        for i in range(9):                                   # 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)      # 物理速度，始终带半格修正
 
 
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})")

输出如下。

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)

把通道实验的格式按钮、$\tau$、半格修正和 prefactor 调到上面输出里的 both、shift_only、gain_only。

四个格式按钮都按过，青绿色曲线仍贴在虚线抛物线上。只关一个开关，$\tau=0.6$ 时幅值会跳向 1.83 或 0.17。

$\tau=1.8$ 时与预测差了 0.7%，而同样条件下格式本身的离散化误差是 0.74%。量级相同。这说明偏差的来源不是记账，而是格子。

只留半格修正时，力会变成1.83倍#

$\tau=0.6$ 的 shift_only 是 1.8316。预测 $1+\frac{1}{2\tau}=1.8333$，到小数点后第三位一致。只抄 Guo 的力项、速度仍用 $\mathbf{j}/\rho$，就会落在这一格：修正开着，prefactor 关着。网格加倍，这个比例也不降。$\tau$ 收到 1 附近，比例靠近 1.5；收到 0.5，比例变大。这是 $1+\frac{1}{2\tau}$ 的指纹。

五种格式的一阶矩都是 $\mathbf{F}$。分叉在二阶，目标里带着 $(1-\frac{1}{2\tau})$。定常通道没有把这份 $xx$ 偏差送进剪切平衡的通路，所以五列停在 0.0875%、0.0306%、0.7358% 同一行上。流动非定常，或力在界面上陡变时，再去看 M2。二阶之上把三阶也填满，误差反而变大，记在力项的三阶矩界限里。