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

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

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

从连续方程里截出力项

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

\partial_t f + \boldsymbol{\xi}\cdot\nabla_{\mathbf{x}} f + \frac{\mathbf{F}}{\rho}\cdot\nabla_{\boldsymbol{\xi}} f = -\frac{1}{\lambda}\left(f - f^{\rm eq}\right)

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

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

\frac{\mathbf{F}}{\rho}\cdot\nabla_{\boldsymbol{\xi}} f \simeq -\frac{(\boldsymbol{\xi}-\mathbf{u})\cdot\mathbf{F}}{\rho c_s^2}\,f^{\rm eq}

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

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

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

F_i = w_i\left[\frac{\mathbf{c}_i-\mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot\mathbf{u})}{c_s^4}\mathbf{c}_i\right]\cdot\mathbf{F}

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

先把三个矩确认下来。

\sum_i F_i = 0,\qquad \sum_i \mathbf{c}_i F_i = \mathbf{F},\qquad \sum_i \mathbf{c}_i\mathbf{c}_i F_i = \mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u}

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

梯形积分留下的那一半

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

f_i(\mathbf{x}+\mathbf{c}_i\delta t,\ t+\delta t) - f_i(\mathbf{x},t) = -\frac{\delta t}{2\lambda}\left[(f_i-f_i^{\rm eq})^{n} + (f_i-f_i^{\rm eq})^{n+1}\right] + \frac{\delta t}{2}\left[F_i^{\,n} + F_i^{\,n+1}\right]

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

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

\bar f_i = f_i + \frac{\delta t}{2\lambda}\left(f_i - f_i^{\rm eq}\right) - \frac{\delta t}{2}F_i

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

\bar f_i(\mathbf{x}+\mathbf{c}_i\delta t,\ t+\delta t) = \bar f_i - \frac{1}{\tau}\left(\bar f_i - f_i^{\rm eq}\right) + \delta t\left(1-\frac{1}{2\tau}\right)F_i

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

\rho\mathbf{u} = \sum_i \bar f_i \mathbf{c}_i + \frac{\delta t}{2}\mathbf{F}

这两处是一体的。它们从同一次变换里一起掉出来。只顾上一处、漏掉另一处，账就对不上。差多少，下面来量。

四种格式对照表

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

格式	力放进去的位置	平衡态所用速度	二阶矩 $\sum\mathbf{c}_i\mathbf{c}_i F_i$
plain	$F_i = w_i(\mathbf{c}_i\cdot\mathbf{F})/c_s^2$	$\mathbf{j}/\rho$	0 — $\mathbf{u}\mathbf{F}$ 项整个缺失
Shan–Chen (1993)	没有显式项	$\mathbf{j}/\rho + \tau\mathbf{F}/\rho$	$\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u} + \tau\mathbf{F}\mathbf{F}/\rho$
He (1998)	$(\mathbf{c}_i-\mathbf{u})\cdot\mathbf{F}\,f_i^{\rm eq}/(\rho c_s^2)$	$(\mathbf{j}+\delta t\mathbf{F}/2)/\rho$	$\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u}$ （在平衡态截断误差之内）
EDM (2004)	$f_i^{\rm eq}(\rho,\mathbf{u}+\delta\mathbf{u}) - f_i^{\rm eq}(\rho,\mathbf{u})$	$\mathbf{j}/\rho$	$\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u} + \rho\,\delta\mathbf{u}\delta\mathbf{u}$
Guo (2002)	上面的 $F_i$ 再乘以 $(1-\frac{1}{2\tau})$	$(\mathbf{j}+\delta t\mathbf{F}/2)/\rho$	$(1-\frac{1}{2\tau})(\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u})$

把这五行竖着读，共同点就出来了。在乘上系数之前，力项的一阶矩全都精确等于 $\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}$ 相比），或者像多相流那样力在界面处急剧变化的位置。

在下面的图里直接动手试试。

τ1.00F angle θ0°|u|0.080

不管怎么转力方向的旋钮，M0（质量）柱条都贴在 0 上。力不制造质量，本该如此。M1（动量）的目标不是 $\mathbf{F}$ ，而是 $(1-\frac{1}{2\tau})\mathbf{F}$ 。因为量的是乘上系数之后的力项。缺掉的 $\mathbf{F}/(2\tau)$ 由平移半格的平衡态还回来，最下面一行 $\Delta j$ 正好落在 $\mathbf{F}$ 上，就是这件事的确认。

格式分道扬镳的只有 M2（应力）这一行。把 $|u|$ 滑块推上去，plain 的 M2 停在 0 上，够不着目标。错的是形状而不是大小，系数救不回来。把 $\tau$ 推向 2，Shan–Chen 的偏差会明显长大。一关掉 prefactor 开关，M1、M2 和 $\Delta j$ 会同时走样。

定常 Poiseuille 分辨不出它们#

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

u(y) = \frac{F_x}{2\nu}\,y\,(H-y),\qquad \nu = c_s^2\left(\tau-\tfrac{1}{2}\right)

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

在下面的模拟里直接动手试试。

τ1.00fx8.00e-6

把四个格式按钮全按一遍。虚线抛物线上的青绿色曲线纹丝不动。把 $\tau$ 从 0.6 推到 2.0 也一样。接着只关掉下面两个开关中的一个。就在那一刻，抛物线的幅值整个换了一档。

数字上也一样。取 33 个流体节点、 $F_x = 10^{-5}$ ，跑 40,000 步之后测得的最大相对误差如下。

$\tau$	plain	Shan–Chen	He	EDM	Guo
0.60	0.0875%	0.0875%	0.0875%	0.0875%	0.0875%
1.00	0.0306%	0.0306%	0.0306%	0.0306%	0.0306%
1.80	0.7358%	0.7358%	0.7358%	0.7358%	0.7358%

列与列之间没什么可比的。数值相同。剩下的误差不是来自格式，而是来自 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(\rho u) = \frac{1}{\tau}\cdot\frac{s\,\delta t F}{2} + g\,\delta t F

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

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

$\tau$	两个都开	只开修正	预测 $1+\frac{1}{2\tau}$	只开系数	预测 $1-\frac{1}{2\tau}$
0.60	0.9991	1.8316	1.8333	0.1666	0.1667
0.80	0.9995	1.6240	1.6250	0.3751	0.3750
1.00	1.0003	1.5002	1.5000	0.5005	0.5000
1.40	1.0030	1.3609	1.3571	0.6452	0.6429
1.80	1.0074	1.2867	1.2778	0.7280	0.7222

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

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

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

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

值得记住的点

力项前面的 $(1-\frac{1}{2\tau})$ 与速度上的 $+\delta t\mathbf{F}/2\rho$ ，是梯形积分的变量变换里一起出来的一对。只用其中一个，力就会偏成 $1\pm\frac{1}{2\tau}$ 倍。 $\tau=0.6$ 时超出 83%。
定常单向 Poiseuille 分辨不出 forcing 格式。五种格式给出的误差到第四位都相同。格式的比较得在非定常流动或非均匀力下做。
无论哪种格式，力项本身的一阶矩都被调成了 $\mathbf{F}$ 。差别在二阶矩，也就是应力上，而且目标值里含有 $(1-\frac{1}{2\tau})$ 。审视一个新格式时，先看 M2 最快。