τ=0.6时流速只有六分之一 — Guo forcing漏掉半步速度F/2的后果

施加同样的力,为什么流速只有六分之一

我在周期边界的盒子里施加了正弦形的体积力。解析解在教科书里只有一行。然而格子玻尔兹曼方法(LBM,通过碰撞和迁移分布函数来求解流体的方法)程序给出的振幅只有解析解的0.16倍。改变黏性,这个比值也跟着变。τ=1时恰好是一半。

原因只有一行。计算宏观速度时没有加上 $\mathbf{F}/2$。本文说明这一行为什么恰好造成 $(1-1/2\tau)$ 倍的误差。我们用Kolmogorov流(由正弦体积力驱动的周期剪切流)直接测量,最后给出一个30秒内在自己程序中找到这个错误的方法。

一个时间步里动量变化三次

LBM的一步是碰撞加迁移。迁移只搬运动量。动量只在碰撞中产生。因此只需看一个均匀流体节点的账本。

$\mathbf{m}=\sum_i f_i\mathbf{c}_i$ 是碰撞前动量,$\mathbf{m}'$ 是碰撞后动量,$\rho\mathbf{u}$ 是放进平衡分布的动量,$\tau$ 是松弛时间,$\mathbf{F}$ 是格子单位的体积力。

第二项是松弛项。若平衡速度等于碰撞前速度,它为零。第三项是Guo源项的一阶矩。往 $\rho\mathbf{u}$ 里放什么,决定了一切。

第一行中,松弛项补上了 $\mathbf{F}/(2\tau)$,恰好注入 $\mathbf{F}$。第二行中这份补偿消失了。只剩源项的 $(1-1/2\tau)$ 系数。

在下面的模拟中亲手调一调。

把τ滑块拖向0.51。红线的斜率几乎变为零,右侧柱状图中深色的松弛部分整块消失。

半步速度:碰撞前后的正中间

为什么偏偏是 $\mathbf{F}/2$?碰撞前动量是 $\mathbf{m}$,完全受力之后是 $\mathbf{m}+\mathbf{F}$。物理速度是这段时间的平均,也就是中点。这与梯形公式(用区间两端值的平均做积分,二阶精度)的道理相同。

Guo等(2002)把源项写成

$w_i$ 是格子权重,$\mathbf{c}_i$ 是离散速度,$c_s^2=1/3$ 是格子声速的平方。它满足 $\sum_i S_i\mathbf{c}_i = \mathbf{F}$ 和 $\sum_i S_i = 0$。

$(1-1/2\tau)$ 系数来自Chapman–Enskog展开(按Knudsen数展开分布函数以导出宏观方程的方法)。离散时间的Taylor展开会在源项中留下 $-\tfrac{1}{2}S_i$ 形式的残余。这个残余与半步黏性修正同根。要消去这个系数,就必须在平衡速度中加入 $\mathbf{F}/2$,使一阶矩平衡。两者从一开始就是成对推导出来的。

代码:Kolmogorov流与三种速度定义#

把D2Q9格子做成 $1\times N$ 的周期柱。力为 $F_x = F_0\sin(ky)$,$k=2\pi/N$。定常解析解为

$\nu$ 是格子单位的运动黏度。没有壁面,所以不会混入反弹边界的误差。只测量力的处理方式。

运行三种情形。guo 在平衡态和输出中都加入 $\mathbf{F}/2$。nohalf 保留Guo源项,但速度中去掉 $\mathbf{F}/2$。rawout 计算正确,只在输出时用 $\sum f_i\mathbf{c}_i/\rho$。振幅取 $\sin(ky)$ 的傅里叶系数,再除以解析值。

import math
 
# D2Q9 lattice: velocities, weights, sound speed squared
C = [(0,0),(1,0),(0,1),(-1,0),(0,-1),(1,1),(-1,1),(-1,-1),(1,-1)]
W = [4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36]
CS2 = 1/3
 
def f_equil(rho, ux):
    # uy = 0 everywhere: the flow is u_x(y) only
    out = []
    for (cx, cy), wi in zip(C, W):
        cu = cx*ux
        out.append(wi*rho*(1 + cu/CS2 + 0.5*(cu/CS2)**2 - 0.5*ux*ux/CS2))
    return out
 
def guo_term(ux, Fx, tau):
    out = []
    for (cx, cy), wi in zip(C, W):
        cu = cx*ux
        out.append((1 - 0.5/tau)*wi*((cx - ux)/CS2 + cu*cx/CS2**2)*Fx)
    return out
 
def kolmogorov_amp(tau, mode, N=16, F0=1e-6):
    """mode 'guo'    : u = (sum f c + F/2)/rho in equilibrium and output
       mode 'nohalf' : Guo source, but u = sum f c / rho everywhere
       mode 'rawout' : correct update, but report sum f c / rho"""
    k = 2*math.pi/N
    Fx = [F0*math.sin(k*(j + 0.5)) for j in range(N)]
    nu = (tau - 0.5)/3
    steps = int(12/(nu*k*k))                    # ~12 viscous decay times
    f = [f_equil(1.0, 0.0) for _ in range(N)]
    for _ in range(steps):
        post = []
        for j in range(N):
            fj = f[j]
            rho = sum(fj)
            mx = sum(fi*cx for fi, (cx, _) in zip(fj, C))
            ux = mx/rho if mode == 'nohalf' else (mx + 0.5*Fx[j])/rho
            feq = f_equil(rho, ux)
            src = guo_term(ux, Fx[j], tau)
            post.append([fi - (fi - fe)/tau + s for fi, fe, s in zip(fj, feq, src)])
        # streaming: only the y component moves anything (one cell in x, periodic)
        f = [[post[(j - C[i][1]) % N][i] for i in range(9)] for j in range(N)]
    amp = 0.0
    for j in range(N):
        rho = sum(f[j]); mx = sum(fi*cx for fi, (cx, _) in zip(f[j], C))
        u = mx/rho if mode in ('nohalf', 'rawout') else (mx + 0.5*Fx[j])/rho
        amp += 2*u*math.sin(k*(j + 0.5))/N
    return amp/(F0/(nu*k*k)), steps
 
print(" tau  steps   guo     nohalf  rawout   nohalf/guo  1-1/(2tau)")
for tau in [0.55, 0.6, 0.8, 1.0, 1.5]:
    r1, n = kolmogorov_amp(tau, 'guo')
    r2, _ = kolmogorov_amp(tau, 'nohalf')
    r3, _ = kolmogorov_amp(tau, 'rawout')
    print(f"{tau:4.2f} {n:6d}  {r1:.4f}  {r2:.4f}  {r3:.4f}   {r2/r1:.4f}      {1-0.5/tau:.4f}")

 tau  steps   guo     nohalf  rawout   nohalf/guo  1-1/(2tau)
0.55   4668  0.9746  0.0885  0.9734   0.0908      0.0909
0.60   2334  0.9754  0.1621  0.9729   0.1662      0.1667
0.80    778  0.9836  0.3660  0.9759   0.3721      0.3750
1.00    466  1.0001  0.4936  0.9872   0.4936      0.5000
1.50    233  1.0772  0.7010  1.0515   0.6508      0.6667

把上面的代码直接放到浏览器里跑一跑。

τ设为0.6时,红色 nohalf 振幅停在0.16附近,贴着红色虚线 $1-1/2\tau$。把N调到32以上,绿色 guo 更接近1,而 nohalf/guo 比值不变。

0.166和0.973说明了什么#

nohalf/guo 一列与 $1-1/2\tau$ 重合:τ=0.55时0.0908对0.0909,τ=0.6时0.1662对0.1667。动量账本的结论在真实流动中原样出现。定常状态下,黏性扩散与每步注入的动量相平衡。注入量是 $(1-1/2\tau)$ 倍,振幅也就是这个倍数。

τ越大,比值越偏离理论值:τ=1.5时0.6508对0.6667。$N=16$ 时 $k\approx0.39$,每个波长的格点很少。这里会混入 $k^2$ 阶的格子截断误差。guo 本身为1.0772而不是1,也是同样的原因。

rawout 与 guo 几乎相同。输出速度只是滞后 $F/(2\rho)$。相对误差约为 $\nu k^2/2$,τ=0.55时约0.1%。这是调试时容易混淆的地方。输出端漏掉 $\mathbf{F}/2$ 误差很小,平衡端漏掉误差很大。

最危险的是高Reynolds数计算。为降低黏性把τ压向0.5,有效力也随之趋于零。若因为流速上不去就加大力,隐藏的系数会残留在所有结果中。

F/2 应放在伪代码的哪里#

常见的LBM循环顺序是:初始化 → 计算宏观量 → 平衡分布 → 力 → 碰撞 → 迁移。错误正是这个顺序造成的。在宏观量步骤求 $\mathbf{u}$ 时,力还没有算出来。于是很容易停在 $\mathbf{u}=\sum f_i\mathbf{c}_i/\rho$。

解决办法是稍微调整顺序。

1. 先求 $\rho = \sum_i f_i$。
2. 计算力 $\mathbf{F}$。若像重力那样只需要 $\rho$,这一步就完成了。
3. 用 $\rho\mathbf{u} = \sum_i f_i\mathbf{c}_i + \mathbf{F}/2$ 求速度。
4. 用同一个 $\mathbf{u}$ 构造 $f^{eq}$ 和 $S_i$。输出也用这个 $\mathbf{u}$。

若力依赖于速度(例如拖曳力),第3步就成为隐式方程。线性拖曳力可以代数求解。其他情况做一两次不动点迭代就够了。

删掉 $(1-1/2\tau)$ 来"修正"也不是办法。一阶矩对了,但二阶矩里会留下依赖τ的应力。这个问题在力项三阶矩一文和Shan–Chen速度平移forcing一文中讨论过。

在自己的LBM循环里找F/2的单步测试#

建一个周期盒子,填充 $\rho=1$、$\mathbf{u}=0$。施加均匀力 $\mathbf{F}$,只运行一步。然后用 $\sum f_i\mathbf{c}_i$ 的变化量除以 $\mathbf{F}$。

结果为1说明配对一致。若为 $1-1/2\tau$,说明平衡速度中漏了 $\mathbf{F}/2$。取两个τ值各跑一次更有把握。比值随τ变化就是这个错误。运行任何验证算例之前,先检查这一步。