九支箭里的Maxwell–Boltzmann — LBM平衡分布为何是多项式

1986年,Frisch、Hasslacher、Pomeau把粒子放在格子上,让它们模仿Navier–Stokes方程。这就是LGCA(Lattice Gas Cellular Automaton)。它没撑过五年 — 每个格子方向上"0或1个粒子"的整数条件产生了无法抹去的噪声。1992年McNamara和Zanetti把整数换成实数,用粒子分布函数 $f_i$ 填满格子。噪声消失了 — 但一个新问题接踵而至:9个实数怎么可能替代无穷维的Maxwell–Boltzmann分布?本文追寻这一答案 — Hermite多项式的二阶截断造就了LBM的平衡分布,而同样的截断把LBM拴在低Mach数流动上。

从整数到实数 — 1986年LGCA之死与LBM的诞生#

LGCA的想法简洁优雅。正方格子、6个粒子方向、一条碰撞规则。宏观平均之下,Navier–Stokes重新浮现。但布尔粒子(0或1)放大了统计噪声。要靠时间平均压制噪声,就得让平均窗口足够长,而你本想看的transient流动也就一同消失了。

McNamara和Zanetti把整数粒子换成了实值分布 $f_i(\mathbf{x}, t)$ ,放在每个格子节点上。 $f_i$ 是在时刻 $t$ 、位置 $\mathbf{x}$ 沿格子方向 $\mathbf{c}_i$ 运动的粒子的 期望数量。噪声从结构上被消除。每一步分两段 — collision(碰撞)与streaming(传播)。

f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f_i(\mathbf{x}, t) - \frac{f_i(\mathbf{x}, t) - f_i^{eq}(\rho, \mathbf{u})}{\tau}

$\tau$ 是relaxation time(松弛时间), $f_i^{eq}$ 即 平衡分布函数。碰撞项以 $\tau$ 的时间尺度削去非平衡部分,streaming把削好的分布移到相邻格子。问题归结到 — $f_i^{eq}$ 究竟是什么,使9个实数就够用?

Maxwell–Boltzmann怀里的定时炸弹#

连续动理论中,平衡分布是Maxwell–Boltzmann分布。

f^{MB}(\mathbf{c}) = \rho \, (2\pi R T)^{-D/2} \exp\!\left( -\frac{(\mathbf{c} - \mathbf{u})^2}{2 R T} \right)

$\mathbf{c}$ 是粒子的分子速度, $\mathbf{u}$ 是宏观流速, $D$ 是空间维数, $T$ 为温度, $R$ 为气体常数。这个公式的炸弹就是 $\mathbf{c}$ 是 连续变量。在格子上只留9或19个离散方向,就必须巧妙地"瘦身"这一连续分布。

直觉尝试 — 把Gaussian积分用中点法离散在固定网格上。失败。Navier–Stokes所需的几个矩 — $\rho$ 、 $\rho \mathbf{u}$ 、动量张量 $\rho \mathbf{u}\mathbf{u}$ 、能量 — 与格子求和对不上。守恒一破,LBM当场作废。

正解是函数展开。把Maxwell–Boltzmann展成 正交多项式级数,保留所需的矩,丢掉其余。这一族正交多项式就是Hermite。

Hermite多项式把分布切成碎片#

概率论意义下的Hermite多项式 $He_n(c)$ 在标准Gaussian权重 $w(c) = (2\pi)^{-1/2} e^{-c^2/2}$ 下正交。

\int He_n(c)\, He_m(c)\, w(c)\, dc = n!\, \delta_{nm}

前几个是 $He_0 = 1$ 、 $He_1 = c$ 、 $He_2 = c^2 - 1$ 、 $He_3 = c^3 - 3c$ 。把1D Maxwell–Boltzmann取为等温( $RT = 1$ 无量纲化),指数生成函数立刻给出下式。

f^{MB}(c) = w(c)\, \rho \exp\!\left( u c - \frac{u^2}{2} \right) = w(c)\, \rho \sum_{n=0}^{\infty} \frac{u^n}{n!} He_n(c)

每一项抓住物理世界的一片。 $n=0$ 是静止Gaussian; $n=1$ 是微小drift; $n=2$ 是平均动能修正。更高的 $n$ 是越来越不对称的尾部效应。LBM的整个目标就是把这个级数在 有限阶数 上截断。

Maxwell–Boltzmann과 그 Hermite 절단의 비교. 표의 세 모멘트(밀도·운동량·에너지)가 같아지는 차수가 바로 LBM이 멈추는 자리다.

drift velocity u = 0.40 (= |u|/c_s)truncation order N = 2

N=2면 ρ·ρu·에너지 세 모멘트가 정확히 일치한다. u가 |u|/c_s ≈ 1을 넘어가면 절단오차가 폭주하며 격자가 깨진다 — LBM이 저속(낮은 마하수)에 묶이는 이유.

在下面的模拟里亲自试一下。把drift $u$ 推到 0.4 并选择截断阶 $N=2$ — 前三个矩(密度、动量、能量)与Maxwell–Boltzmann精确吻合。再把 $u$ 推过 1.0,虚线立刻偏离实线。这正是LBM的Mach数约束 $|\mathbf{u}|/c_s \lesssim 0.3$ 的来源,一张图就能看清。

二阶截断的含义 — 被保住的矩#

要复刻一阶粘性Navier–Stokes,需要哪些矩?答案恰好是 — $\rho$ 、 $\rho \mathbf{u}$ ,以及动量通量 $\Pi_{\alpha\beta} = \rho u_\alpha u_\beta + p \delta_{\alpha\beta}$ ,共 三阶矩。因为 $\rho \mathbf{u}\mathbf{u}$ 是 $u$ 的二阶,所以分布函数里必须有 $u^2$ 项。

因此把Hermite级数在 $N=2$ 处截断。

f^{eq}(c) = w(c)\, \rho \left[ 1 + uc + \frac{u^2}{2}(c^2 - 1) \right]

这是LBM平衡分布的 连续形式。 $N=3$ 及以上对非粘性的高阶矩(热通量)有意义,但对标准isothermal LBM是冗余。被丢掉的项造就了LBM两个广为人知的限制 — 等温假设与Mach数0.3的上限。

D2Q9是如何造出来的#

要把连续积分 $\int f^{eq}(c) \phi(c)\, dc$ 替换为格子求和 $\sum_i w_i f_i^{eq} \phi(\mathbf{c}_i)$ ,需要quadrature(数值求积公式)。关键事实是 — 对 $N$ 阶截断而言,只要quadrature能精确积出 $N$ 阶多项式就够了。

在2D中,精确到 $u^2$ 项要求quadrature对Gaussian权重下五阶多项式精确(动量张量会产生 $c^2 \times u^2$ 项)。三点Gauss–Hermite quadrature恰好胜任这件事。

1D下有三个节点 $c \in \{-\sqrt{3}, 0, +\sqrt{3}\}$ ,权 $\{1/6, 2/3, 1/6\}$ 。2D中作张量积得到 $3 \times 3 = 9$ 个节点 — 这就是D2Q9。重新换算到格子单位, $c_s = 1/\sqrt{3}$ ,格子速度向量都是整数坐标 $(0,0), (\pm 1, 0), (0, \pm 1), (\pm 1, \pm 1)$ 。格子权重 $w_i$ 是把Gaussian因子吸收进Gauss–Hermite权后的结果。

编号	$\mathbf{c}_i$	$w_i$
0	$(0, 0)$	$4/9$
1–4	$(\pm 1, 0), (0, \pm 1)$	$1/9$
5–8	$(\pm 1, \pm 1)$	$1/36$

把 $c_s^2 = 1/3$ 代入连续形式,就得到格子上的平衡分布。

f_i^{eq} = w_i\, \rho \left[ 1 + 3 (\mathbf{c}_i \cdot \mathbf{u}) + \frac{9}{2}(\mathbf{c}_i \cdot \mathbf{u})^2 - \frac{3}{2}\, \mathbf{u} \cdot \mathbf{u} \right]

这一行是所有isothermal LBM代码的心脏。表面上只是一个简单多项式,里面却封存着Maxwell–Boltzmann的前三阶矩。

Python — 一次BGK碰撞#

在一个格子节点上跑一步碰撞。

import numpy as np
 
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])
 
def equilibrium(rho, ux, uy):
    """D2Q9平衡 — 9个多项式"""
    ciu = cx * ux + cy * uy            # 点积
    usq = ux * ux + uy * uy
    return w * rho * (1 + 3 * ciu + 4.5 * ciu**2 - 1.5 * usq)
 
def bgk_collision(f, tau):
    rho = f.sum()
    ux  = (cx * f).sum() / rho
    uy  = (cy * f).sum() / rho
    feq = equilibrium(rho, ux, uy)
    return f - (f - feq) / tau, rho, ux, uy
 
# 非平衡初始分布
f = w * 1.0
f[1] += 0.05; f[3] -= 0.04          # 略向东偏
 
for step in range(20):
    f, rho, ux, uy = bgk_collision(f, tau=0.8)
    err = np.linalg.norm(f - equilibrium(rho, ux, uy))
    print(f"step {step:2d}: rho={rho:.4f}, u=({ux:+.4f},{uy:+.4f}), ||delta||={err:.2e}")

跑大约二十步,残差落到 $10^{-6}$ 以下,9个 $f_i$ 沿着平衡多项式各就各位。每步保留下来的是什么 — 密度 $\rho$ 与动量 $\rho \mathbf{u}$ 。碰撞只削去非平衡矩。

D2Q9 격자 위에서 9개 분포 함수가 BGK 충돌을 통해 평형으로 끌려가는 과정. 노란 화살표는 거시 속도 u, 가는 화살표는 각 격자 방향 c_i의 분포 크기.

u_x = 0.12u_y = 0.06τ (이완 시간) = 0.80

τ→0.5에 가까울수록 한 스텝에 평형으로 끌려가지만, 점성 ν = c_s²(τ−1/2)도 작아져 안정성이 흔들린다.

把上方的 $\tau$ 滑到 0.55 附近,一两步内残差就归零。粘性 $\nu = c_s^2 (\tau - 1/2)$ ,当 $\tau \to 1/2$ 时粘性消失,格子立刻失稳。反向把 $\tau$ 调大,粘性涨,松弛速度变慢。这一trade-off是每个LBM从业者都要调的旋钮。

三行带走#

LBM的平衡分布 $f_i^{eq}$ 是Maxwell–Boltzmann在Hermite多项式上的展开,截断到二阶 — 这一截断精确保住三个矩: $\rho$ 、 $\rho \mathbf{u}$ 、动量通量。
D2Q9 推导自 二阶平衡分布所需的最小精度quadrature(Gauss–Hermite三点)的二维张量积。9不是随便的数字,而是正交求积的副产品。
截断阶决定了一切边界。等温假设、 $|\mathbf{u}|/c_s \lesssim 0.3$ 上限,都从同一句话"我们在 $N=2$ 处截断"自然延伸而来。