9本の矢に閉じ込めたMaxwell–Boltzmann — LBM平衡分布が多項式である理由

エルミート2次切り捨てがD2Q9を作り、同じ切り捨てがLBMを低速流に縛る

1986年、Frisch・Hasslacher・Pomeauが格子の上に粒子を放ってNavier–Stokesを真似させました。LGCA(Lattice Gas Cellular Automaton)です。5年も持ちませんでした — 格子1セルに粒子が0個か1個という整数条件が、どうしても消えないノイズを作り出したからです。1992年、McNamaraとZanettiは整数を実数に置き換えます。粒子の分布関数 $f_i$ で格子を埋める。ノイズは消えました — その代わり、新しい問いが残ります。9個の実数で、無限次元のMaxwell–Boltzmannをどう真似られるのか。本記事はその答えを追います — エルミート多項式の2次切り捨てがLBMの平衡分布を作り、同じ切り捨てがLBMを低速流に縛り付けるという事実までを。

整数から実数へ — 1986年のLGCAの死とLBMの誕生#

LGCAの発想は単純で美しいものでした。正方格子、6方向の粒子、たった1つの衝突則。マクロに平均すればNavier–Stokesが戻ってきます。しかしブール粒子(0または1)は統計的ノイズを増幅しました。時間平均を長く取ればノイズは死にますが、それと同時に、見たかったtransient流れも消えてしまいます。

McNamara・Zanettiは整数粒子の代わりに実数分布 $f_i(\mathbf{x}, t)$ を各格子点に置きます。 $f_i$ は時刻 $t$ 、位置 $\mathbf{x}$ で格子方向 $\mathbf{c}_i$ に向かって動く粒子の 期待個数。ノイズは構造的に消えます。1ステップの進行は、collision(衝突)とstreaming(伝送)の2段階で完結します。

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が削られた分布を格子1セル隣へ運びます。しかし、 $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個の離散方向だけを残すには、どうにかしてこの連続分布を間引かなければなりません。

素朴な試み — ガウシアン積分を矩形和で近似する。失敗します。Navier–Stokesが必要とするモーメント — $\rho$ 、 $\rho \mathbf{u}$ 、運動量テンソル $\rho \mathbf{u}\mathbf{u}$ 、エネルギー — が格子和とずれます。保存が崩れた瞬間にLBMは使い物になりません。

正解は関数展開です。Maxwell–Boltzmannを 直交多項式級数 で展開し、必要なモーメントだけ残してその先を切り捨てる。その直交多項式がHermite(エルミート)です。

エルミート多項式が分布を細切れにする#

確率論的エルミート多項式 $He_n(c)$ は標準ガウシアン重み $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$ です。1次元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$ は静止ガウシアン。 $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$ を選ぶと、最初の3つのモーメント(密度・運動量・エネルギー)がMaxwell–Boltzmannと正確に一致します。 $u$ を 1.0 を超えて押し上げると、破線が一気に外れていきます。これがLBMのマッハ数制約 $|\mathbf{u}|/c_s \lesssim 0.3$ の出どころを、1枚の図で示しています。

2次切り捨ての意味 — 保存されるモーメント#

Navier–Stokesの1次粘性応力まで再現するには、どのモーメントが必要でしょうか。正確な答えは — $\rho$ 、 $\rho \mathbf{u}$ 、そして運動量フラックス $\Pi_{\alpha\beta} = \rho u_\alpha u_\beta + p \delta_{\alpha\beta}$ までの 3つのモーメント次数。 $\rho \mathbf{u}\mathbf{u}$ が $u$ の2次なので、分布関数に $u^2$ 項が入っている必要があります。

したがってエルミート級数を $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の2つの有名な限界を作ります — 等温仮定と、マッハ数 0.3 の天井です。

D2Q9はどうやって作られたか#

連続積分 $\int f^{eq}(c) \phi(c)\, dc$ を格子和 $\sum_i w_i f_i^{eq} \phi(\mathbf{c}_i)$ に置き換えるにはquadrature(直交積分公式)が要ります。核心は — $N$ 次切り捨てまで正確なquadratureは、その次数の多項式を正確に積分すれば十分 という点です。

2次元で $u^2$ 項まで正確にするには、ガウシアン重みに対して5次多項式まで正確なquadratureが必要です(運動量テンソルが $c^2 \times u^2$ 項を生むため)。Gauss–Hermite 3点 quadratureが、ちょうどその仕事をします。

1Dで点が3つ $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$ はGauss–Hermite重みにガウシアン因子を吸収させた値です。

index	$\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]

この1行があらゆるisothermal LBMコードの心臓です。単純な多項式ですが、その中にMaxwell–Boltzmannの最初の3モーメントが封印されています。

Python — BGK衝突の1ステップ#

1つの格子セルで衝突を1ステップ回してみます。

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}")

20ステップほど回すと残差は $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 近くまで下げると、1〜2ステップで残差が0に達します。粘性は $\nu = c_s^2 (\tau - 1/2)$ なので、 $\tau \to 1/2$ では粘性が0になり、格子が不安定になります。逆に $\tau$ を大きくすると粘性は増えますが、平衡へ引き寄せられる速度は遅くなります。このtrade-offがLBMソルバー調整の核心です。

覚えておく3行#

LBMの平衡分布 $f_i^{eq}$ は、Maxwell–Boltzmannをエルミート多項式級数で展開して2次で切った形です — その切り捨てが $\rho$ ・ $\rho \mathbf{u}$ ・運動量フラックスの3モーメントを正確に保存します。
D2Q9格子は、2次切り捨ての平衡分布を正確に積分するGauss–Hermite 3点quadratureの2Dテンソル積から導出されます。9は自然数ではなく、直交積分の産物です。
切り捨て次数が平衡分布の限界を決めます — 等温仮定、マッハ数 $\lesssim 0.3$ の天井は、すべて「 $N=2$ で切った」という1行から自動的に従います。