速度调高后扩散少了27% — LBM对流扩散模型中多出来的通量

格子玻尔兹曼里，扩散系数由一个松弛时间定死。$D = c_s^2(\tau - 1/2)\delta t$。一行。没有速度。流动打开之后，还会是同一个 $D$ 吗。均匀对流速度 $u$ 下，它降到 $D(1 - u^2/c_s^2)$。$u = 0.3$ 时少掉 27%。Chapman–Enskog 展开把这一项漏在哪一步，一个源项又怎么把它找回来，按展开的顺序往下对。

疼的是 phase-field 多相流。用格子玻尔兹曼解 Cahn–Hilliard 或 Allen–Cahn，界面厚度就被迁移率拴住。迁移率差 27%，厚度就差 27%，表面张力系数跟着偏。

扩散系数标定好了，界面却越来越薄

先用眼睛看。下面是一维对流扩散方程

用 D1Q3 格子玻尔兹曼求解的结果。初始条件是一个高斯分布，精确解也是高斯分布，宽度按 $\sigma^2(t) = \sigma_0^2 + 2Dt$ 增长。把 u 从 0 推到 0.30，峰比虚线更窄，tau 怎么转，这个比例都不动。

u = 0.30 时实线（计算）比虚线（精确解）更尖。下方 $\sigma^2$ 够不到虚线斜率 $2D$。tau 放在 0.6 或 1.5，亏损比例停在同一格。

D1Q3 需要满足的矩只有三个#

与求解 Navier–Stokes 的 LBM 不同，用于标量输运的分布函数 $g_i$ 需要满足的矩要少得多。Guo 在2009年非线性对流扩散方程模型中使用的平衡分布是这样的：

其中 $w_i$ 是格子权重，$\mathbf{c}_i$ 是格子速度，$c_s^2 = 1/3$ 是格子声速的平方。这个分布满足三个矩条件。

和流动平衡分布分开的地方在第三个：这里没有 $\phi\mathbf{u}\mathbf{u}$ 项。速度的二次项从一开始就没放进去。标量方程没有应力张量，二阶矩只要各向同性部分就够，格子也可以用 D2Q5 代替 D2Q9。扩散系数前面 $\tau - 1/2$ 里的那半步，和离散 LBE 的二阶展开留在黏性公式里的是同一项。

看上去够用，而且 $u = 0$ 时确实完美。问题在于，截掉的那部分不是免费的。

Chapman-Enskog 留下的第四项#

给 BGK 格子玻尔兹曼方程加上源项 $S_i$：

按 $g_i = g_i^{(0)} + \varepsilon g_i^{(1)} + \varepsilon^2 g_i^{(2)}$ 展开，并把时间导数拆成 $\partial_t = \varepsilon\partial_{t_1} + \varepsilon^2\partial_{t_2}$。$\varepsilon$ 阶的零阶矩原样给出目标方程的对流部分。

同一个 $\varepsilon$ 阶方程的一阶矩才是核心，因为 $g^{(1)}$ 携带的通量在这里被确定下来。

括号里第二项 $c_s^2\nabla\phi$ 正是我们想要的扩散通量。可第一项也跟着挤了进来。在流动 LBM 中，平衡分布二阶矩里的 $\rho\mathbf{u}\mathbf{u}$ 会把它抵消掉大半，而标量模型没有那一项，所以它留了下来。

把 $\varepsilon^2$ 阶的零阶矩也合并进来，就得到最终形式。

$D = c_s^2(\tau - 1/2)\delta t$ 一如预期。右端第二项则是谁也没有订购的多余通量。在真正的定常态下，$u$ 不随时间变化、$\phi$ 也不再变化时，它会消失。但只要 $\phi$ 还在动，也就是只要计算还在进行，$\partial_t(\phi\mathbf{u})$ 就不是零。

均匀 $u$ 下这一项的真身是负扩散#

取最简单的情形：$u$ 在空间和时间上都是常数。于是 $\partial_t(\phi u) = u\,\partial_t\phi$，而在主导阶上 $\partial_t\phi \simeq -u\,\partial_x\phi$。代入得

这是一个与扩散项符号相反的扩散项。与物理扩散项合并后

亏损的形状就是一个平方。格子声速的平方是三分之一，所以速度 0.3 的平方 0.09 乘以三，得到 0.27。松弛时间为 1 时，标定扩散是三分之一乘二分之一，也就是 0.166667。剩下的是 0.166667 乘 0.73，等于 0.121667。差额 0.045 正好是标定扩散乘 0.27。输出里 $u = 0.30$ 那一行就是这道减法。低速时平方变小。速度 0.05 的平方 0.0025 乘以三是 0.0075，亏损只有 0.75%。峰贴在精确解旁边，眼睛很容易漏掉。

比例不跟松弛时间走，原因在同一行里。松弛时间减去二分之一，这个因子同时乘在物理扩散和多余项上。两者相除，因子就消失了。所以比例与松弛时间无关。把松弛时间从 0.6 抬到 1.5，标定扩散从 0.033333 变成 0.333333，差了十倍。速度 0.25 时，0.0625 的三倍是 0.1875，用 1 去减，剩下 0.8125，四行都钉在这个数上。绝对亏损只是从 0.006 附近一起涨到 0.062 附近。把扩散调大，没有地方把误差埋掉。

速度碰到格子声速，也就是三分之一的平方根的倒数，大约 0.577，括号里就变成 0。再往外，有效扩散是负的。峰不再铺开，而是自己变尖，分布发散。

$1 - 1/(2\tau)$ 又出现了#

现在来定源项 $S_i$ 的系数。要让上面最终方程中的多余项与源项互相抵消，需要

在满足这个条件的同时还保证 $\sum_i S_i = 0$ 的最简单形式只有一个。

$1 - 1/(2\tau)$ 又来了。这个系数和力项在三阶矩上不能越过的那条线在同一个位置。在离散时间的格子上，源项会作用两次：一次通过 $g^{(1)}$，一次通过 Taylor 展开的二阶项。系数 $1/2$ 的差别就产生于此。

如果漏掉这个系数，直接取 $\sum_i \mathbf{c}_i S_i = \partial_t(\phi\mathbf{u})$ 会怎样？在 $\tau = 1$ 时恰好加倍，原本亏损27%的扩散变成过量27%。误差的绝对值没变，所以在对数坐标上测收敛阶也很难察觉。

$\partial_t(\phi\mathbf{u})$ 在代码里只需保存上一时间步的 $\phi u$ 并做后向差分即可。多一个数组就是全部代价。

用60行 Python 量出的 $D_{\text{eff}}$#

不要停在嘴上，量一量。让高斯分布随流走，用最小二乘拟合其二阶矩 $\sigma^2$ 的增长斜率，这个斜率就是 $2D_{\text{eff}}$。

import numpy as np
 
CS2 = 1.0 / 3.0
C = np.array([0, 1, -1])
W = np.array([2 / 3, 1 / 6, 1 / 6])
 
 
def d1q3_equilibrium(phi, u):
    """g_i^eq = w_i phi (1 + c_i u / cs^2) — 只满足三个矩的平衡分布"""
    return np.stack([W[i] * phi * (1.0 + C[i] * u / CS2) for i in range(3)])
 
 
def gaussian_moments(x, phi):
    m0 = phi.sum()
    mean = (x * phi).sum() / m0
    return mean, (((x - mean) ** 2) * phi).sum() / m0
 
 
def run_cde_lbm(L, steps, tau, u, sigma0, x0, corrected):
    x = np.arange(L, dtype=float)
    phi = np.exp(-((x - x0) ** 2) / (2 * sigma0**2))
    g = d1q3_equilibrium(phi, u)
    phi_old = phi.copy()
    hist = []
    for n in range(steps + 1):
        if n % 100 == 0:
            hist.append((n, gaussian_moments(x, phi)[1]))
        src = np.zeros_like(g)
        if corrected and n > 0:
            # S_i = w_i (1 - 1/(2 tau)) c_i d_t(phi u) / cs^2,  dt = 1
            dt_phiu = (1.0 - 1.0 / (2 * tau)) * u * (phi - phi_old)
            for i in range(3):
                src[i] = W[i] * C[i] * dt_phiu / CS2
        geq = d1q3_equilibrium(phi, u)
        g = g - (g - geq) / tau + src          # 碰撞
        for i in range(3):
            g[i] = np.roll(g[i], C[i])         # 迁移
        phi_old = phi
        phi = g.sum(axis=0)
    return np.array(hist)
 
 
def fit_diffusivity(hist):
    """从 sigma^2 = sigma0^2 + 2 D_eff t 的斜率读出 D_eff"""
    return np.polyfit(hist[:, 0], hist[:, 1], 1)[0] / 2.0
 
 
L, STEPS, SIG0, X0 = 800, 1600, 10.0, 80.0
 
tau = 1.0
D = CS2 * (tau - 0.5)
print(f"tau = {tau},  D = cs^2 (tau-1/2) = {D:.6f},  cs^2 = {CS2:.6f}")
print(f"{'u':>6} {'u^2/cs^2':>9} | {'D_eff (no src)':>14} {'ratio':>7} {'1-u^2/cs^2':>11} |"
      f" {'D_eff (src)':>12} {'ratio':>7}")
for u in [0.05, 0.10, 0.20, 0.30]:
    d_raw = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, False))
    d_fix = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, True))
    print(f"{u:>6.2f} {u * u / CS2:>9.4f} | {d_raw:>14.6f} {d_raw / D:>7.4f} {1 - u * u / CS2:>11.4f} |"
          f" {d_fix:>12.6f} {d_fix / D:>7.4f}")
 
u = 0.25
print(f"\nu = {u} fixed, tau sweep   (theory: ratio = 1 - u^2/cs^2 = {1 - u * u / CS2:.4f}, tau-independent)")
print(f"{'tau':>6} {'D':>10} | {'D_eff (no src)':>14} {'ratio':>7} | {'D_eff (src)':>12} {'ratio':>7}")
for tau in [0.6, 0.8, 1.0, 1.5]:
    D = CS2 * (tau - 0.5)
    d_raw = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, False))
    d_fix = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, True))
    print(f"{tau:>6.1f} {D:>10.6f} | {d_raw:>14.6f} {d_raw / D:>7.4f} | {d_fix:>12.6f} {d_fix / D:>7.4f}")

输出如下。

tau = 1.0,  D = cs^2 (tau-1/2) = 0.166667,  cs^2 = 0.333333
     u  u^2/cs^2 | D_eff (no src)   ratio  1-u^2/cs^2 |  D_eff (src)   ratio
  0.05    0.0075 |       0.165417  0.9925      0.9925 |     0.166666  1.0000
  0.10    0.0300 |       0.161667  0.9700      0.9700 |     0.166666  1.0000
  0.20    0.1200 |       0.146667  0.8800      0.8800 |     0.166663  1.0000
  0.30    0.2700 |       0.121667  0.7300      0.7300 |     0.166658  0.9999
 
u = 0.25 fixed, tau sweep   (theory: ratio = 1 - u^2/cs^2 = 0.8125, tau-independent)
   tau          D | D_eff (no src)   ratio |  D_eff (src)   ratio
   0.6   0.033333 |       0.027096  0.8129 |     0.033345  1.0004
   0.8   0.100000 |       0.081258  0.8126 |     0.100006  1.0001
   1.0   0.166667 |       0.135417  0.8125 |     0.166661  1.0000
   1.5   0.333333 |       0.270794  0.8124 |     0.333275  0.9998

把 u 调大，只有红色的多余通量按 $u^2$ 长大。把 alpha 从 0 推到 1，源项就把这份宽度抵掉。

alpha = 0 时，红色波瓣把蓝色物理通量上下翻过来。alpha 到 1，绿色盖住红色，黄色的合计回到蓝线。源项打开那一列的 ratio 贴上 1.0000，就是这一幕。

第一张表的 ratio 列与 1-u^2/cs^2 列小数点后四位完全一致。这说明预测不是近似，而是精确的主导阶。打开源项后，四个速度全部回到1.0000。

第二张表更刺人。$\tau$ 从0.6扫到1.5，$D$ 变了十倍，亏损比例却只从0.8129挪到0.8124。想靠把扩散取大来掩盖误差是行不通的：$D$ 放大十倍，消失的量也放大十倍。

格子速度的预算早已被占用

$u^2/c_s^2$ 就是格子马赫数的平方。$u = 0.1$ 时 $0.01 \times 3 = 0.03$，扩散少 3%。$u = 0.2$ 时 $0.04 \times 3 = 0.12$，少 12%。0.05 及以下只有 0.75%，埋在网格误差里。因为压缩性误差把 $u$ 拴在 0.1 附近的习惯，还没看见标量那一侧的账单。界面厚度或 Sherwood 数要当成数字报出去，$u$ 一过 0.1 就该打开源项。再低一点，省一个数组，表里的 0.9925 还站得住。

phase-field 不能只信 $u^2/c_s^2$。多余项是完整的 $\partial_t(\phi\mathbf{u})$，平方形式只在 $u$ 均匀且定常时成立。界面附近 $u$ 会跳，$\partial_t u$ 也不为零，$\phi\,\partial_t u$ 那一半就活下来。$\partial_t(\phi\mathbf{u})$ 用上一时刻的 $\phi u$ 做后向差分放进去。$1 - 1/(2\tau)$ 里的 $\tau$ 是一阶矩的松弛时间。速度参考系一移，黏性怎么掉，和平衡分布缺掉的三阶矩是同一类账单。迁移率的算术怎么复核都对，界面却在变薄，那就不是 $\tau$ 错了。流动先把扩散的预算拿走了。