在BGK发散的低粘性区域中存活 — MRT格子玻尔兹曼碰撞

只把s_ν推到1.99、粘性接近0时，留下来的s_e和s_q柱条决定棋盘格出不出现。

同一张格子、同一个平衡，一边平滑地转，另一边碎成棋盘格。BGK用一个弛豫时间推九个分布函数。MRT(Multiple Relaxation Time，多重弛豫时间)先变成矩，再以各自的速率推。

把 $s_\nu$ 柱条推到1.99，粘性读数往0掉。 $s_e$ 和 $s_q$ 柱条留在原地。

D2Q9 moment-space relaxation rates · S = diag(s₀…s₈)

sₑ (bulk) = 1.64

s_ε = 1.54

s_q (ghost) = 1.70

s_ν (shear) = 1.20

shear viscosity ν = cs²(1/s_ν − 1/2) = 0.1111 (lattice units)

✓ all rates in (0, 2) — linearly stable region

守恒量柱条 $\rho$ 、 $j_x$ 、 $j_y$ 钉在0上。BGK把这些柱条绑成同一个数。低粘性下棋盘格出现的地方，就是这根绑绳。

一拍弛豫所有量的BGK弱点#

格子玻尔兹曼方法(LBM)重复两个步骤:分布函数 $f_i$ 在格子上迁移(streaming),在每个点向平衡弛豫(collision)。最常见的BGK碰撞写作

f_i(\mathbf{x}+\mathbf{c}_i,\, t+1) = f_i(\mathbf{x}, t) - \frac{1}{\tau}\bigl(f_i - f_i^{\text{eq}}\bigr)

其中 $f_i$ 是沿速度 $\mathbf{c}_i$ 的分布函数, $\tau$ 是弛豫时间, $f_i^{\text{eq}}$ 是局部平衡分布。粘性仅由 $\tau$ 决定。

\nu = c_s^2\left(\tau - \tfrac{1}{2}\right), \qquad c_s^2 = \tfrac{1}{3}

问题从这里开始。要降低粘性(提高雷诺数),必须把 $\tau$ 逼近 $1/2$ 。可是当 $\tau \to 1/2$ 时,与粘性毫无关系的高阶模几乎不弛豫,在格子上四处游走。这些模在格子分辨率极限处放大幅度时,就出现棋盘振荡,计算发散。在BGK中,调粘性的旋钮和调稳定性的旋钮被绑在同一个手柄上。

从矩空间审视碰撞

不要直接弛豫九个分布函数。先把它们变成有物理意义的矩,再弛豫。矩是分布函数的线性组合。

\mathbf{m} = M\,\mathbf{f}

这里 $\mathbf{f} = (f_0,\dots,f_8)^\top$ 是分布函数向量, $M$ 是 $9\times 9$ 变换矩阵, $\mathbf{m}$ 是矩向量。碰撞在矩空间向平衡矩 $\mathbf{m}^{\text{eq}}$ 弛豫,再变换回速度空间。

\mathbf{f}^{*} = \mathbf{f} - M^{-1}\,S\,\bigl(M\mathbf{f} - \mathbf{m}^{\text{eq}}\bigr)

$S = \mathrm{diag}(s_0,\dots,s_8)$ 是装着各矩弛豫率的对角矩阵。若把所有 $s_i$ 设为同一值 $1/\tau$ ,此式就精确退回BGK。MRT所做的不过是把那些对角元逐矩单独求解。

D2Q9的九个矩与变换矩阵M#

在D2Q9格子(二维、九速度)上,矩有以下九个:密度 $\rho$ 、能量 $e$ 、能量平方 $\varepsilon$ 、动量 $j_x, j_y$ 、能量通量 $q_x, q_y$ ,以及对应应力的 $p_{xx}, p_{xy}$ 。

其中 $\rho, j_x, j_y$ 是碰撞不改变的守恒量,因此其弛豫率固定为 $s=0$ 。其余六个是可调旋钮。Lallemand–Luo标准变换矩阵被设计成各行正交,因此逆矩阵就是按行范数缩放的转置。

平衡矩仅由密度和动量闭合。

\begin{aligned} e^{\text{eq}} &= -2\rho + 3(j_x^2 + j_y^2), &\quad \varepsilon^{\text{eq}} &= \rho - 3(j_x^2 + j_y^2) \\ q_x^{\text{eq}} &= -j_x, &\quad q_y^{\text{eq}} &= -j_y \\ p_{xx}^{\text{eq}} &= j_x^2 - j_y^2, &\quad p_{xy}^{\text{eq}} &= j_x j_y \end{aligned}

这里 $j_x = \rho u_x$ 、 $j_y = \rho u_y$ 是动量分量。应力矩 $p_{xx}, p_{xy}$ 的平衡与速度梯度耦合,因此它们的弛豫率决定粘性。

每个矩不同的弛豫率 — 与粘性无关的自由度

决定剪切粘性的是应力矩的弛豫率 $s_\nu$ 。

\nu = c_s^2\left(\frac{1}{s_\nu} - \frac{1}{2}\right)

其余弛豫率 $s_e$ (体积)、 $s_\varepsilon$ 、 $s_q$ (能量通量)不进粘性。它们是自由参数。把 $s_\nu$ 贴到 $2$ ,高阶模仍可单独放在稳定区间 $1.5\sim1.8$ 。BGK绑在一个旋钮上的东西,在这里分开。

代码默认粘性是 $\nu=0.004$ 。 $s_\nu = 1/(3\nu + 1/2)$ ,所以 $3\times 0.004 + 0.5 = 0.512$ , $s_\nu = 1/0.512 = 1.953125$ 。前面的柱条推到1.99时, $\nu = (1/3)(1/1.99 - 1/2)$ 。 $1/1.99 \approx 0.5025$ 减去0.5再除以3, $\nu \approx 8.4\times 10^{-4}$ 。靠近0,不是0。这个粘性上BGK把 $s_e$ 、 $s_\varepsilon$ 、 $s_q$ 也一起推到1.99。代码里的MRT把这三个位置固定在 $1.64$ 、 $1.54$ 、 $1.7$ 。跑起来的数组是 [0, 1.64, 1.54, 0, 1.7, 0, 1.7, s_nu, s_nu],等于0的只有三个守恒位。

平衡矩只由密度和动量闭合,是因为分布已经被低阶截断。那一刀在九个速度里的 Maxwell–Boltzmann 那边。

碰撞在矩空间,迁移在速度空间

算法流程如下。

由 $f_i$ 计算宏观量 $\rho, u_x, u_y$
用 $\mathbf{m} = M\mathbf{f}$ 变换到矩
计算平衡矩 $\mathbf{m}^{\text{eq}}$ (由密度与动量)
逐矩弛豫: $m_a \leftarrow m_a - s_a(m_a - m_a^{\text{eq}})$
用 $\mathbf{f}^{*} = M^{-1}\mathbf{m}$ 复原到速度空间
迁移:把 $f_i^{*}$ 移到沿 $\mathbf{c}_i$ 的邻居
施加边界条件后回到第1步

只有碰撞在矩空间发生,迁移照常在速度空间进行。额外开销仅为每个单元两次 $9\times 9$ 矩阵乘法。

Python — 用Taylor–Green涡测量粘性衰减#

验证使用Taylor–Green涡。该流动的粘性衰减精确解已知,因此可以用数值核对实现是否正确。动能应按 $\exp(-4\nu k^2 t)$ 衰减( $k$ 为涡的波数)。

import numpy as np
 
# D2Q9 格子
cx = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
cy = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
 
# Lallemand-Luo 变换矩阵 (行: rho,e,eps,jx,qx,jy,qy,pxx,pxy)
M = np.array([
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [-4,-1,-1,-1,-1, 2, 2, 2, 2],
    [ 4,-2,-2,-2,-2, 1, 1, 1, 1],
    [ 0, 1, 0,-1, 0, 1,-1,-1, 1],
    [ 0,-2, 0, 2, 0, 1,-1,-1, 1],
    [ 0, 0, 1, 0,-1, 1, 1,-1,-1],
    [ 0, 0,-2, 0, 2, 1, 1,-1,-1],
    [ 0, 1,-1, 1,-1, 0, 0, 0, 0],
    [ 0, 0, 0, 0, 0, 1,-1, 1,-1],
], dtype=float)
Minv = np.linalg.inv(M)
 
def d2q9_equilibrium(rho, ux, uy):
    eq = np.zeros((9,) + rho.shape)
    for i in range(9):
        cu = cx[i] * ux + cy[i] * uy
        eq[i] = w[i] * rho * (1 + 3*cu + 4.5*cu**2 - 1.5*(ux**2 + uy**2))
    return eq
 
def relaxation_rates(nu):
    s_nu = 1.0 / (3*nu + 0.5)          # 决定粘性的弛豫率
    # 守恒(rho,jx,jy)=0,高阶模为稳定区间的固定值
    return np.array([0, 1.64, 1.54, 0, 1.7, 0, 1.7, s_nu, s_nu])
 
def taylor_green_init(N, U0):
    x = np.arange(N)
    X, Y = np.meshgrid(x, x, indexing='ij')
    k = 2*np.pi / N
    ux = -U0 * np.cos(k*X) * np.sin(k*Y)
    uy =  U0 * np.sin(k*X) * np.cos(k*Y)
    rho = 1 - 0.75 * U0**2 * (np.cos(2*k*X) + np.cos(2*k*Y))
    return d2q9_equilibrium(rho, ux, uy), k
 
def macros(f):
    rho = f.sum(axis=0)
    ux = (cx[:, None, None] * f).sum(axis=0) / rho
    uy = (cy[:, None, None] * f).sum(axis=0) / rho
    return rho, ux, uy
 
def mrt_collide(f, s):
    rho, ux, uy = macros(f)
    jx, jy = rho*ux, rho*uy
    m = np.tensordot(M, f, axes=([1], [0]))        # 进入矩空间
    meq = np.zeros_like(m)
    meq[0] = rho
    meq[1] = -2*rho + 3*(jx**2 + jy**2)
    meq[2] = rho - 3*(jx**2 + jy**2)
    meq[3], meq[4] = jx, -jx
    meq[5], meq[6] = jy, -jy
    meq[7] = jx**2 - jy**2
    meq[8] = jx*jy
    m -= s[:, None, None] * (m - meq)               # 逐矩弛豫
    return np.tensordot(Minv, m, axes=([1], [0]))    # 复原到速度空间
 
def stream(f):
    for i in range(9):
        f[i] = np.roll(f[i], (cx[i], cy[i]), axis=(0, 1))
    return f
 
def run_mrt_lbm(N=64, nu=0.004, U0=0.04, steps=2000):
    f, k = taylor_green_init(N, U0)
    s = relaxation_rates(nu)
    e0 = None
    for t in range(steps):
        f = mrt_collide(f, s)
        f = stream(f)
        if t % 400 == 0:
            _, ux, uy = macros(f)
            E = np.mean(ux**2 + uy**2)
            e0 = E if e0 is None else e0
            print(f"t={t:5d}  KE/KE0={E/e0:.4f}  analytic={np.exp(-4*nu*k**2*t):.4f}")
 
run_mrt_lbm()

print 并排放的解析值是 $\exp(-4\nu k^2 t)$ 。默认 $N=64$ 、 $\nu=0.004$ 、 $k=2\pi/N$ ,于是 $4\nu k^2 = 0.016\,(\pi/32)^2 \approx 1.54\times 10^{-4}$ 。 $t=400$ 时指数约 $-0.0617$ ,解析因子约0.940。 $t=1600$ 时约 $-0.247$ ,因子约0.781。range(2000) 使得最后印出的 $t$ 是1600。说到小数点后两位相合,比的就是这个因子和数值动能。同一格子、同一 $\nu$ 上跑 BGK,会在更低的粘性下先发散。

在相同Re下并排比较BGK与MRT#

粘性滑块把衰减放慢,BGK/MRT 按钮决定同一个 $\nu$ 下格子还站不站得住。

viscosity ν = 0.0060 (Re ≈ 267)

mode: MRT

step: 0

KE / KE₀: 1.000

analytic e⁻⁴ᵛᵏ²ᵗ: 1.000

Color = vorticity (red ↺ / blue ↻). Taylor–Green vortex on a 40×40 periodic lattice.

把粘性降到 $0.002$ 附近并停在 BGK,格子上会爬出粗糙花纹并显示 "diverged"。同一粘性下切到 MRT 再 reset,涡就平滑衰减。 $\nu=0.002$ 时 $s_\nu=1/(3\times 0.002+0.5)=1/0.506\approx 1.976$ 。BGK 把高阶模也推到这个数。代码里的 $1.64$ 、 $1.54$ 、 $1.7$ 留在旁边。

只把 s_ν 贴向 2，高阶模就会留下#

MRT 做的事就是把粘性旋钮 $s_\nu$ 和稳定性旋钮 $s_e$ 、 $s_q$ 在物理上分开。把 $\nu=0.004$ 的 $s_\nu=1.953$ 抄到每个模上,那是 BGK。留下 $1.64$ 、 $1.54$ 、 $1.7$ ,只抬 $s_\nu$ ,那是 MRT。额外开销是每格两次 $9\times 9$ 矩阵乘。

没有这道分离,低粘性的 Taylor–Green 会死成棋盘格。用熵限制拦住同一种失控,是另一条路：把 H 定理放进碰撞的格子玻尔兹曼。衰减对不对,仍是 $\exp(-4\nu k^2 t)$ 这一行。