フレームを0.2で押したら粘性が5.97%減った — LBM平衡分布の3次モーメント欠損

LBMの速度上限を決めるのは安定性ではなく、平衡分布が合わせられるモーメントの次数です。

静止した水と流れる水で粘性は同じであるべきか

同じ流体を二度測るとします。一度は静止した箱の中で、もう一度はその箱を一定速度で押しながら。粘性係数は二度とも同じ値になるはずです。ガリレイ不変性、つまり等速で動く二つの座標系で物理法則が同じ形になるという性質がそう要求します。

D2Q9の格子ボルツマン法(LBM)でこの実験をそのまま行うと、二度目の値が5.97%小さく出ます。格子を細かくしても変わりません。緩和時間 $\tau$ を動かしてもこの比率はほとんど動きません。この記事では、その5.97%を平衡分布のモーメント一行まで追いかけます。結論を先に言えば、安定性の問題ではなく代数の問題です。格子速度の集合が、欠けている項をそもそも作れないのです。

格子上の一行をTaylor展開すると#

LBMが解く式は一行です。

f_i(\mathbf{x} + \mathbf{c}_i \Delta t,\; t + \Delta t) - f_i(\mathbf{x}, t) = -\frac{\Delta t}{\tau}\left[ f_i - f_i^{eq} \right]

$f_i$ は格子速度 $\mathbf{c}_i$ に乗る分布関数、 $f_i^{eq}$ はその場の平衡分布、 $\tau$ は緩和時間です。

左辺を $\Delta t$ で展開すると、1次の項だけでなく2次の項が残ります。

\Delta t \left( \partial_t + \mathbf{c}_i \cdot \nabla \right) f_i + \frac{\Delta t^2}{2} \left( \partial_t + \mathbf{c}_i \cdot \nabla \right)^2 f_i + O(\Delta t^3) = -\frac{\Delta t}{\tau} \left( f_i - f_i^{eq} \right)

この2次の項が、Chapman–Enskog展開(分布関数をクヌーセン数のべき級数に分解する多重スケール手法)で有名な $\tau - 1/2$ を生みます。そして同じ展開の次の段階で、得られる連続体方程式がNavier–Stokesになるための条件が出てきます。平衡分布が四つのモーメントを正確に満たすことです。

\sum_i f_i^{eq} = \rho, \qquad \sum_i f_i^{eq} c_{i\alpha} = \rho u_\alpha, \qquad \sum_i f_i^{eq} c_{i\alpha} c_{i\beta} = \rho u_\alpha u_\beta + \rho c_s^2 \delta_{\alpha\beta}

\sum_i f_i^{eq} c_{i\alpha} c_{i\beta} c_{i\gamma} = \rho c_s^2 \left( u_\alpha \delta_{\beta\gamma} + u_\beta \delta_{\gamma\alpha} + u_\gamma \delta_{\alpha\beta} \right) + \rho u_\alpha u_\beta u_\gamma

前の三つが質量・運動量・圧力を決めます。四番目の3次モーメントは、粘性応力の時間微分を置き換えるときに使われます。ここがずれると、連続の式もEulerレベルも健全なまま、粘性項だけが汚染されます。

下のシミュレーションで3次モーメントが合っているかを直接確かめてみましょう。

velocity u : 0.180 (Ma = 0.31)

Σ f c³ = 0.180000  |  Maxwell = 0.185832  |  gap = 0.005832  |  gap / u³ = 1.0000

速度スライダーを動かすと、黄色の曲線(Maxwell分布が要求する値)と青の直線(D2Q9が実際に出す値)が離れていきます。sweep u で自動往復させると、開く幅が $u$ の奇数乗をたどることが見えます。speeds ±2 を押すと棒が五本に増え、隙間がゼロに閉じます。

$c_x^3 = c_x$ — 格子が作れない項#

D2Q9の $x$ 成分は $-1, 0, 1$ の三つだけです。この三つはどれも三乗すると自分自身に戻ります。したがって任意の分布 $f_i$ に対して

\sum_i f_i c_{ix}^3 = \sum_i f_i c_{ix} = \rho u_x

が常に成り立ちます。平衡分布をどう設計しても同じです。九つの枠に乱数を入れても成り立ちます。3次モーメントはすでに1次モーメントに縛られています。

一方、Maxwell分布が要求する値は $\rho(u_x^3 + 3 c_s^2 u_x)$ で、D2Q9では $c_s^2 = 1/3$ なので $3c_s^2 = 1$ です。要求値は $\rho(u_x + u_x^3)$ に畳まれます。二つの差はちょうど $\rho u_x^3$ 一つです。おおよそその程度という意味ではなく、余りなくその項一つです。

欠損項が運動量方程式に残す顔

欠けた3次モーメントは粘性応力にそのまま乗ります。Chapman–Enskog展開を最後まで進めると、応力に項が一つ加わります。

\sigma_{\alpha\beta}^{\text{err}} = -\left( \tau - \tfrac{1}{2} \right) \Delta t \, \partial_\gamma \left( \rho\, u_\alpha u_\beta u_\gamma \right)

ここで平均流 $U$ が $x$ 方向に流れ、その上に小さな擾乱 $u'$ が乗っているとします。上の項から $u'$ の1次だけを残すと $-(\tau - 1/2)\Delta t\, U^2 \partial_x^2 (\rho u'_\alpha)$ となり、これは粘性項と同じ形です。つまり $x$ 方向の粘性係数が丸ごと書き換わります。

\nu_{xx}^{\text{eff}} = \left( \tau - \tfrac{1}{2} \right) \Delta t \left( c_s^2 - U^2 \right)

相対誤差は $-U^2 / c_s^2$ です。ここで $\tau$ が分子と分母から同時に消えます。緩和をどれだけ調整しても割合は変わらないという意味です。そして誤差は速度の二乗に乗ります。マッハ数で書けば $-\mathrm{Ma}^2$ です。

Pythonで同じ渦を二つのフレームから測る#

Taylor–Green渦を64×64の周期格子に載せ、そこに一定速度 $U_0$ を加えます。振幅の減衰率を対数の傾きで測れば粘性係数が逆算できます。変えるのは $U_0$ だけで、残りはすべて固定します。

import numpy as np
 
CX  = np.array([0, 1, 0, -1,  0, 1, -1, -1,  1], dtype=float)
CY  = np.array([0, 0, 1,  0, -1, 1,  1, -1, -1], dtype=float)
W   = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
CS2 = 1.0 / 3.0
 
def f_equilibrium(rho, ux, uy):
    # 2次まで展開した標準的な平衡分布
    cu = CX[:, None, None] * ux + CY[:, None, None] * uy
    usq = ux * ux + uy * uy
    return W[:, None, None] * rho * (1 + cu / CS2 + cu * cu / (2 * CS2**2) - usq / (2 * CS2))
 
print("[1] third moment audit:  sum f_i c_ix^3   vs   Maxwell rho*(u^3 + 3*cs2*u)")
print("   u        lattice      Maxwell        gap        gap/u^3")
for u in (0.05, 0.10, 0.20, 0.30):
    feq = f_equilibrium(np.ones((1, 1)), np.full((1, 1), u), np.zeros((1, 1)))
    lat = float((feq * (CX**3)[:, None, None]).sum())
    exact = u**3 + 3 * CS2 * u
    print(f"  {u:4.2f}  {lat:11.8f}  {exact:11.8f}  {exact-lat:11.8f}   {(exact-lat)/u**3:8.5f}")
 
rng = np.random.default_rng(7)
frand = rng.random((9, 1, 1))
d = float((frand * (CX**3)[:, None, None]).sum() - (frand * CX[:, None, None]).sum())
print(f"  any random f:  sum f c_x^3 - sum f c_x = {d:.3e}   (c_x^3 = c_x, so always 0)")
 
def taylor_green_run(U0, tau=0.8, N=64, steps=900, amp=0.04, warmup=300, sample=25):
    # 同じ渦に一定速度U0を載せ、減衰率から粘性を測り直す
    k = 2 * np.pi / N
    x = np.arange(N)[:, None] * np.ones(N)[None, :]
    y = np.ones(N)[:, None] * np.arange(N)[None, :]
    ux = U0 - amp * np.cos(k * x) * np.sin(k * y)
    uy = amp * np.sin(k * x) * np.cos(k * y)
    rho = np.ones((N, N))
    f = f_equilibrium(rho, ux, uy)
    ts, logs = [], []
    for n in range(steps + 1):
        rho = f.sum(axis=0)
        ux = (f * CX[:, None, None]).sum(axis=0) / rho
        uy = (f * CY[:, None, None]).sum(axis=0) / rho
        if n % sample == 0:
            up, vp = ux - ux.mean(), uy - uy.mean()   # 平均流を引いた変動成分
            ts.append(n)
            logs.append(0.5 * np.log(2 * np.mean(up * up + vp * vp)))
        f += (f_equilibrium(rho, ux, uy) - f) / tau   # BGK衝突
        for i in range(9):                            # ストリーミング
            f[i] = np.roll(np.roll(f[i], int(CX[i]), axis=0), int(CY[i]), axis=1)
    ts, logs = np.array(ts, dtype=float), np.array(logs)
    m = ts >= warmup
    slope = np.polyfit(ts[m], logs[m], 1)[0]
    return -slope / (2 * k * k)
 
print("\n[2] same vortex measured in a uniformly moving frame  (tau=0.8, nu_theory=0.100000)")
nu_th = CS2 * (0.8 - 0.5)
print("   U0      nu_eff      rel.err     rel.err/U0^2")
for U0 in (0.00, 0.05, 0.10, 0.15, 0.20):
    nu = taylor_green_run(U0)
    e = nu / nu_th - 1
    tail = f"{e/U0**2:10.4f}" if U0 > 0 else "         -"
    print(f"  {U0:4.2f}  {nu:10.7f}  {e:+9.4%}  {tail}")
 
print("\n[3] does the error depend on tau?  (U0=0.15 fixed)")
print("   tau     nu_theory    nu_eff      rel.err")
for tau in (0.6, 0.8, 1.0):
    nu = taylor_green_run(0.15, tau=tau)
    th = CS2 * (tau - 0.5)
    print(f"  {tau:4.2f}  {th:10.7f}  {nu:10.7f}  {nu/th-1:+9.4%}")

[1] third moment audit:  sum f_i c_ix^3   vs   Maxwell rho*(u^3 + 3*cs2*u)
   u        lattice      Maxwell        gap        gap/u^3
  0.05   0.05000000   0.05012500   0.00012500    1.00000
  0.10   0.10000000   0.10100000   0.00100000    1.00000
  0.20   0.20000000   0.20800000   0.00800000    1.00000
  0.30   0.30000000   0.32700000   0.02700000    1.00000
  any random f:  sum f c_x^3 - sum f c_x = 0.000e+00   (c_x^3 = c_x, so always 0)
 
[2] same vortex measured in a uniformly moving frame  (tau=0.8, nu_theory=0.100000)
   U0      nu_eff      rel.err     rel.err/U0^2
  0.00   0.1000175   +0.0175%           -
  0.05   0.0996433   -0.3567%     -1.4268
  0.10   0.0985207   -1.4793%     -1.4793
  0.15   0.0966495   -3.3505%     -1.4891
  0.20   0.0940296   -5.9704%     -1.4926
 
[3] does the error depend on tau?  (U0=0.15 fixed)
   tau     nu_theory    nu_eff      rel.err
  0.60   0.0333333   0.0321975   -3.4076%
  0.80   0.1000000   0.0966495   -3.3505%
  1.00   0.1666667   0.1611502   -3.3099%

誤差は $\tau$ ではなく $U^2$ に乗っている#

[1]の最後の列がすべて1.00000です。欠損が $\rho u^3$ 一つだという主張が小数点以下五桁まで合っています。

[2]では静止フレームの誤差が $+0.0175\%$ です。測定方法そのもののノイズがその程度だという基準線になります。 $U_0$ を上げると誤差は $-0.36\%$ 、 $-1.48\%$ 、 $-3.35\%$ 、 $-5.97\%$ と大きくなります。最後の列は誤差を $U_0^2$ で割った値で、 $U_0$ が小さくなるほど $-1.5$ に寄ります。

この $-1.5$ が前節の予測です。相対誤差 $-U^2/c_s^2$ は $x$ 成分だけにかかりますが、Taylor–Greenモードは $k_x = k_y$ の対角モードなので、減衰率は二方向の平均になります。だから半分、つまり $-U^2/(2c_s^2) = -1.5 U^2$ です。 $U_0 = 0.2$ での予測は $-6.00\%$ 、測定値は $-5.97\%$ でした。

さらに重要なのが[3]です。 $\tau$ を0.6から1.0まで上げて粘性を五倍にしても、相対誤差は $-3.41\%$ から $-3.31\%$ までしか動きません。粘性を増やして誤差を埋めることはできないという意味です。

frame velocity U : 0.18 (Ma = 0.31)relaxation τ : 0.80

ν = 0.10000  |  νeff = 0.09514  |  error -4.86 %  |  step 0  |  amplitude ratio 1.000

frame velocity U を上げると、右の渦が流れていきながら左より遅く薄れます。同じ流体なのに右がなかなか消えません。 $\tau$ スライダーを揺らすと、下のerror値がほとんど動かないことも同時に見えます。

直す三つの方法と、それぞれの請求書

速度を下げる。 誤差が $U^2$ なので、格子速度を半分にすれば誤差は四分の一になります。代わりに同じ物理時間を再現するためのステップ数が増えます。LBMで「マッハ数を0.1以下に」という経験則は安定性の話として紹介されがちですが、実際に先に効くのはこの欠損項です。

補正項を加える。 $\partial_\gamma(\rho u_\alpha u_\beta u_\gamma)$ を差分で計算し、衝突段階に逆符号で入れます。追加コストは勾配一回分で、格子ボルツマンの長所である完全に局所的な衝突が少し崩れます。それでも回転座標系のように平均流が大きい問題では、いちばん安い選択肢です。

速度集合を広げる。 $\pm 2$ を含む多速度格子は $c^3 \neq c$ なので、3次モーメントを正確に合わせられます。最初のvizの speeds ±2 ボタンがまさにその計算です。対価はステンシルが広がること、メモリが増えること、境界処理が複雑になることです。一部の速度の重みが負になり得るので、安定性も見直しが必要になります。

この欠損が実際の計算で顔を出す場所

平均流が大きい場所がすべて候補です。回転機械の回転座標系、スライディングメッシュ、速い一様流に乗った乱流、動く物体に貼り付いたフレーム。静止流のベンチマークでは健全だったコードが、こうした問題でだけ有効レイノルズ数を外すなら、これが症状です。

LBM格子細分化の非平衡リスケールではレベルごとに $\tau$ を取り直す話をしましたが、この欠損はその操作では消えません。相対誤差が $\tau$ に依存しないからです。熱LBMに沈んでいる二つの定数で見た構造とも同じです。物性として入力していない値が、格子の構造から自動的に決まってしまいます。

診断は安く済みます。同じ問題に平均流を載せてもう一度回し、減衰率か抗力係数を比べるだけです。差が $U^2$ に比例して大きくなるなら、この欠損項です。STLボクセル化の偶奇判定のように境界側を疑う必要はありません。

次にフレームごとに粘性が変わって見えたら

格子ボルツマンで速度の上限を決めるのは安定性だけではありません。平衡分布が合わせられるモーメントの次数が先に上限を決めます。D2Q9は3次で $\rho u^3$ を落とし、その代価は $-\mathrm{Ma}^2$ の粘性誤差として請求されます。

数値がおかしいとき、格子を細かくするか $\tau$ をいじるかが先に思い浮かびます。この誤差はそのどちらでも減りません。減るのは $u/c_s$ を下げるか、欠損項を直接戻すか、速度集合を広げるときだけです。