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

格子气在1986年倒了。每格粒子不是0就是1，噪声靠平均也抹不掉。1992年，那些整数换成了实数分布 $f_i$。噪声没了。问题还在。九个实数凭什么代替无穷维的Maxwell–Boltzmann。

拨动图上的 drift $u$ 和截断阶 $N$，实线的Maxwell–Boltzmann会和虚线多项式叠在一起，或者分开。

$N=2$ 且 $u=0.4$ 时，密度、动量、能量叠在一起。$u$ 超过1，只有虚线离开实线。马赫上限 $|\mathbf{u}|/c_s \lesssim 0.3$ 并不是这个分开的位置。数字留到后面再除。

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

LGCA（Lattice Gas Cellular Automaton）是正方格子、六个方向、一条碰撞规则。宏观平均回到Navier–Stokes。布尔粒子（0或1）把统计噪声放大了。把时间平均拉长，噪声会死，想看的瞬态（transient）也一起死。

McNamara和Zanetti在每个节点放上实数分布 $f_i(\mathbf{x}, t)$。$f_i$ 是时刻 $t$、位置 $\mathbf{x}$ 上沿方向 $\mathbf{c}_i$ 运动的粒子期望个数。一步分成两截。碰撞（collision），然后是传播（streaming）。

$\tau$ 是松弛时间（relaxation time）。$f_i^{eq}$ 是平衡分布。碰撞按 $\tau$ 的尺度削掉非平衡部分，streaming把这个分布推过一格。九个数够用，是因为 $f_i^{eq}$ 的形状。

Maxwell–Boltzmann怀里的定时炸弹#

连续动理论里的平衡是Maxwell–Boltzmann。

$\mathbf{c}$ 是分子速度，$\mathbf{u}$ 是宏观速度，$D$ 是维数，$T$ 是温度，$R$ 是气体常数。炸弹在于 $\mathbf{c}$ 是连续的。只留九个或十九个方向，就得把这个分布削薄。

把高斯积分换成矩形和，会失败。Navier–Stokes要的矩，$\rho$、$\rho\mathbf{u}$、动量张量 $\rho\mathbf{u}\mathbf{u}$、能量，和格子对不上。守恒一破，LBM就停在那里。

留下来的办法是展开。把Maxwell–Boltzmann展成正交多项式，只留需要的矩。那族多项式是Hermite。

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

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

前几项是 $He_0=1$、$He_1=c$、$He_2=c^2-1$、$He_3=c^3-3c$。把一维平衡取成等温，也就是 $RT=1$，生成函数给出下面这一式。

$n=0$ 是静止的高斯。$n=1$ 是小的 drift。$n=2$ 修正平均动能。更高的 $n$ 是不对称的尾巴。LBM把这个和在有限阶截断。

等温展开里 $RT=1$，热速度也是1。图上的 $u$ 已经除以这个热速度。写成上限的0.3，就是这把尺子上的 $|u|$。$u=0.4$ 比上限高出0.1，$N=2$ 的前三个矩仍然重合。虚线分开的 $|u|>1$ 已经过了马赫1。0.3和1不是同一次失败。

级数的 $n=3$ 项是 $(u^3/6)\,He_3(c)$，$He_3=c^3-3c$。只看系数，$u=0.3$ 时 $0.3^3/6=0.027/6=0.0045$。$u=0.4$ 时 $0.4^3/6=0.064/6\approx 0.0107$，大约是上限处的2.37倍。$u=1$ 时 $1/6\approx 0.1667$，大约是0.0045的37倍。这个比就是 $1^3/0.3^3=1/0.027$。$N=2$ 把这一项整个丢掉。所以 $u=0.4$ 时密度、动量、能量仍然重合，丢掉的三阶大小却已经超过上限值的两倍。虚线分开的 $u=1$，是这个系数爬到0.17附近的地方。马赫上限0.3并不是曲线明显离开的那个点。它是三阶系数越过0.0045时事先画下的一条线。

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

一阶粘性Navier–Stokes要的是三层。$\rho$、$\rho\mathbf{u}$，以及动量通量 $\Pi_{\alpha\beta}=\rho u_\alpha u_\beta+p\delta_{\alpha\beta}$。$\rho\mathbf{u}\mathbf{u}$ 是 $u^2$，所以分布里必须有 $u^2$ 项。

级数在 $N=2$ 截断。

这是平衡分布的连续形式。$N=3$ 及以上进热通量。标准等温LBM不留它们。丢掉的项同时留下等温假设和马赫数0.3的上限。把那个三阶加回去，等温这个假设会先破。热通量是该项搬走的那一块，所以停在二阶的代码，没法用这个多项式接一个温度在流动的问题。马赫上限和等温假设是同一次截断的两面。就算把速度放在0.173以下，同一张网格上放一面热壁和一面冷壁，这个平衡也不运热。

D2Q9是如何造出来的#

把连续积分 $\int f^{eq}(c)\,\phi(c)\,dc$ 换成格子求和 $\sum_i w_i f_i^{eq}\phi(\mathbf{c}_i)$，需要求积（quadrature）。要对到 $N$ 阶，把该阶多项式精确积出来就够了。

二维里要对到 $u^2$，高斯权重下的五阶多项式必须精确。动量张量会造出 $c^2\times u^2$。三点Gauss–Hermite求积负责这个积分。

一维节点是 $c\in\{-\sqrt{3},0,+\sqrt{3}\}$，权重是 $\{1/6,2/3,1/6\}$。这三个根为什么落在那里，把六阶矩折进九个速度的Hermite求积另外在数。二维张量积是 $3\times 3=9$。这就是D2Q9。换到格子单位，$c_s=1/\sqrt{3}$，速度落在整数坐标 $(0,0),(\pm 1,0),(0,\pm 1),(\pm 1,\pm 1)$ 上。权重 $w_i$ 是把高斯因子吸进Gauss–Hermite权重之后的数。

格子上的平衡，是连续形式代入 $c_s^2=1/3$ 之后的结果。

$1/c_s^2=3$，$1/(2c_s^4)=9/2$，$1/(2c_s^2)=3/2$。式子里的3、$9/2$、$3/2$ 就是这三个数。$\mathbf{u}=0$ 时括号等于1，$f_i^{eq}=w_i\rho$。权重相加是 $4/9+4\cdot(1/9)+4\cdot(1/36)=1$。静止质量的 $4/9$ 留在0号方向上。表里的三个权重，是一维权重 $\{1/6,\ 2/3,\ 1/6\}$ 的乘积。静止方向是 $(2/3)\times(2/3)=4/9$。轴向是 $(2/3)\times(1/6)=1/9$。对角是 $(1/6)\times(1/6)=1/36$。四个轴向加四个对角，和再次闭上，$4/9+4\cdot(1/9)+4\cdot(1/36)=1$。九这个数目和三个分母，出自同一次张量积。

这一行就是等温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}")

初始分布从 $f=w$ 出发，给 $f_1$ 加上0.05，从 $f_3$ 减去0.04。质量是 $1+0.05-0.04=1.01$。$x$ 方向动量是 $1\cdot 0.05+(-1)\cdot(-0.04)=0.09$。所以 $u_x=0.09/1.01\approx 0.0891$，$u_y=0$。除以 $c_s=1/\sqrt{3}$，马赫数约为0.154。上限0.3对应的格子速度是 $0.3/\sqrt{3}\approx 0.173$，0.0891在它下面。文中说大约二十步之后残差降到 $10^{-6}$ 以下，指的就是贴到这个被守恒的 $(1.01,\ 0.0891)$ 上的平衡。碰撞每一步留下的只有 $\rho$ 和 $\rho\mathbf{u}$。

把 $\tau$ 拨向0.55，残差一两步就没了。粘性是 $\nu=c_s^2(\tau-1/2)$。代码里的 $\tau=0.8$ 给出 $\nu=(1/3)\cdot 0.3=0.1$。$\tau=0.55$ 给出 $\nu=(1/3)\cdot 0.05=1/60\approx 0.0167$，是代码粘性的六分之一。再把 $\tau$ 推到 $1/2$，粘性变成0，格子失稳。把 $\tau$ 加大，粘性上升，被拉向平衡的速度变慢。非平衡部分每一步乘上 $1-1/\tau$。代码里的 $\tau=0.8$ 给出 $1/\tau=1.25$，倍率就是 $-0.25$。符号会翻过来，这是过松弛。绝对值 $0.25$ 乘20次是 $0.25^{20}\approx 9.1\times 10^{-13}$。偏差就算是1的量级，20步之后也远低于 $10^{-6}$。文中把20步当作 $10^{-6}$ 的刻度，单凭这个倍率已经留了余量。$\tau=1$ 时倍率为0，一步就落到平衡上。代码里的0.8，是从那一步平衡挪开一格的值。

过了马赫数0.3，就不要信这个多项式#

$N=2$ 只对准 $\rho$、$\rho\mathbf{u}$ 和动量通量。再往上的都丢了。格子速度一过大约0.173，丢掉的项就不再小。把 $\tau$ 降下，让粘性从 $0.1$ 变到 $1/60$，平衡多项式本身不动。把速度抬过上限，和只把粘性降下，是两次不同的失败。

静止参考系里调好的粘性，参考系横移之后还在不在，这个式子答不了。被截掉的三阶拿着那个答案。参考系一动，缺失的三阶矩会不会改掉粘性仍然在这个多项式外面。