τ=0.51时把黏性读大了51倍 — 离散LBE二阶展开留下的半个时间步

格子LBE的黏性响应的是τ−1/2而不是τ。那半个时间步是离散展开的二阶项留下的。

连续方程给出的黏性和格子返回的黏性相差51倍#

只看连续BGK方程，缓和时间 $\tau$ 与运动黏性系数 $\nu$ 的关系是 $\nu = c_s^2 \tau$ 。把 $\tau$ 取成0.51， $\nu$ 应当是0.17。可是把同一个 $\tau$ 送进D2Q9代码，测量Taylor–Green涡的衰减率，得到的是0.00333。差了51倍。

这不是bug，也不是符号弄反。是缓和时间正好少了 $\Delta t/2$ 。半步从两条路走出来：特征线上的二阶Taylor展开，以及碰撞积分的中点法则。两条路在同一处汇合。把 $\tau$ 滑块推向0.51，橙色虚线 $c_s^2\tau$ 就离开格子上的振幅。

tau = 0.600  run

The blue curve is the vortex actually living on the lattice at the left. Green is cs2(tau-1/2) = 3.333e-2, orange is cs2*tau = 2.000e-1 — 1.0x apart. Fitted from the blue curve: nu = 0.000e+0. Push tau toward 0.51 and the orange line stops resembling anything on screen.

蓝色曲线是32×32格子上还活着的涡的振幅。绿色虚线是 $c_s^2(\tau - 1/2)$ 。 $\tau = 0.51$ 时绿线在0.00333，橙线在0.17。

格子上的那一行不是微分方程

我们真正在跑的是这一行。

f_i(\mathbf{x} + \mathbf{c}_i \Delta t,\, t + \Delta t) - f_i(\mathbf{x}, t) = -\frac{\Delta t}{\tau}\,\Omega_i(\mathbf{x}, t)

其中 $f_i$ 是沿离散速度 $\mathbf{c}_i$ 的分布函数， $\Omega_i = f_i - f_i^{\rm eq}$ 是非平衡部分（偏离平衡的量）， $\Delta t$ 是一个时间步。

左边不是导数。是跨越有限间隔 $\Delta t$ 实际跳过去的差。要把它化成导数就得展开，而展开里有二阶项。CFD里常见的做法是只留一阶项，把其余的叫作"截断误差"。这里那个二阶项不是误差，它改写了一个物理系数。

数一下量级就知道为什么不能丢。二阶项的大小是 $\Delta t^2 D^2 f / 2$ ，而右边的碰撞项被 $\Delta t / \tau$ 除过。在格子单位下 $\tau$ 是0.5到2之间的 $O(1)$ 数， $\Delta t = 1$ 。也就是说二阶项与一阶项的比值不会变小。把本该趋于无穷小的 $\Delta t$ 钉在1上，正是格子气方法的起点，而这一项就是代价。

沿特征线展开到二阶会剩下什么

记 $D_i \equiv \partial_t + \mathbf{c}_i \cdot \nabla$ ，左边就成了 $D_i$ 的指数函数。保留到二阶并除以 $\Delta t$ 。

D_i f_i + \frac{\Delta t}{2} D_i^2 f_i + O(\Delta t^2) = -\frac{1}{\tau} \Omega_i

问题就在 $D_i^2$ 这一项。丢掉它得到 $\nu = c_s^2\tau$ ，留下它得到另一个值。

做一次Chapman–Enskog展开。取 $f_i = f_i^{(0)} + \varepsilon f_i^{(1)} + \varepsilon^2 f_i^{(2)}$ 、 $\partial_t = \varepsilon \partial_{t_1} + \varepsilon^2 \partial_{t_2}$ 、 $\nabla = \varepsilon \nabla_1$ ，且 $f_i^{(0)} = f_i^{\rm eq}$ 。在 $\varepsilon$ 阶得到

D_{i1} f_i^{(0)} = -\frac{1}{\tau} f_i^{(1)}

把这个关系原样搬进 $\varepsilon^2$ 阶，是整件事的关键。因为 $D_{i1}^2 f_i^{(0)} = -\tfrac{1}{\tau} D_{i1} f_i^{(1)}$ ，于是

\partial_{t_2} f_i^{(0)} + \left(1 - \frac{\Delta t}{2\tau}\right) D_{i1} f_i^{(1)} = -\frac{1}{\tau} f_i^{(2)}

$D_{i1} f_i^{(1)}$ 前面多了一个括号。黏性应力由 $f_i^{(1)}$ 的二阶矩乘 $\tau$ 得到，所以那个 $\tau$ 变成了 $\tau\left(1 - \frac{\Delta t}{2\tau}\right) = \tau - \frac{\Delta t}{2}$ 。

\nu = c_s^2 \left(\tau - \frac{\Delta t}{2}\right)

在格子单位 $\Delta t = 1$ 下就是 $c_s^2(\tau - 1/2)$ 。因为 $c_s^2 = 1/3$ ， $\tau = 0.51$ 给出 $\nu = 0.00333$ 。上面模拟里的绿色虚线正是这个数。

把剩下的项推进分布函数里

还有一条不用展开就能到达同一答案的路。沿特征线积分连续BGK方程，但采用梯形法则。把缓和时间记作 $\lambda$ 。

f_i(\mathbf{x} + \mathbf{c}_i\Delta t, t + \Delta t) - f_i(\mathbf{x},t) = -\frac{\Delta t}{2\lambda}\left[\Omega_i\big|_{\mathbf{x},t} + \Omega_i\big|_{\mathbf{x}+\mathbf{c}_i\Delta t,\, t+\Delta t}\right]

右边含有未来时刻，所以这是隐式的。就这样是跑不起来的。于是换变量。

\bar f_i = f_i + \frac{\Delta t}{2\lambda}\,\Omega_i

用 $\bar f_i$ 重写，未来时刻的项相互抵消而消失。代入 $\bar\Omega_i = \bar f_i - f_i^{\rm eq} = \left(1 + \frac{\Delta t}{2\lambda}\right)\Omega_i$ ，剩下的就是这个。

\bar f_i(\mathbf{x}+\mathbf{c}_i\Delta t, t+\Delta t) - \bar f_i(\mathbf{x},t) = -\frac{\Delta t}{\lambda + \Delta t/2}\,\bar\Omega_i(\mathbf{x},t)

完全显式，而且和我们在用的那一行连字符都一样——但只有在 $\tau = \lambda + \Delta t/2$ 时成立。也就是说代码里的 $\tau$ 不是物理缓和时间，物理的那个是 $\lambda = \tau - \Delta t/2$ 。这不是舍入问题，而是变量替换开出的账单。梯形积分把这半步留下来的地方，接在缓和时间上的 $\Delta t/2$ 。

把同一行当作中点法则重读

梯形法则是 $\theta$ 加权法则中 $\theta = 1/2$ 的情形。扫一遍 $\theta$ 就能看出为什么偏偏是一半。单步放大率是 $A_\theta(z) = \frac{1 - (1-\theta)z}{1 + \theta z}$ ，正确答案是 $e^{-z}$ ，其中 $z = \Delta t/\lambda$ 。

$\theta$ 一旦离开1/2，单步放大率的误差就按 $(\theta - \tfrac12)z^2$ 走，落在斜率2上。只有 $\theta = 1/2$ 时才掉到斜率3。这个一半不是可调参数，而是与 $e^{-z}$ 二阶系数 $1/2$ 咬合的唯一取值。

用六个τ值让Taylor–Green涡衰减#

论证到此为止。接下来直接问格子。二维Taylor–Green涡是量 $\nu$ 的尺子，因为它的振幅按 $e^{-2\nu k^2 t}$ 下降。

import numpy as np
 
C = np.array([[0,0],[1,0],[0,1],[-1,0],[0,-1],[1,1],[-1,1],[-1,-1],[1,-1]])
W = np.array([4/9] + [1/9]*4 + [1/36]*4)
CS2 = 1/3
 
def equilibrium_d2q9(rho, ux, uy):
    cu = C[:,0,None,None]*ux + C[:,1,None,None]*uy
    u2 = ux*ux + uy*uy
    return W[:,None,None]*rho*(1 + cu/CS2 + cu*cu/(2*CS2**2) - u2/(2*CS2))
 
def taylor_green_init(n, u0):
    k = 2*np.pi/n
    x, y = np.meshgrid(np.arange(n), np.arange(n), indexing='ij')
    ux = -u0*np.cos(k*x)*np.sin(k*y)
    uy =  u0*np.sin(k*x)*np.cos(k*y)
    return np.ones((n,n)), ux, uy, k
 
def advance_lattice(f, tau):
    rho = f.sum(axis=0)
    ux = (C[:,0,None,None]*f).sum(axis=0)/rho
    uy = (C[:,1,None,None]*f).sum(axis=0)/rho
    f = f - (f - equilibrium_d2q9(rho, ux, uy))/tau      # 碰撞
    for i in range(9):                                    # 迁移
        f[i] = np.roll(f[i], (C[i,0], C[i,1]), axis=(0,1))
    return f, ux, uy
 
def measure_decay_rate(n, tau, u0=0.01, nstep=400, skip=100):
    rho, ux, uy, k = taylor_green_init(n, u0)
    f = equilibrium_d2q9(rho, ux, uy)
    amp = []
    for s in range(nstep):
        f, ux, uy = advance_lattice(f, tau)
        amp.append(np.sqrt((ux**2 + uy**2).mean()))
    t = np.arange(skip, nstep, dtype=float)
    a = np.log(np.array(amp[skip:]))                      # 丢弃初始过渡段
    slope = np.polyfit(t, a, 1)[0]
    return -slope/(2*k*k), k
 
def viscosity_scan(n=32):
    print(f"{'tau':>6} {'nu_measured':>12} {'cs2*(tau-1/2)':>14} {'cs2*tau':>10} "
          f"{'err_half[%]':>11} {'err_naive[%]':>12}")
    for tau in (0.51, 0.55, 0.6, 0.8, 1.0, 1.5):
        nu, k = measure_decay_rate(n, tau)
        half = CS2*(tau - 0.5)
        naive = CS2*tau
        print(f"{tau:6.2f} {nu:12.6f} {half:14.6f} {naive:10.6f} "
              f"{100*abs(nu-half)/half:11.2f} {100*abs(nu-naive)/naive:12.2f}")
 
viscosity_scan()

   tau  nu_measured  cs2*(tau-1/2)    cs2*tau err_half[%] err_naive[%]
  0.51     0.003333       0.003333   0.170000        0.01        98.04
  0.55     0.016667       0.016667   0.183333        0.00        90.91
  0.60     0.033337       0.033333   0.200000        0.01        83.33
  0.80     0.100115       0.100000   0.266667        0.12        62.46
  1.00     0.167205       0.166667   0.333333        0.32        49.84
  1.50     0.337688       0.333333   0.500000        1.31        32.46

把 theta 从1/2拉开，右边的斜率就贴在2上，只有1/2才落到3。

theta = 0.500 dt/tau = 0.800  run

Only theta = 1/2 puts the pink line where the blue area demands, so its error rides slope 3 instead of slope 2 — that single half is the whole difference. Realised relaxation: tau_eff/tau = 0.9442. The lattice Boltzmann update you write is theta = 1/2 in f-bar, which is why the viscosity answers to tau - 1/2 and not tau. Cursor sits at n = 0.

左边蓝色区域是 $e^{-s/\lambda}$ 的面积，粉色虚线是 $\theta$ 法则用来代替它的那一个高度。只有 $\theta = 1/2$ 与 $e^{-z}$ 的二阶系数 $1/2$ 咬合。表里半步误差0.01%就是这个中点对上时剩下的残差。

$\tau - 1/2$ 一侧的误差是0.01%到1.31%，正是格子离散化留下的高阶误差水平。 $\tau$ 一侧是32%到98%，这个量级已经不能叫误差了。

初始振幅取0.01是为了抹掉压缩性。马赫数在0.02量级时， $f^{\rm eq}$ 被截断的三阶项不会混进衰减率。丢掉前100步也是同样的理由。用平衡分布给初始条件时 $f^{(1)}$ 还是0，非平衡成分需要几十步才安定下来。把那一段算进去再拟合斜率，黏性会比实际小。

$\tau$ 越大 $\tau - 1/2$ 一侧的误差越往上爬，也值得一读。到 $\tau = 1.5$ 时是1.31%，因为展开中丢掉的 $O(\Delta t^2)$ 项随 $\tau$ 增长。半步修正不是终点，它只是二阶为止的故事。

把表按 $\tau/(\tau - 1/2)$ 再读一遍，51倍就是一行。 $\tau = 0.51$ 时 $0.51/0.01 = 51$ ， $0.170000/0.003333$ 也是51.0。 $\tau = 0.55$ 时 $0.55/0.05 = 11$ ， $0.183333/0.016667$ 也是11。 $\tau = 0.60$ 是6， $\tau = 1.0$ 是2， $\tau = 1.5$ 是1.5。朴素误差98.04%、90.91%、83.33%、62.46%、49.84%、32.46%就是这个比走向1时留下的痕迹。

半步那一列安静得多。 $\tau = 0.51$ 时测量和预测都是0.003333，误差0.01%。 $\tau = 0.80$ 测量0.100115，预测0.100000，差0.000115，也就是0.12%。 $\tau = 1.50$ 是0.337688对0.333333，差0.004355，也就是1.31%。丢掉的 $O(\Delta t^2)$ 随 $\tau$ 变大，就写在这三格里。初始振幅0.01、丢掉前100步、格子32，是为了把这份残差从可压缩性里剥出来。马赫数0.02时，被截掉的三阶项进不了衰减率。 $\tau = 0.55$ 的朴素误差90.91%等于 $1 - 0.016667/0.183333$ 。那个比是1/11，按连续式 $\nu = c_s^2\tau$ 去读，等于在代码里写进十一倍的黏性。即便 $\tau = 1.5$ ，比仍是1.5，朴素误差32.46%算不上小。

边界条件和体力项再次开出同一张半步账单

这半个时间步不止出现一次。

bounce-back边界把壁面放在两个节点中间而不是节点上，原因相同：分布函数反射的时刻在时间步的中间。 $\tau$ 一旦离开0.5，这个有效壁面位置就随 $\tau$ 漂移，这正是著名的二阶精度失效。

Guo型体力项带一个 $\left(1 - \frac{1}{2\tau}\right)$ 因子，原因也相同。体力沿同一条特征线积分，受同一个梯形法则，而 $\bar f$ 替换只作用在 $\Omega_i$ 上，于是留下一个系数。漏掉这个括号，本该静止的流体会缓慢流动，或者重力下的静水压分布会略微倾斜。

把格子间距减半时不能简单地把 $\tau$ 乘2，根源也在这里。要保持 $\nu$ 不变，得解 $\tau_{\rm fine} - 1/2 = 2(\tau_{\rm coarse} - 1/2)$ 。细化循环里这一行放在哪里，会把误差分开。

当 $\tau$ 逐节点变化时，下界就是这半个时间步。因为 $\tau \to 0.5$ 时 $\nu \to 0$ ，黏性在局部趋近于零的问题都会撞上格子本身砌的墙。

0.5不是下界，而是坐标原点#

把 $\tau > 0.5$ 当成稳定性条件背下来，会漏掉原因。物理缓和时间是 $\tau - 1/2$ 。 $\tau = 0.51$ 时 $c_s^2\tau = 0.170000$ ， $c_s^2(\tau - 1/2) = 0.003333$ ，比是51。负的缓和时间不是物理。

该写进代码的是 $\nu$ ，不是 $\tau$ 。把 $\tau = 1/2 + \nu/c_s^2$ 放在输入端一处。体力项上的 $\left(1 - \frac{1}{2\tau}\right)$ 也是同一张半步账单。把这个括号和力的三阶矩放在一起时的界限，接在力项必须守住的那一条线。衰减若慢了51倍，就是有人把 $\tau$ 当成 $\nu$ 在读。