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

黏性交给β，稳定性交给那个让离散H在碰撞前后保持不变的非平凡根α。

把熵只当成“丢失的信息”，解释不了这组计算为什么不炸。格子玻尔兹曼方法(Lattice Boltzmann Method, LBM)里，为了压低黏性而把 $\tau$ 贴向 $1/2$ ，BGK就会走过平衡。分布掉成负数。Karlin这一支的熵LBM(Entropic LBM, ELBM)在原地把离散H绑在碰撞前的值上。代价是每个节点一个非线性根。

把非平衡扰动滑块向右推，黄点会离开虚线 $\alpha=2$ ， $G(\alpha)=0$ 的第二个根向左移动。

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.

扰动小时，点贴在2上。虚线和黄点之间的空隙，就是BGK不肯丢掉的过松弛宽度。

当τ趋近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的二阶截断是Hermite展开在二阶截断的多项式平衡，它不是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，把过松弛削掉。

把黏性公式摊成数字，过松弛离墙壁有多近就清楚了。 $c_s^2 = 1/3$ ， $\nu = c_s^2(\tau - 1/2)$ 。 $\tau = 0.51$ 时差是 $0.01$ ，所以 $\nu = 1/300 \approx 0.00333$ 。松弛率 $\omega = 1/0.51 \approx 1.961$ 。偏离平衡的长度1被走了1.961倍，于是落在平衡另一侧0.961处。再贴到 $\tau = 0.501$ ， $\omega \approx 1.996$ 。离镜像 $\omega = 2$ 只剩0.004。

代码里的 $\beta = 0.95$ 是同一个关系。 $\tfrac{1}{2\beta}-\tfrac{1}{2} = \tfrac{1}{1.9}-\tfrac{1}{2} = \tfrac{1}{38}$ ，因此 $\nu = 1/114 \approx 0.00877$ 。同样黏性的BGK对应 $\tau = 1/(2\beta) \approx 0.526$ 。平衡附近若 $\alpha \approx 2$ ，有效松弛率 $\alpha\beta = 1.9$ ，黏性与这次BGK相同。非平衡变大、 $\alpha$ 降到2以下时， $\alpha\beta$ 也降到1.9以下。越过平衡的长度 $\alpha\beta-1$ 比0.9更短。黏性仍是 $\beta$ 定下的那个值。被削掉的只有过松弛。

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的修正量。

两行扰动可以读成矩。 $f_1$ 加 $0.03$ 、 $f_3$ 减 $0.03$ ，质量和不变。方向1和3的格子速度是 $+1$ 和 $-1$ ，所以 $x$ 动量增加 $0.03-(-0.03)=0.06$ 。起初 $\rho=1$ 、 $u_x=0.05$ ，扰动刚加上时 $u_x=0.11$ 。乘积因子是 $(2u+\sqrt{1+3u^2})/(1-u)$ 。 $u=0.05$ 时约1.162， $u=0.11$ 时约1.391。无量纲速度0.3时约2.467，分母 $1-u$ 只剩0.7。 $|u| \lesssim 0.1$ 选的是这个因子还比较缓的一段。

二分法的上界也来自负号。 $d_i$ 为负的分量碰到0的位置是 $\alpha=-f_i/d_i$ ，其中最小的就是 $a_{\max}$ 。搜索停在 $\min(0.999\,a_{\max},\,4)$ 。超过4、或者会改分布符号的根，不进候选。

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

前面那条曲线就是这个 $G$ 。黄点是非平凡根 $\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。红灯在那时亮。H的这次上升，就是低黏性格子裂开的起点。

实现中遇到的三件事

每个节点都做非线性求解，成本会比黏性项更先成为负担。平衡附近， $\|f - f^{eq}\|$ 很小的时候取 $\alpha \approx 2$ ，只在非平衡越过阈值的格点上打开二分法。成本大约降到一成。

乘积形式的平衡只在 $|u| < 1$ 有定义。分母里有 $1 - u$ 。无量纲速度超过0.3，压力张量开始按方向扭曲。即便是ELBM，放在 $|u| \lesssim 0.1$ 里也更稳。

平衡若仍用二阶截断多项式，H取最小这件事就没有保证了。稳定子的好处要拿全，平衡本身也得是熵乘积。多项式上只加一个 $\alpha$ ，只剩一半。

只抬β、把α留在2，H的上界就松开了#

平衡附近直接取 $\alpha \approx 2$ ，只是 $\|f - f^{eq}\|$ 很小时的省成本办法。把 $\beta$ 往1推，又把这个捷径固定在全部节点上，有效松弛率 $\alpha\beta$ 会重新贴上2。直到镜像的过松弛回来了，离散H可以比碰撞前更大。黏性降了，根却没有解。

碰撞可以把H压平，可一旦迁移把邻点上已经为负的分布带过来，这个根就晚了。符号若在对流阶段被打破，保正的通量限制器是否仍要单独拦住，还是熵条件自己就够，是下一处分叉。