エントロピーが暴走を止める — Karlinのエントロピック格子ボルツマン

エントロピーを失われた情報のラベルだけで覚えると、この計算がなぜ壊れないのか説明がつきません。格子ボルツマン法(Lattice Boltzmann Method, LBM)で粘性を下げるために $\tau$ を $1/2$ へ寄せると、BGKは平衡を通り過ぎます。分布は負に落ちます。Karlin系のエントロピックLBM(Entropic LBM, ELBM)は、その場で離散Hを衝突前と同じ値に縛ります。代償はノードごとの非線形な根ひとつです。

非平衡摂動のスライダーを右へ動かすと、黄色い点が破線 $\alpha=2$ から離れ、$G(\alpha)=0$ の第二の根が左へ動きます。

摂動が小さいうち、点は2の上に張り付いています。破線と黄色い点の間隔が、BGKが捨てない過緩和の幅です。

τが0.5に近づくとBGKは崩れる#

標準的なLBMはBGK衝突(単一緩和時間近似)を使います。分布 $f_i$ を平衡 $f_i^{eq}$ の方へ $\omega$ だけ引き寄せます。

$f_i$ は格子方向 $i$ の分布、$f_i^{eq}$ は平衡分布、$\omega = 1/\tau$ は緩和率です。

粘性は $\nu = c_s^2(\tau - 1/2)$ で決まります。$c_s^2 = 1/3$ は格子音速の二乗です。高Reynolds、つまり低粘性の流れでは $\tau \to 1/2$、すなわち $\omega \to 2$ まで押し込みます。

$\omega$ が2に近づくほど、衝突は平衡の反対側に着地します。一度負になった $f_i$ は次のステップでより大きく揺れます。格子が壊れます。BGKにはこの過緩和(over-relaxation)を止める仕組みがありません。

ボルツマンのH定理 — 時間の矢#

1872年、ボルツマンは次の量が衝突だけでは決して増えないことを証明しました。

$f$ は速度分布、$\mathbf{c}$ は分子速度です。$dH/dt \le 0$ で、等号は平衡(Maxwell-Boltzmann分布)だけで成り立ちます。

Hは負のエントロピーです。Hの単調減少は、エントロピーの単調増加と同じ文です。不可逆性そのものです。時間の矢でもあります。この不等式は安定性の保証書です。Hが増えないスキームでは振幅が際限なく増えません。増えるにはHが一緒に増える必要があるからです。

離散H関数とエントロピック平衡#

LBMは連続速度を9方向(D2Q9)に切ります。Hもそれに合わせて離散化します。

$w_i$ は格子重みです。D2Q9では中心 $4/9$、面 $1/9$、対角 $1/36$ です。重みで割った比の対数です。

エントロピック平衡は、この $H$ を保存量の制約のもとで最小化した分布です。質量 $\sum_i f_i = \rho$ と運動量 $\sum_i \mathbf{c}_i f_i = \rho\mathbf{u}$ を固定したままHを下げます。ラグランジュ乗数の問題です。D2Q9には閉じた積(product)形があります。

$u_\alpha$ は速度成分、$c_{i\alpha}$ は格子速度成分(−1, 0, 1)です。標準LBMの2次切断はHermite展開を2次で切った多項式平衡であり、Hの最小点ではありません。最小点は上の積です。

エントロピック安定化子 — αを毎ステップ解く

ELBMの衝突はBGKと同じ形です。緩和率だけを二つに分けます。

$\beta \in (0,1)$ が粘性を決めます。$\nu = c_s^2(\tfrac{1}{2\beta} - \tfrac{1}{2})$ で、$\beta\to1$ なら無粘性です。$\alpha$ は安定性の側です。

$\alpha$ は定数ではありません。毎ステップ、各格子点で次式の非自明な根として解きます。

衝突後のHを衝突前と正確に同じにする距離だけ、平衡へ跳ばせます。Hが増える余地はありません。$\alpha=0$ は自明な根です。要るのは平衡の向こうにある第二の根です。平衡の近くでは $\alpha \approx 2$ となり、BGKと一致します。非平衡が大きいときだけ $\alpha$ が2から外れ、過緩和を削ります。

粘性の式を数値にすると、過緩和が壁にどれだけ張り付いているかが見えます。$c_s^2 = 1/3$、$\nu = c_s^2(\tau - 1/2)$ です。$\tau = 0.51$ なら差は $0.01$ なので $\nu = 1/300 \approx 0.00333$ です。緩和率は $\omega = 1/0.51 \approx 1.961$ です。平衡から外れた長さ1を1.961だけ進むので、平衡を0.961過ぎた反対側に着きます。$\tau = 0.501$ まで寄せると $\omega \approx 1.996$ です。鏡像 $\omega = 2$ までの残りは0.004だけです。

コードに書いた $\beta = 0.95$ も同じ関係です。$\tfrac{1}{2\beta}-\tfrac{1}{2} = \tfrac{1}{1.9}-\tfrac{1}{2} = \tfrac{1}{38}$ なので $\nu = 1/114 \approx 0.00877$ です。同じ粘性のBGKなら $\tau = 1/(2\beta) \approx 0.526$ です。平衡付近で $\alpha \approx 2$ なら有効緩和率は $\alpha\beta = 1.9$ で、このBGKと粘性が一致します。非平衡が大きくなって $\alpha$ が2を下回ると、$\alpha\beta$ も1.9より小さくなります。平衡を通り越す長さ $\alpha\beta-1$ は0.9より短くなります。粘性は $\beta$ が決めたままです。削られるのは過緩和だけです。

Python — エントロピックαを二分法で#

積形の平衡、H関数、そして二分法で $\alpha$ を解く一ステップです。

import numpy as np
 
# D2Q9 格子重みと速度 (lattice weights & velocities)
W  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
CX = np.array([0, 1, 0, -1,  0, 1, -1, -1,  1])
CY = np.array([0, 0, 1,  0, -1, 1,  1, -1, -1])
 
def h_function(f):
    # 離散H関数: H = Σ f_i ln(f_i / w_i)  (ボルツマンHの格子版)
    return np.sum(f * np.log(f / W))
 
def entropic_equilibrium(rho, ux, uy):
    # 正確なAnsumali–Karlin積(product)形エントロピック平衡
    def factor(u):
        s = np.sqrt(1 + 3 * u * u)
        return (2 - s), (2 * u + s) / (1 - u)
    bx, rx = factor(ux)
    by, ry = factor(uy)
    return rho * W * bx * by * (rx ** CX) * (ry ** CY)
 
def entropic_alpha(f, feq):
    # H(f + α(feq − f)) = H(f) の非自明な根 α* を二分法で
    d = feq - f
    H0 = h_function(f)
    G = lambda a: h_function(f + a * d) - H0
    amax = np.min(np.where(d < 0, -f / d, np.inf))   # 正値を保証する上限
    hi = min(0.999 * amax, 4.0) if np.isfinite(amax) else 4.0
    if G(hi) <= 0:
        return 2.0                                   # 平衡近傍の安全値
    lo = 1.0
    while G(lo) > 0 and lo > 1e-3:
        lo *= 0.5
    for _ in range(60):
        mid = 0.5 * (lo + hi)
        if G(mid) > 0: hi = mid
        else:          lo = mid
    return 0.5 * (lo + hi)
 
# --- 一ノード、一ステップ ---
feq = entropic_equilibrium(1.0, 0.05, 0.0)
f = feq.copy()
f[1] += 0.03; f[3] -= 0.03            # streaming が運んできた非平衡 (non-equilibrium)
rho = f.sum(); ux = (CX * f).sum() / rho; uy = (CY * f).sum() / rho
feq = entropic_equilibrium(rho, ux, uy)
 
alpha = entropic_alpha(f, feq)
beta  = 0.95                          # ν = c_s²(1/(2β) − 1/2), β→1 で低粘性
f_post = f + alpha * beta * (feq - f)
 
print(f"エントロピック α = {alpha:.4f}   (BGK固定値 2.0)")
print(f"H 衝突前         = {h_function(f):.8f}")
print(f"H 衝突後         = {h_function(f_post):.8f}")
print(f"最小分布         = {f_post.min():.3e}  (> 0)")

出力の $\alpha$ は2に近いものの、正確に2ではありません。その隙間がHを保存する補正量です。

摂動の二行はモーメントとして読めます。$f_1$ に $0.03$ を足し、$f_3$ から $0.03$ を引いても質量の和は変わりません。方向1と3の格子速度は $+1$ と $-1$ なので、$x$ 運動量は $0.03-(-0.03)=0.06$ だけ増えます。初めは $\rho=1$、$u_x=0.05$ だったので、摂動の直後は $u_x=0.11$ です。積の因子は $(2u+\sqrt{1+3u^2})/(1-u)$ です。$u=0.05$ で約1.162、$u=0.11$ で約1.391です。無次元速度0.3では約2.467まで上がり、分母 $1-u$ は0.7です。$|u| \lesssim 0.1$ は、この因子がまだ緩い側を選んだ範囲です。

二分法の上限も符号から来ます。$d_i$ が負の成分が0になる位置は $\alpha=-f_i/d_i$ で、その最小が $a_{\max}$ です。探索は $\min(0.999\,a_{\max},\,4)$ で切れます。4を超える根や、分布の符号を変える根は候補に入りません。

インタラクティブ — G(α)の第二の根#

先の曲線がこの $G$ です。黄色い点は非自明な根 $\alpha^*$ です。摂動を大きくすると2から離れ、BGKはその位置を2に固定します。差がステップごとに積もると格子は発散します。

インタラクティブ — Hが増える瞬間を防ぐ#

同じ非平衡状態を二つの衝突則に並べて入れると、負になった分布だけが赤い棒に変わります。

エントロピック側は常に $H_{ent} = H_{pre}$ です。そう組んであります。BGK側は非平衡が大きくなると、$H_{BGK}$ が衝突前のHを超える瞬間が来ます。そこで赤い警告が点きます。そのHの増加が、低粘性で格子が裂ける起点です。

実装でぶつかった三つのこと

全ノードで非線形求解を回すと、粘性項より先にコストが重くなります。平衡の近く、$\|f - f^{eq}\|$ が小さいときは $\alpha \approx 2$ とし、非平衡が閾値を超えた格子だけ二分法を起動します。コストはおおよそ一割まで下がります。

積形平衡は $|u| < 1$ でしか定義されません。分母に $1 - u$ があるからです。無次元速度が0.3を超えると、圧力テンソルが方向ごとに歪み始めます。ELBMでも $|u| \lesssim 0.1$ に置く方が安全です。

2次切断の多項式を平衡のまま使うと、H最小という保証が抜けます。安定化子の利得を全部取るには、平衡そのものがエントロピック積でなければなりません。多項式の上に $\alpha$ だけ載せると半分です。

βだけ上げてαを2のままにするとHの上限が外れる#

平衡付近で $\alpha \approx 2$ と飛ばすのは、$\|f - f^{eq}\|$ が小さいときのコスト削減です。$\beta$ を1へ寄せながらその省略を全ノードに固定すると、有効緩和率 $\alpha\beta$ は再び2に張り付きます。鏡像までの過緩和が戻り、離散Hは衝突前より大きくなり得ます。粘性だけ下げて、根は解いていない選択です。

衝突がHを平らに押さえても、隣のノードで既に負になった分布をストリーミングが運んでくると、その根は遅れます。符号が移流の段階で壊れるとき、正値を保つフラックス制限がまだ要るのか、エントロピー条件だけで足りるのかは、次の分かれ目として残っています。