只改了 τ，静止气体的密度就变了 — Shan–Chen 速度平移 forcing 的 (τ−½)² 应力

把平衡速度平移 τF/ρ 的 forcing 会在动量通量里加上 (τ−½)²F²/ρ，因此静止状态的密度分布随松弛时间改变。

只改了黏性，静止状态却变了

静止的气体并不知道自己的黏性。速度为零，黏性应力也为零。所以力与压力平衡的静止状态应当与黏性无关。在格子 Boltzmann 方法（LB, Lattice Boltzmann）中，黏性由松弛时间 $\tau$ 决定。因此改变 $\tau$ ，静止状态的密度分布应当不变。

但有一类 forcing 做不到这一点。这就是 Shan–Chen 类多相模型使用的速度平移（velocity shift）forcing。同样的力、同样的格子，只把 $\tau$ 从 0.8 改到 3，密度分布就会移动。本文把三种 forcing 并排比较，找出差别的来源。读完后你会得到一行结论：每种格式悄悄加进动量通量的应力系数。Guo 是 0，Shan–Chen 是 $(\tau-\tfrac12)^2$ ，EDM 是 $\tfrac14$ 。

比较标准一：平衡速度被移到哪里

三种格式每一步都加入动量 $F\,\delta t$ 。区别在于把这份动量放进哪个分布。下面 $u'=\sum_i f_i c_i/\rho$ 是直接从分布读出的速度。三种格式的物理速度都定义为 $u = u' + F/(2\rho)$ （半步修正）。

格式	平衡中使用的速度	附加源项
Guo	$u$	$(1-\tfrac{1}{2\tau})S_i(u,F)$
Shan–Chen	$u' + \tau F/\rho$	无
EDM (Kupershtokh)	$u'$	$f_i^{eq}(u'+F/\rho)-f_i^{eq}(u')$

静止时 $u=0$ ，所以 $u'=-F/(2\rho)$ 。于是 Shan–Chen 放进平衡的速度为

u_{eq} = u' + \frac{\tau F}{\rho} = \left(\tau - \tfrac12\right)\frac{F}{\rho}

$u_{eq}$ 是构造平衡分布所用的速度， $\tau$ 是无量纲松弛时间， $F$ 是单位体积力。气体静止不动，平衡却以正比于 $\tau$ 的速度在运动。

平衡分布的二阶矩（动量通量）是速度的抛物线：

\Pi_{eq}(u_{eq}) = \sum_i f_i^{eq} c_i c_i = c_s^2\rho + \rho\, u_{eq}^2

$c_s^2=1/3$ 是格子声速的平方。 $u_{eq}$ 一旦偏离零，通量就被抬高 $\rho u_{eq}^2$ 。在下面的模拟中亲手操作一下。光标沿气柱移动，用该处的力 $F(y)$ 把各格式的平衡点标在抛物线上。

τ = 1.50 force amplitude G = 0.120

y = 0.0 · F = 0.0000 · SC u_eq = (τ−½)F/ρ = 0.0000 · SC extra 0.00e+0 · EDM extra 0.00e+0 · Guo 0

把 $\tau$ 滑块从 0.5 附近拉到 3。橙色点（Shan–Chen）沿抛物线上爬，绿色点（Guo）停在底部。紫色点（EDM）位于 $\pm F/(2\rho)$ 两处，高度 $F^2/(4\rho)$ 与 $\tau$ 无关。

比较标准二：这份抬高能否保留到静止状态

平衡被抬高，并不意味着宏观应力立刻被抬高。真正输运动量的是碰撞后分布的二阶矩 $\Pi^*$ 。在定常状态下令碰撞前后的二阶矩相等，得到

\Pi^* = \Pi_{eq} + \tau\, S_\Pi

$S_\Pi$ 是附加源项的二阶矩。Guo 的源项在 $u=0$ 时二阶矩为零。EDM 的源项为 $\rho(u'+F/\rho)^2-\rho u'^2 = 0$ 。所以三种格式的差别只来自放进平衡的速度。静止气柱的力平衡变为

\frac{d}{dy}\left(c_s^2\rho + b\,\frac{F^2}{\rho}\right) = F, \qquad b_{Guo}=0,\quad b_{SC}=\left(\tau-\tfrac12\right)^2,\quad b_{EDM}=\tfrac14

$b$ 是虚假应力的系数。连续介质的静水压平衡对应 $b=0$ 。

为什么教科书推导里看不到这一项？标准的 Chapman–Enskog 展开（用小参数 $\epsilon$ 把分布按阶分解的多尺度展开）把 $F$ 当作一阶小量。于是 $F^2/\rho$ 是 $\epsilon^2$ 项，在 Navier–Stokes 层次被舍弃。可一旦系数按 $\tau^2$ 增长，"小"就没有保证了。源笔记引用的 Revised Chapman–Enskog analysis for a class of forcing schemes in the lattice Boltzmann method 正是重新审视力项阶数分配的论文。力该以压力梯度还是体积力的形式加入，已在非理想 LBM forcing 的两种形式中讨论过。本文关注的是同一个力如何加载到分布上。

代码：周期气柱、三种 forcing、五个 τ#

在 32 格的 D1Q3（速度 0, ±1）周期气柱上施加 $F = G\sin(ky)$ 。初始密度均匀。声波被黏性耗散后，剩下静止状态。对收敛后的密度用最小二乘拟合 $c_s^2\rho = \lambda P(y) - bF^2/\rho + c$ ，测出 $b$ 。 $P(y)=-\tfrac{G}{k}\cos(ky)$ 是满足 $P'=F$ 的静水压势。 $\lambda$ 吸收格子离散对振幅造成的倍率。

import math
 
# D1Q3 格子: 速度 0, +1, -1 ; cs^2 = 1/3
C = (0.0, 1.0, -1.0)
W = (2/3, 1/6, 1/6)
CS2 = 1/3
 
def d1q3_eq(rho, u):
    return [W[i]*rho*(1 + C[i]*u/CS2 + 0.5*(C[i]*u)**2/CS2**2 - 0.5*u*u/CS2)
            for i in range(3)]
 
def collide_node(fn, F, tau, scheme):
    rho = sum(fn)
    mom = fn[1] - fn[2]
    if scheme == "guo":                 # Guo: u = (sum f c + F/2)/rho, 源项乘以 (1 - 1/2tau)
        u = (mom + 0.5*F)/rho
        feq = d1q3_eq(rho, u)
        out = []
        for i in range(3):
            S = W[i]*((C[i] - u)/CS2 + C[i]*u*C[i]/CS2**2)*F
            out.append(fn[i] - (fn[i] - feq[i])/tau + (1 - 0.5/tau)*S)
        return out
    if scheme == "sc":                  # Shan-Chen: 平衡速度平移 tau F / rho
        feq = d1q3_eq(rho, mom/rho + tau*F/rho)
        return [fn[i] - (fn[i] - feq[i])/tau for i in range(3)]
    up = mom/rho                        # EDM: 加上 feq(u'+F/rho) - feq(u')
    feq = d1q3_eq(rho, up)
    fsh = d1q3_eq(rho, up + F/rho)
    return [fn[i] - (fn[i] - feq[i])/tau + fsh[i] - feq[i] for i in range(3)]
 
def run_column(scheme, tau, G, N=32, tol=1e-12):
    k = 2*math.pi/N
    F = [G*math.sin(k*(j + 0.5)) for j in range(N)]
    f = [d1q3_eq(1.0, 0.0) for _ in range(N)]
    for it in range(1, 100001):
        post = [collide_node(f[j], F[j], tau, scheme) for j in range(N)]
        new = [[post[j][0], post[j-1][1], post[(j+1) % N][2]] for j in range(N)]
        if it % 100 == 0 and max(abs(sum(new[j]) - sum(f[j])) for j in range(N)) < tol:
            return F, [sum(fn) for fn in new], it
        f = new
    return F, [sum(fn) for fn in f], it
 
def stress_fit(F, rho, G):
    # 拟合  cs^2 rho = lam*P(y) - b*F^2/rho + c,  P' = F  (静水压时 b = 0)
    N = len(F); k = 2*math.pi/N
    X = [[F[j]**2/rho[j], -G/k*math.cos(k*(j + 0.5)), 1.0] for j in range(N)]
    r = [CS2*x for x in rho]
    A = [[sum(X[n][a]*X[n][b] for n in range(N)) for b in range(3)] for a in range(3)]
    v = [sum(X[n][a]*r[n] for n in range(N)) for a in range(3)]
    for i in range(3):
        for j in range(i + 1, 3):
            m = A[j][i]/A[i][i]
            A[j] = [A[j][c] - m*A[i][c] for c in range(3)]
            v[j] -= m*v[i]
    x = [0.0]*3
    for i in (2, 1, 0):
        x[i] = (v[i] - sum(A[i][c]*x[c] for c in range(i + 1, 3)))/A[i][i]
    return -x[0], x[1]
 
G = 0.02
print(f"periodic column N=32, F = {G} sin(ky), rho0 = 1")
print(" tau | b_Guo  b_SC   (tau-1/2)^2  b_EDM | lam   | max|rho_SC-rho_Guo|  iters")
for tau in [0.8, 1.0, 1.5, 2.0, 3.0]:
    out = {s: run_column(s, tau, G) for s in ("guo", "sc", "edm")}
    b = {s: stress_fit(out[s][0], out[s][1], G) for s in out}
    drho = max(abs(p - q) for p, q in zip(out["sc"][1], out["guo"][1]))
    bg = round(b["guo"][0], 3) + 0.0
    print(f"{tau:4.1f} | {bg:5.3f} {b['sc'][0]:6.3f}   {(tau-0.5)**2:6.3f}    {b['edm'][0]:5.3f} |"
          f" {b['guo'][1]:.4f} |  {drho:.2e}            {out['sc'][2]}")

periodic column N=32, F = 0.02 sin(ky), rho0 = 1
 tau | b_Guo  b_SC   (tau-1/2)^2  b_EDM | lam   | max|rho_SC-rho_Guo|  iters
 0.8 | 0.000  0.090    0.090    0.250 | 0.9968 |  5.49e-05            6100
 1.0 | 0.000  0.250    0.250    0.250 | 0.9968 |  1.53e-04            3600
 1.5 | 0.000  1.000    1.000    0.250 | 0.9968 |  6.11e-04            1900
 2.0 | 0.000  2.250    2.250    0.250 | 0.9968 |  1.38e-03            1300
 3.0 | 0.000  6.250    6.250    0.250 | 0.9968 |  3.84e-03            800

0、(τ−½)² 和 1/4 说明了什么#

在浏览器里直接运行上面的代码。三根气柱在同样的力下同时推进，每一帧重新拟合 $b$ 。

τ = 1.50 G = 0.020 steps/frame 10

step 0 · b: Guo 0.000 · SC 0.000 (theory 1.000) · EDM 0.000 · max|ρ_SC−ρ_Guo| 0.00e+0 · Δρ/step 1.0e+0

默认 $\tau=1.5$ 、 $G=0.02$ 时，橙色曲线贴在 1.000，紫色 0.250，绿色 0.000，与输出块第三行一致。把 $\tau$ 调到 3 再按 Reset，只有橙色升到 6.25。改变 $G$ 时点的幅度按 $G^2$ 变化，但 $b$ 不变。

输出中首先要看 $b_{SC}$ 列和 $(\tau-\tfrac12)^2$ 列。五行全部吻合到小数点后三位。说明推导中没有遗漏。Guo 的 0.000 同样重要：静水压平衡一直成立到 $F^2$ 量级。

$\tau=1$ 时 Shan–Chen 和 EDM 同为 0.250。这并非巧合。 $\tau=1$ 时碰撞后的分布就是平衡分布。Shan–Chen 变为 $f^{eq}(u'+F/\rho)$ ，EDM 也变为 $f^{eq}(u') + f^{eq}(u'+F/\rho) - f^{eq}(u')$ ，两者相同。两种格式在 $\tau=1$ 处交于一点。

$\lambda=0.9968$ 与 forcing 无关。它是 32 格格子对正弦振幅的离散倍率，数值与 $1-k^2/12$ 一致。加密网格后趋于 1。拟合的设计使 $b$ 与这个倍率分离。

密度差本身看上去很小。 $\tau=3$ 时为 $3.84\times10^{-3}$ ，是密度振幅 0.30 的 1.3%。但这个比例随 $G$ 成正比增大。主效应正比于 $F$ ，而虚假应力正比于 $F^2$ 。

一张表看三种格式

	Guo	Shan–Chen	EDM
平衡速度（静止时）	$0$	$(\tau-\tfrac12)F/\rho$	$-F/(2\rho)$
虚假应力系数 $b$	$0$	$(\tau-\tfrac12)^2$	$\tfrac14$
对 $\tau$ 的依赖	无	二次	无
静止状态与黏性无关	是	否	是
符合静水压平衡	是	仅在 $\tau=\tfrac12$ 附近	否（固定 $\tfrac14$ ）

EDM 与 $\tau$ 无关，但偏差固定为 $\tfrac14$ 。偏差是常数，修正起来容易。Shan–Chen 只在 $\tau=0.5$ 附近才小，而那正是 LB 最不稳定的区域。

在多相界面上，F 很大#

本实验中的力光滑且弱。伪势（pseudopotential）多相模型的力不同。相互作用力 $F=-\mathcal{G}\psi\nabla\psi$ （ $\mathcal{G}$ 为相互作用强度， $\psi$ 为有效密度）集中在界面附近几格，而且很大。在那里， $bF^2/\rho$ 可以增长到与状态方程给出的压力同一量级。

因此用 Shan–Chen forcing 建立平直界面时，共存密度（coexistence density，两相平衡时各相的密度）会随 $\tau$ 移动。本想调 $\tau$ 改变黏性，液相密度却变了。这个气柱实验把该现象在没有界面的情况下单独分离出来。虚假应力会原样进入界面的力学平衡条件。

反过来，也有方法有意利用这一项。把 $b$ 当作可调参数，它就成了让共存密度满足热力学一致性的旋钮：在 EDM 的 $\tfrac14$ 或 Guo 的 0 之上按需叠加。无论哪种情况，都得先知道自己代码里的 $b$ 是多少。用同一套 Chapman–Enskog 框架找出热模型隐藏项的例子，见 2-population LB 的准平衡 g*。

引入新 forcing 前先跑的 32 格#

往代码里加入新的 forcing 格式时，不要先跑多相问题。在界面上各种原因会混在一起。先跑上面的周期气柱。用两个 $\tau$ 值测量 $b$ 。如果 $b$ 随 $\tau$ 变化，说明平衡里藏着速度平移。如果 $b$ 是非零常数，这个常数就是将叠加到界面压力上的量。32 格、几千步即可完成的检查。