BGKが発散する低粘性域を生き抜く方法 — MRT格子ボルツマンの衝突

s_νだけを1.99へ上げて粘性を0付近に置くと、残ったs_eとs_qのバーがチェッカーボードを分けます。

同じ格子、同じ平衡なのに、一方は滑らかに回り、もう一方はチェッカーボードで破綻します。BGKは9個の分布関数を一つの緩和時間で押します。MRT(Multiple Relaxation Time、多重緩和時間)はモーメントに変えてからそれぞれ別の速さで押します。

$s_\nu$ のバーを1.99まで押すと、粘性の数字が0の方へ落ちます。 $s_e$ と $s_q$ のバーはその場に残ります。

D2Q9 moment-space relaxation rates · S = diag(s₀…s₈)

sₑ (bulk) = 1.64

s_ε = 1.54

s_q (ghost) = 1.70

s_ν (shear) = 1.20

shear viscosity ν = cs²(1/s_ν − 1/2) = 0.1111 (lattice units)

✓ all rates in (0, 2) — linearly stable region

保存量のバー $\rho$ , $j_x$ , $j_y$ は0に付いたままです。BGKはこのバーを一つの値に縛ります。低粘性でチェッカーボードが出る場所がその束ねです。

すべてを一拍で緩和させるBGKの弱点#

格子ボルツマン法(LBM)は、分布関数 $f_i$ が格子上を移動(streaming)し、各点で平衡へ緩和(collision)する2段階を繰り返します。最も一般的なBGK衝突はこう書けます。

f_i(\mathbf{x}+\mathbf{c}_i,\, t+1) = f_i(\mathbf{x}, t) - \frac{1}{\tau}\bigl(f_i - f_i^{\text{eq}}\bigr)

ここで $f_i$ は速度 $\mathbf{c}_i$ 方向の分布関数、 $\tau$ は緩和時間、 $f_i^{\text{eq}}$ は局所平衡分布です。粘性は $\tau$ 一つで決まります。

\nu = c_s^2\left(\tau - \tfrac{1}{2}\right), \qquad c_s^2 = \tfrac{1}{3}

問題はここから始まります。粘性を下げる(レイノルズ数を上げる)には $\tau$ を $1/2$ に近づける必要があります。ところが $\tau \to 1/2$ になると、粘性と何の関係もない高次モードまでほとんど緩和されずに格子を動き回ります。これらのモードが格子解像度の限界で振幅を増すと、チェッカーボード振動が生じ、解析は発散します。BGKでは粘性の調整と安定性の調整が一つのつまみに縛られています。

衝突をモーメント空間で眺める

9個の分布関数をそのまま緩和させません。まず物理的な意味を持つモーメントに変換してから緩和させます。モーメントは分布関数の線形結合です。

\mathbf{m} = M\,\mathbf{f}

$\mathbf{f} = (f_0,\dots,f_8)^\top$ は分布関数ベクトル、 $M$ は $9\times 9$ の変換行列、 $\mathbf{m}$ はモーメントベクトルです。衝突はモーメント空間で平衡モーメント $\mathbf{m}^{\text{eq}}$ へ緩和させ、再び速度空間へ戻します。

\mathbf{f}^{*} = \mathbf{f} - M^{-1}\,S\,\bigl(M\mathbf{f} - \mathbf{m}^{\text{eq}}\bigr)

$S = \mathrm{diag}(s_0,\dots,s_8)$ はモーメントごとの緩和率を並べた対角行列です。すべての $s_i$ を同じ値 $1/\tau$ にすれば、この式は正確にBGKへ戻ります。MRTはその対角成分をモーメントごとに別々に解くだけです。

D2Q9の9つのモーメントと変換行列M#

D2Q9格子(2次元、9速度)でのモーメントは次の9つです。密度 $\rho$ 、エネルギー $e$ 、エネルギー二乗 $\varepsilon$ 、運動量 $j_x, j_y$ 、エネルギーフラックス $q_x, q_y$ 、そして応力に対応する $p_{xx}, p_{xy}$ です。

このうち $\rho, j_x, j_y$ は衝突で変化しない保存量なので、緩和率を $s=0$ に固定します。残りの6つが調整可能なつまみです。Lallemand–Luoの標準変換行列は各行が直交するように設計されており、逆行列はノルムで割った転置として簡単に得られます。

平衡モーメントは密度と運動量だけで閉じます。

\begin{aligned} e^{\text{eq}} &= -2\rho + 3(j_x^2 + j_y^2), &\quad \varepsilon^{\text{eq}} &= \rho - 3(j_x^2 + j_y^2) \\ q_x^{\text{eq}} &= -j_x, &\quad q_y^{\text{eq}} &= -j_y \\ p_{xx}^{\text{eq}} &= j_x^2 - j_y^2, &\quad p_{xy}^{\text{eq}} &= j_x j_y \end{aligned}

$j_x = \rho u_x$ 、 $j_y = \rho u_y$ は運動量成分です。応力モーメント $p_{xx}, p_{xy}$ の平衡が速度勾配と結びつくため、これらの緩和率が粘性を決めます。

モーメントごとに異なる緩和率 — 粘性と無関係な自由度

せん断粘性を決めるのは応力モーメントの緩和率 $s_\nu$ です。

\nu = c_s^2\left(\frac{1}{s_\nu} - \frac{1}{2}\right)

残りの緩和率 $s_e$ (体積)、 $s_\varepsilon$ 、 $s_q$ (エネルギーフラックス)は粘性に入りません。自由パラメータです。 $s_\nu$ を $2$ に寄せても、高次モードは安定区間 $1.5\sim1.8$ に別に置けます。BGKが一つのつまみに縛っていたものが、ここで分かれます。

コードの既定粘性は $\nu=0.004$ です。 $s_\nu = 1/(3\nu + 1/2)$ なので $3\times 0.004 + 0.5 = 0.512$ 、 $s_\nu = 1/0.512 = 1.953125$ です。先のバーを1.99まで押すと $\nu = (1/3)(1/1.99 - 1/2)$ です。 $1/1.99 \approx 0.5025$ から0.5を引き、3で割ると $\nu \approx 8.4\times 10^{-4}$ です。0に近いだけで0ではありません。その粘性でBGKは $s_e$ 、 $s_\varepsilon$ 、 $s_q$ まで1.99に一緒に押します。コードのMRTはその三席を $1.64$ 、 $1.54$ 、 $1.7$ に固定します。回る配列は [0, 1.64, 1.54, 0, 1.7, 0, 1.7, s_nu, s_nu] で、0は保存量の三席だけです。

平衡モーメントが密度と運動量だけで閉じるのは、分布を低い次数で切った後の話です。その切り方は 9速度の中のMaxwell–Boltzmann の側にあります。

衝突はモーメント空間、ストリーミングは速度空間

アルゴリズムの流れは次の通りです。

分布 $f_i$ から巨視量 $\rho, u_x, u_y$ を計算
$\mathbf{m} = M\mathbf{f}$ でモーメントへ変換
平衡モーメント $\mathbf{m}^{\text{eq}}$ を計算(密度・運動量から)
モーメントごとの緩和: $m_a \leftarrow m_a - s_a(m_a - m_a^{\text{eq}})$
$\mathbf{f}^{*} = M^{-1}\mathbf{m}$ で速度空間へ復元
ストリーミング: $f_i^{*}$ を $\mathbf{c}_i$ 方向の隣へ移動
境界条件を適用して1へ

衝突だけがモーメント空間で起こり、ストリーミングはそのまま速度空間で進みます。追加コストはセルごとの $9\times 9$ 行列積2回だけです。

Python — Taylor–Green渦で粘性減衰を測る#

検証にはTaylor–Green渦を使います。この流れは粘性減衰の厳密解が知られているので、実装が正しいかを数値で確認できます。運動エネルギーは $\exp(-4\nu k^2 t)$ で減るはずです( $k$ は渦の波数)。

import numpy as np
 
# D2Q9 格子
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])
 
# Lallemand-Luo 変換行列 (行: rho,e,eps,jx,qx,jy,qy,pxx,pxy)
M = np.array([
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [-4,-1,-1,-1,-1, 2, 2, 2, 2],
    [ 4,-2,-2,-2,-2, 1, 1, 1, 1],
    [ 0, 1, 0,-1, 0, 1,-1,-1, 1],
    [ 0,-2, 0, 2, 0, 1,-1,-1, 1],
    [ 0, 0, 1, 0,-1, 1, 1,-1,-1],
    [ 0, 0,-2, 0, 2, 1, 1,-1,-1],
    [ 0, 1,-1, 1,-1, 0, 0, 0, 0],
    [ 0, 0, 0, 0, 0, 1,-1, 1,-1],
], dtype=float)
Minv = np.linalg.inv(M)
 
def d2q9_equilibrium(rho, ux, uy):
    eq = np.zeros((9,) + rho.shape)
    for i in range(9):
        cu = cx[i] * ux + cy[i] * uy
        eq[i] = w[i] * rho * (1 + 3*cu + 4.5*cu**2 - 1.5*(ux**2 + uy**2))
    return eq
 
def relaxation_rates(nu):
    s_nu = 1.0 / (3*nu + 0.5)          # 粘性を決める緩和率
    # 保存(rho,jx,jy)=0、高次モードは安定区間の固定値
    return np.array([0, 1.64, 1.54, 0, 1.7, 0, 1.7, s_nu, s_nu])
 
def taylor_green_init(N, U0):
    x = np.arange(N)
    X, Y = np.meshgrid(x, x, indexing='ij')
    k = 2*np.pi / N
    ux = -U0 * np.cos(k*X) * np.sin(k*Y)
    uy =  U0 * np.sin(k*X) * np.cos(k*Y)
    rho = 1 - 0.75 * U0**2 * (np.cos(2*k*X) + np.cos(2*k*Y))
    return d2q9_equilibrium(rho, ux, uy), k
 
def macros(f):
    rho = f.sum(axis=0)
    ux = (cx[:, None, None] * f).sum(axis=0) / rho
    uy = (cy[:, None, None] * f).sum(axis=0) / rho
    return rho, ux, uy
 
def mrt_collide(f, s):
    rho, ux, uy = macros(f)
    jx, jy = rho*ux, rho*uy
    m = np.tensordot(M, f, axes=([1], [0]))        # モーメント空間へ
    meq = np.zeros_like(m)
    meq[0] = rho
    meq[1] = -2*rho + 3*(jx**2 + jy**2)
    meq[2] = rho - 3*(jx**2 + jy**2)
    meq[3], meq[4] = jx, -jx
    meq[5], meq[6] = jy, -jy
    meq[7] = jx**2 - jy**2
    meq[8] = jx*jy
    m -= s[:, None, None] * (m - meq)               # モーメントごとの緩和
    return np.tensordot(Minv, m, axes=([1], [0]))    # 速度空間へ復元
 
def stream(f):
    for i in range(9):
        f[i] = np.roll(f[i], (cx[i], cy[i]), axis=(0, 1))
    return f
 
def run_mrt_lbm(N=64, nu=0.004, U0=0.04, steps=2000):
    f, k = taylor_green_init(N, U0)
    s = relaxation_rates(nu)
    e0 = None
    for t in range(steps):
        f = mrt_collide(f, s)
        f = stream(f)
        if t % 400 == 0:
            _, ux, uy = macros(f)
            E = np.mean(ux**2 + uy**2)
            e0 = E if e0 is None else e0
            print(f"t={t:5d}  KE/KE0={E/e0:.4f}  analytic={np.exp(-4*nu*k**2*t):.4f}")
 
run_mrt_lbm()

printが並べる解析値は $\exp(-4\nu k^2 t)$ です。既定は $N=64$ 、 $\nu=0.004$ 、 $k=2\pi/N$ なので $4\nu k^2 = 0.016\,(\pi/32)^2 \approx 1.54\times 10^{-4}$ です。 $t=400$ では指数が約 $-0.0617$ 、解析因子は約0.940です。 $t=1600$ では約 $-0.247$ 、因子は約0.781です。range(2000) なので最後に出る $t$ は1600です。小数点以下2桁まで合うというのは、この因子と数値の運動エネルギーの比較です。同じ格子・同じ $\nu$ でBGKを回すと、より低い粘性で先に発散します。

BGKとMRTを同じReで並べる#

粘性スライダーを下げると減衰が遅くなり、BGK/MRTボタンは同じ $\nu$ で格子が持つかを分けます。

viscosity ν = 0.0060 (Re ≈ 267)

mode: MRT

step: 0

KE / KE₀: 1.000

analytic e⁻⁴ᵛᵏ²ᵗ: 1.000

Color = vorticity (red ↺ / blue ↻). Taylor–Green vortex on a 40×40 periodic lattice.

粘性を $0.002$ あたりまで下げてBGKのままにすると、格子に粗い模様が入り「diverged」が表示されます。同じ粘性でMRTに切り替えてresetすると、渦は滑らかに減衰します。 $\nu=0.002$ では $s_\nu=1/(3\times 0.002+0.5)=1/0.506\approx 1.976$ です。BGKは高次モードまでこの値で押します。コードの $1.64$ 、 $1.54$ 、 $1.7$ はその横に残ります。

s_νだけを2に寄せると高次モードが残る#

粘性のつまみ $s_\nu$ と安定性のつまみ $s_e$ 、 $s_q$ を物理的に分離したのがMRTの中身です。 $\nu=0.004$ の $s_\nu=1.953$ を全モードにコピーするとBGKです。 $1.64$ 、 $1.54$ 、 $1.7$ を残して $s_\nu$ だけを上げるのがMRTです。追加コストはセルあたり $9\times 9$ の行列積2回です。

その分離がないと、低粘性のTaylor–Greenはチェッカーボードで死にます。エントロピー制限で同じ暴走を止める別路は H定理を衝突に入れた格子ボルツマンです。減衰の確認は $\exp(-4\nu k^2 t)$ の1行で足ります。