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

エントロピーはふつう「失われた情報」のラベルで習います。無秩序が増える向き、不可逆な損失の尺度です。ところが格子ボルツマン法(Lattice Boltzmann Method, LBM)では、まさにそのエントロピーが暴走する計算を捕まえて止めます。この記事では、Karlin系のエントロピックLBM(Entropic LBM, ELBM)がどのように離散化の段階でH定理を強制し、低粘性領域の発散を防ぐのかを扱います。読み終えると、衝突演算子の一行のために毎ステップ非線形方程式を解くコストを払う理由が腑に落ちます。

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

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

f_i^* = f_i + \omega\,(f_i^{eq} - f_i)

ここで $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年、ボルツマンは次の量が衝突だけでは決して増えないことを証明しました。

H = \int f \ln f \, d\mathbf{c}

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

Hは負のエントロピーです。Hが単調減少することは、エントロピーが単調増加することと同じです。これが不可逆性、すなわち時間の矢です。要点は、この不等式が安定性の保証書だということです。Hが決して増えない数値スキームは発散できません。振幅が無限に増えるにはHが増える必要があるからです。

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

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

H(f) = \sum_i f_i \ln\frac{f_i}{w_i}

$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)形が存在します。

f_i^{eq} = \rho\, w_i \prod_{\alpha=x,y}\left(2 - \sqrt{1+3u_\alpha^2}\right)\left(\frac{2u_\alpha + \sqrt{1+3u_\alpha^2}}{1 - u_\alpha}\right)^{c_{i\alpha}}

$u_\alpha$ は速度成分、 $c_{i\alpha}$ は格子速度成分(−1, 0, 1)です。標準LBMが使う2次切断多項式平衡と違い、これはHの真の最小点です。

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

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

f_i^* = f_i + \alpha\beta\,(f_i^{eq} - f_i)

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

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

H\big(f + \alpha(f^{eq} - f)\big) = H(f)

意味は明快です。衝突後のHを衝突前と正確に同じにするぶんだけ平衡へ跳ぶ。Hが増えることは構造的に禁じられます。 $\alpha=0$ は自明な根で、求めるのは平衡の向こうの第二の根です。平衡近傍では $\alpha \approx 2$ となり、BGKと一致します。非平衡が大きいときだけ $\alpha$ が2から外れ、過緩和を削ります。

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を保存する補正量です。

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

下のシミュレーションを直接操作してみましょう。非平衡摂動を大きくすると曲線がどう曲がるかを見ます。

Non-equilibrium perturbation: 0.060

The yellow dot is the non-trivial root α* where the discrete H-function returns to its pre-collision value. Near equilibrium α* ≈ 2, which is exactly standard BGK. Push the perturbation up and α* drifts away from 2 — that gap is the entropic correction BGK ignores.

黄色い点が求める $\alpha^*$ です。摂動が小さいときは破線( $\alpha=2$ )のすぐ上に貼り付いています。スライダーを右へ押すと $\alpha^*$ が2から離れます。BGKはこの移動を無視し、つねに2を使います。差が積み重なると発散します。

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

同じ非平衡状態を二つの衝突則で処理した結果を並べて見ます。赤い棒は負になった分布です。

Non-equilibrium shear: 0.070

Both panels show post-collision populations. The entropic collision lands exactly on the equal-entropy point, so H_ent = H_pre. Fixed-α BGK misses it; once the non-equilibrium part grows, BGK over-shoots and H actually increases — a violated discrete H-theorem, which is the seed of low-viscosity blow-up.

エントロピック側はつねに $H_{ent} = H_{pre}$ です。そう設計されています。BGK側は非平衡が大きくなると $H_{BGK}$ が衝突前を上回る瞬間が来ます。そのとき赤い警告が点きます。まさにそのHの増加が、低粘性で格子を破裂させる種です。

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

第一に、各格子点での非線形解法は高価です。そこで実務では平衡近傍( $\|f - f^{eq}\|$ が小さいとき)は $\alpha \approx 2$ とし、非平衡が閾値を超える格子だけ二分法を起動します。コストは約一割に下がります。

第二に、積形平衡は $|u| < 1$ でしか定義されません。分母に $1 - u$ があります。無次元速度が0.3を超えると圧力テンソルが非等方的に歪み始めるので、ELBMでも $|u| \lesssim 0.1$ の領域で使うのが安全です。

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

3行まとめ#

BGKが低粘性で破綻するのは、過緩和が離散Hを増やすからです。ELBMは衝突後のHを衝突前と等しくする $\alpha$ を毎ステップ解いてこれを防ぎます。
$\alpha$ は $H(f + \alpha(f^{eq}-f)) = H(f)$ の非自明な根で、平衡近傍では $\alpha \approx 2$ として自動的にBGKを復元します。
粘性は $\beta$ が、安定性は $\alpha$ が分担します。コストは毎ステップの非線形解法ですが、平衡近傍の短絡で実務では十分に許容できます。