熵阻止了爆破 — Karlin的熵格子玻尔兹曼方法

熵通常被当作"丢失的信息"来介绍——无序增长的方向，不可逆损失的度量。然而在格子玻尔兹曼方法(Lattice Boltzmann Method, LBM)中，正是这个量抓住了失控的模拟并把它稳住。本文讨论Karlin系的熵LBM(Entropic LBM, ELBM)如何在离散化层面强制H定理，从而消除低黏性区的发散。读完后你会明白，为什么值得为碰撞算子中的一行而在每一步求解一个非线性方程。

当τ趋近0.5时BGK崩溃#

标准LBM使用BGK碰撞(单松弛时间近似)。它把分布 $f_i$ 按因子 $\omega$ 拉向平衡 $f_i^{eq}$ 。

f_i^* = f_i + \omega\,(f_i^{eq} - f_i)

这里 $f_i$ 是格子方向 $i$ 的分布， $f_i^{eq}$ 是平衡分布， $\omega = 1/\tau$ 是松弛率。

黏性由 $\nu = c_s^2(\tau - 1/2)$ 确定，其中 $c_s^2 = 1/3$ 是格子声速的平方。要解高Reynolds(低黏性)流动，就要把 $\tau \to 1/2$ ，即 $\omega \to 2$ 。

麻烦在于， $\omega$ 越接近2，碰撞越会越过平衡而落到另一侧。分布会跌为负值，而一旦某个 $f_i$ 变负，它在下一步会振荡得更厉害。格子炸掉。BGK没有任何东西能阻止这种过松弛(over-relaxation)。

玻尔兹曼的H定理 — 时间之箭#

1872年，玻尔兹曼证明了以下量在仅有碰撞时永不增长。

H = \int f \ln f \, d\mathbf{c}

$f$ 是速度分布， $\mathbf{c}$ 是分子速度。有 $dH/dt \le 0$ ，等号仅在平衡(Maxwell-Boltzmann分布)处成立。

H是负熵。H单调下降与熵单调上升是同一陈述。这就是不可逆性，即时间之箭。关键在于这个不等式是一份稳定性证书。一个H永不增长的数值格式不可能发散，因为振幅无界增长需要H增加。

离散H函数与熵平衡#

LBM把连续速度切成9个方向(D2Q9)。H也相应离散化。

H(f) = \sum_i f_i \ln\frac{f_i}{w_i}

$w_i$ 是格子权重(D2Q9中：中心 $4/9$ ，面 $1/9$ ，对角 $1/36$ )。它是按权重归一后比值的对数。

熵平衡定义为在守恒矩约束下最小化这个 $H$ 的分布。固定质量 $\sum_i f_i = \rho$ 和动量 $\sum_i \mathbf{c}_i f_i = \rho\mathbf{u}$ ，对H做最小化(拉格朗日乘子问题)。在D2Q9中存在闭式的乘积(product)形式。

f_i^{eq} = \rho\, w_i \prod_{\alpha=x,y}\left(2 - \sqrt{1+3u_\alpha^2}\right)\left(\frac{2u_\alpha + \sqrt{1+3u_\alpha^2}}{1 - u_\alpha}\right)^{c_{i\alpha}}

$u_\alpha$ 是速度分量， $c_{i\alpha}$ 是格子速度分量(−1、0、1)。与标准LBM所用的二阶截断多项式平衡不同，这个平衡是H的真正最小点。

熵稳定子 — 每步求解α#

ELBM的碰撞与BGK形状相同，但把松弛率一分为二。

f_i^* = f_i + \alpha\beta\,(f_i^{eq} - f_i)

$\beta \in (0,1)$ 决定黏性( $\nu = c_s^2(\tfrac{1}{2\beta} - \tfrac{1}{2})$ ， $\beta\to1$ 为无黏)， $\alpha$ 决定稳定性。

$\alpha$ 不是常数。在每一步、每个格点上，它作为下式的非平凡根被求解。

H\big(f + \alpha(f^{eq} - f)\big) = H(f)

含义很清楚：只朝平衡松弛到使碰撞后的H与碰撞前完全相同的程度。H的任何增加都被结构性地禁止。 $\alpha=0$ 是平凡根，我们要的是平衡另一侧的第二个根。在平衡附近 $\alpha \approx 2$ ，恰好就是BGK。只有当非平衡部分很大时， $\alpha$ 才偏离2并削掉过松弛。

Python — 用二分法求解熵α#

下面是乘积形平衡、H函数，以及用二分法求解 $\alpha$ 的一步。

import numpy as np
 
# D2Q9 格子权重与速度 (lattice weights & velocities)
W  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
CX = np.array([0, 1, 0, -1,  0, 1, -1, -1,  1])
CY = np.array([0, 0, 1,  0, -1, 1,  1, -1, -1])
 
def h_function(f):
    # 离散H函数: H = Σ f_i ln(f_i / w_i)  (玻尔兹曼H的格子版)
    return np.sum(f * np.log(f / W))
 
def entropic_equilibrium(rho, ux, uy):
    # 精确的Ansumali–Karlin乘积(product)形熵平衡
    def factor(u):
        s = np.sqrt(1 + 3 * u * u)
        return (2 - s), (2 * u + s) / (1 - u)
    bx, rx = factor(ux)
    by, ry = factor(uy)
    return rho * W * bx * by * (rx ** CX) * (ry ** CY)
 
def entropic_alpha(f, feq):
    # 用二分法求 H(f + α(feq − f)) = H(f) 的非平凡根 α*
    d = feq - f
    H0 = h_function(f)
    G = lambda a: h_function(f + a * d) - H0
    amax = np.min(np.where(d < 0, -f / d, np.inf))   # 保证正值的上界
    hi = min(0.999 * amax, 4.0) if np.isfinite(amax) else 4.0
    if G(hi) <= 0:
        return 2.0                                   # 平衡附近的安全值
    lo = 1.0
    while G(lo) > 0 and lo > 1e-3:
        lo *= 0.5
    for _ in range(60):
        mid = 0.5 * (lo + hi)
        if G(mid) > 0: hi = mid
        else:          lo = mid
    return 0.5 * (lo + hi)
 
# --- 单节点，单步 ---
feq = entropic_equilibrium(1.0, 0.05, 0.0)
f = feq.copy()
f[1] += 0.03; f[3] -= 0.03            # 由streaming带入的非平衡 (non-equilibrium)
rho = f.sum(); ux = (CX * f).sum() / rho; uy = (CY * f).sum() / rho
feq = entropic_equilibrium(rho, ux, uy)
 
alpha = entropic_alpha(f, feq)
beta  = 0.95                          # ν = c_s²(1/(2β) − 1/2), β→1 为低黏性
f_post = f + alpha * beta * (feq - f)
 
print(f"熵 α     = {alpha:.4f}   (BGK固定值 2.0)")
print(f"H 碰撞前 = {h_function(f):.8f}")
print(f"H 碰撞后 = {h_function(f_post):.8f}")
print(f"最小分布 = {f_post.min():.3e}  (> 0)")

输出显示 $\alpha$ 接近2但并非恰好为2。那个微小的差距就是守恒H的修正量。

交互 — G(α)的第二个根#

在下面的模拟中亲手操作。增大非平衡扰动，观察曲线如何弯曲。

Non-equilibrium perturbation: 0.060

The yellow dot is the non-trivial root α* where the discrete H-function returns to its pre-collision value. Near equilibrium α* ≈ 2, which is exactly standard BGK. Push the perturbation up and α* drifts away from 2 — that gap is the entropic correction BGK ignores.

黄点就是我们要的 $\alpha^*$ 。扰动小时它紧贴在虚线( $\alpha=2$ )上。把滑块向右推， $\alpha^*$ 就从2脱开。BGK忽略这一移动，始终用2。当差距累积起来，格子就发散。

交互 — 阻止H增长的那一刻#

把同一个非平衡状态用两种碰撞规则并排处理。红色柱是变为负值的分布。

Non-equilibrium shear: 0.070

Both panels show post-collision populations. The entropic collision lands exactly on the equal-entropy point, so H_ent = H_pre. Fixed-α BGK misses it; once the non-equilibrium part grows, BGK over-shoots and H actually increases — a violated discrete H-theorem, which is the seed of low-viscosity blow-up.

熵这一侧始终满足 $H_{ent} = H_{pre}$ ——它就是这么设计的。BGK这一侧，一旦非平衡部分增大，就会出现 $H_{BGK}$ 超过碰撞前的时刻。那时红色警告就亮起。正是这个H的增长，成为低黏性下撕裂格子的种子。

实现中遇到的三件事#

第一，每个节点的非线性求解很昂贵。所以实际中在平衡附近( $\|f - f^{eq}\|$ 很小时)取 $\alpha \approx 2$ ，只在非平衡越过阈值的格点上启动二分法。成本降到约十分之一。

第二，乘积形平衡只在 $|u| < 1$ 时有定义——分母里有 $1 - u$ 。当无量纲速度超过0.3，压力张量开始各向异性地扭曲，所以即便是ELBM，在 $|u| \lesssim 0.1$ 的区域使用也最稳妥。

第三，若直接沿用二阶截断多项式平衡，H最小化的保证就会失效。要完整获得熵稳定子的好处，平衡本身也必须是熵乘积形。在多项式平衡上只加 $\alpha$ 只有一半效果。

三行小结#

BGK在低黏性下崩溃，是因为过松弛增大了离散H。ELBM每步求解使碰撞后H等于碰撞前H的 $\alpha$ 来阻止它。
$\alpha$ 是 $H(f + \alpha(f^{eq}-f)) = H(f)$ 的非平凡根，在平衡附近 $\alpha \approx 2$ 自动恢复BGK。
黏性由 $\beta$ 承担，稳定性由 $\alpha$ 承担。代价是每步一次非线性求解，但平衡附近的捷径让它在实际中可以承受。