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

1986年に格子ガスは倒れました。セルあたりの粒子が0個か1個で、ノイズは平均しても消えません。1992年にその整数は実数分布 $f_i$ に替わりました。ノイズは消えました。問いは残ります。実数9個が、無限次元の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）は正方格子、6方向、衝突則一つでした。マクロ平均はNavier–Stokesに戻ります。ブール粒子（0または1）は統計的ノイズを増やしました。時間平均を延ばすとノイズは死にますが、見たかった過渡（transient）も一緒に死にます。

McNamaraとZanettiは各節点に実数分布 $f_i(\mathbf{x}, t)$ を置きました。$f_i$ は時刻 $t$、位置 $\mathbf{x}$ で方向 $\mathbf{c}_i$ へ向かう粒子の期待個数です。1ステップは二つに分かれます。衝突（collision）と伝送（streaming）です。

$\tau$ は緩和時間（relaxation time）です。$f_i^{eq}$ が平衡分布です。衝突は非平衡部分を $\tau$ の尺度で削り、streamingはその分布を1セル動かします。9個で足りる理由は、この $f_i^{eq}$ の形にあります。

Maxwell–Boltzmannが抱える時限爆弾#

連続運動論の平衡はMaxwell–Boltzmannです。

$\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、エルミートです。

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

確率論的エルミート多項式（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$ です。1次元平衡を等温、つまり $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$ の最初の3モーメントは合います。破線が離れる $|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$ で密度・運動量・エネルギーが合っても、捨てた3次の大きさはすでに天井値の2倍を超えています。破線が離れる $u=1$ は、その係数が0.17付近まで上がった場所です。マッハ天井0.3は、曲線が目に見えて離れる点ではありません。3次の係数が0.0045を超える瞬間に、あらかじめ引いた線です。

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

1次粘性の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の天井を同時に残します。その3次を戻すと、等温という仮定が先に壊れます。熱流束はその項が運ぶ断片なので、2次で止めたコードは、温度が流れる問題をこの多項式では受け取れません。マッハの天井と等温仮定は、一つの切断の両面です。速度を0.173より下に置いても、熱い壁と冷たい壁を同じ格子に載せると、この平衡は熱を運びません。

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

連続積分 $\int f^{eq}(c)\,\phi(c)\,dc$ を格子和 $\sum_i w_i f_i^{eq}\phi(\mathbf{c}_i)$ に変えるには求積（quadrature）が要ります。$N$ 次まで合わせるなら、その次数の多項式を正確に積分すれば足ります。

2次元で $u^2$ まで合わせるには、ガウシアン重みの下で5次多項式まで正確でなければなりません。運動量テンソルが $c^2\times u^2$ を生むからです。Gauss–Hermiteの3点求積がその積分を担います。

1次元の節点は $c\in\{-\sqrt{3},0,+\sqrt{3}\}$、重みは $\{1/6,2/3,1/6\}$ です。その三つの根がなぜその位置なのかは、6次モーメントを9速度へ折ったエルミート求積が別に数えます。2次元のテンソル積は $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次元の重み $\{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本と対角4本を足すと $4/9+4\cdot(1/9)+4\cdot(1/36)=1$ で再び閉じます。9という個数と三つの分母は、同じテンソル積から出ます。

この1行が等温LBMコードの平衡です。多項式ですが、Maxwell–Boltzmannの最初の3モーメントを中に閉じ込めています。

Python — BGK衝突の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}")

初期分布は $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$ より下です。本文が言う20ステップ後の $10^{-6}$ 未満の残差は、この保存された $(1.01,\ 0.0891)$ の上の平衡へ寄る過程です。衝突が残すのは、各ステップの $\rho$ と $\rho\mathbf{u}$ だけです。

$\tau$ を0.55側へ下げると、残差は1、2ステップで消えます。粘性は $\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$ で、コードの粘性の6分の1です。$\tau$ をさらに $1/2$ まで押すと粘性は0になり、格子は不安定になります。$\tau$ を上げると粘性は増え、平衡へ引き寄せられる速さは落ちます。非平衡を1ステップで掛ける倍率は $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で、1ステップで平衡に着きます。コードの0.8は、その1ステップ平衡から一つ外した値です。

マッハ0.3の向こうでは、この多項式を信じません#

$N=2$ は $\rho$、$\rho\mathbf{u}$、運動量フラックスまでを合わせます。その上は捨てました。格子速度が約0.173を超えると、捨てた項はもう小さくありません。$\tau$ を下げて粘性を $0.1$ から $1/60$ にしても、平衡多項式はそのままです。速度を天井の上へ上げることと、粘性だけを下げることは、別の失敗です。

静止枠で合わせた粘性が、枠を横へずらしてもそのままかは、この式は答えません。切り捨てた3次がその答えを持っています。枠を動かしたとき、3次モーメントの欠損が粘性を変えるかは、まだこの多項式の外にあります。