τ=0.51で粘性を51倍も大きく読んでいた — 離散LBEの2次展開が残す半ステップ

連続式が与える粘性と格子が返す粘性が51倍違った#

連続BGK方程式だけを見れば、緩和時間 $\tau$ と動粘性係数 $\nu$ の関係は $\nu = c_s^2 \tau$ です。$\tau$ を0.51にすれば $\nu$ は0.17になるはずです。ところが同じ $\tau$ をD2Q9コードに入れ、Taylor–Green渦の減衰率を測ると0.00333が出ます。51倍の違いです。

これはバグではありません。符号の取り違えでもありません。緩和時間からちょうど $\Delta t/2$ が抜けた値です。半ステップは二通りに出ます。特性曲線上の2次Taylor展開と、衝突積分の中点則です。二つの道は同じ場所で合流します。$\tau$ のスライダーを0.51側へ寄せると、橙の破線 $c_s^2\tau$ が格子の振幅から離れます。

青い曲線は32×32格子の上で生きている渦の振幅です。緑の破線は $c_s^2(\tau - 1/2)$ です。$\tau = 0.51$ で緑は0.00333、橙は0.17に分かれます。

格子上の一行は微分方程式ではない

実際に回しているのはこの一行です。

$f_i$ は離散速度 $\mathbf{c}_i$ 方向の分布関数、$\Omega_i = f_i - f_i^{\rm eq}$ は非平衡部分（平衡からのずれ）、$\Delta t$ は1ステップです。

左辺は微分ではありません。有限の間隔 $\Delta t$ を実際に飛び越えた差です。これを微分に直すには展開が必要で、展開には2次項があります。CFDでよくあるのは1次項だけ残して残りを「打ち切り誤差」と呼ぶことです。ここではその2次項が誤差ではなく物理係数を書き換えます。

無視できない理由は次数を数えれば見えます。2次項の大きさは $\Delta t^2 D^2 f / 2$ です。右辺の衝突項は $\Delta t / \tau$ で割られています。格子単位では $\tau$ は0.5から2の間の $O(1)$ の数で、$\Delta t = 1$ です。つまり2次項と1次項の比は小さくなりません。無限小になるはずの $\Delta t$ を1に固定したことが格子気体法の出発点であり、その代償がこの項です。

特性曲線に沿って2次まで展開すると何が残るか#

$D_i \equiv \partial_t + \mathbf{c}_i \cdot \nabla$ と置くと左辺は $D_i$ の指数関数になります。2次まで残して $\Delta t$ で割ります。

問題の項が $D_i^2$ です。これを捨てると $\nu = c_s^2\tau$ になり、残すと別の値になります。

Chapman–Enskog展開をかけてみましょう。$f_i = f_i^{(0)} + \varepsilon f_i^{(1)} + \varepsilon^2 f_i^{(2)}$、$\partial_t = \varepsilon \partial_{t_1} + \varepsilon^2 \partial_{t_2}$、$\nabla = \varepsilon \nabla_1$、そして $f_i^{(0)} = f_i^{\rm eq}$ です。$\varepsilon$ の次数では

が出ます。この式を $\varepsilon^2$ の次数でそのまま再利用するのが要点です。$D_{i1}^2 f_i^{(0)} = -\tfrac{1}{\tau} D_{i1} f_i^{(1)}$ なので

となります。$D_{i1} f_i^{(1)}$ の前に括弧が一つ付きました。粘性応力は $f_i^{(1)}$ の2次モーメントに $\tau$ を掛けて出るので、その $\tau$ が $\tau\left(1 - \frac{\Delta t}{2\tau}\right) = \tau - \frac{\Delta t}{2}$ に変わります。

格子単位 $\Delta t = 1$ では $c_s^2(\tau - 1/2)$ です。$c_s^2 = 1/3$ なので $\tau = 0.51$ なら $\nu = 0.00333$ です。上のシミュレーションの緑の破線がまさにこの値です。

残った項を分布関数の中へ押し込む

展開を使わずに同じ結果へ至る道もあります。連続BGKを特性曲線上で積分し、台形則を使うのです。緩和時間を $\lambda$ と書きます。

右辺に未来時刻が入っているので陰解法です。このままでは回せません。そこで変数を変えます。

$\bar f_i$ で書き直すと未来時刻の項が打ち消えて消えます。$\bar\Omega_i = \bar f_i - f_i^{\rm eq} = \left(1 + \frac{\Delta t}{2\lambda}\right)\Omega_i$ を代入すると残るのはこれです。

完全に陽解法で、私たちが使っているあの一行と文字まで同じです。ただし $\tau = \lambda + \Delta t/2$ のときだけです。つまりコードの中の $\tau$ は物理的な緩和時間ではありません。物理的なほうは $\lambda = \tau - \Delta t/2$ です。丸めの問題ではありません。変数変換が請求した値です。台形積分がこの半ステップを残す場所は緩和時間に付いた $\Delta t/2$に続きます。

同じ一行を中点則として読み直す

台形則は $\theta$ 重み則の $\theta = 1/2$ の場合です。$\theta$ を動かせば、なぜちょうど半分なのかが見えます。1ステップの増幅率は $A_\theta(z) = \frac{1 - (1-\theta)z}{1 + \theta z}$ で、正解は $e^{-z}$ です。ここで $z = \Delta t/\lambda$ です。

$\theta$ を1/2から外すと、1ステップ増幅率の誤差は $(\theta - \tfrac12)z^2$ として傾き2に乗ります。$\theta = 1/2$ でだけ傾き3に落ちます。半分は調整値ではなく、$e^{-z}$ の2次係数 $1/2$ と噛み合う唯一の値です。

Taylor–Green渦をτ6通りで減衰させた#

主張はここまでです。あとは格子に直接聞きます。2次元Taylor–Green渦は $\nu$ を測る定規になります。振幅が $e^{-2\nu k^2 t}$ で落ちるからです。

import numpy as np
 
C = np.array([[0,0],[1,0],[0,1],[-1,0],[0,-1],[1,1],[-1,1],[-1,-1],[1,-1]])
W = np.array([4/9] + [1/9]*4 + [1/36]*4)
CS2 = 1/3
 
def equilibrium_d2q9(rho, ux, uy):
    cu = C[:,0,None,None]*ux + C[:,1,None,None]*uy
    u2 = ux*ux + uy*uy
    return W[:,None,None]*rho*(1 + cu/CS2 + cu*cu/(2*CS2**2) - u2/(2*CS2))
 
def taylor_green_init(n, u0):
    k = 2*np.pi/n
    x, y = np.meshgrid(np.arange(n), np.arange(n), indexing='ij')
    ux = -u0*np.cos(k*x)*np.sin(k*y)
    uy =  u0*np.sin(k*x)*np.cos(k*y)
    return np.ones((n,n)), ux, uy, k
 
def advance_lattice(f, tau):
    rho = f.sum(axis=0)
    ux = (C[:,0,None,None]*f).sum(axis=0)/rho
    uy = (C[:,1,None,None]*f).sum(axis=0)/rho
    f = f - (f - equilibrium_d2q9(rho, ux, uy))/tau      # 衝突
    for i in range(9):                                    # 移流
        f[i] = np.roll(f[i], (C[i,0], C[i,1]), axis=(0,1))
    return f, ux, uy
 
def measure_decay_rate(n, tau, u0=0.01, nstep=400, skip=100):
    rho, ux, uy, k = taylor_green_init(n, u0)
    f = equilibrium_d2q9(rho, ux, uy)
    amp = []
    for s in range(nstep):
        f, ux, uy = advance_lattice(f, tau)
        amp.append(np.sqrt((ux**2 + uy**2).mean()))
    t = np.arange(skip, nstep, dtype=float)
    a = np.log(np.array(amp[skip:]))                      # 初期の過渡区間は捨てる
    slope = np.polyfit(t, a, 1)[0]
    return -slope/(2*k*k), k
 
def viscosity_scan(n=32):
    print(f"{'tau':>6} {'nu_measured':>12} {'cs2*(tau-1/2)':>14} {'cs2*tau':>10} "
          f"{'err_half[%]':>11} {'err_naive[%]':>12}")
    for tau in (0.51, 0.55, 0.6, 0.8, 1.0, 1.5):
        nu, k = measure_decay_rate(n, tau)
        half = CS2*(tau - 0.5)
        naive = CS2*tau
        print(f"{tau:6.2f} {nu:12.6f} {half:14.6f} {naive:10.6f} "
              f"{100*abs(nu-half)/half:11.2f} {100*abs(nu-naive)/naive:12.2f}")
 
viscosity_scan()

   tau  nu_measured  cs2*(tau-1/2)    cs2*tau err_half[%] err_naive[%]
  0.51     0.003333       0.003333   0.170000        0.01        98.04
  0.55     0.016667       0.016667   0.183333        0.00        90.91
  0.60     0.033337       0.033333   0.200000        0.01        83.33
  0.80     0.100115       0.100000   0.266667        0.12        62.46
  1.00     0.167205       0.166667   0.333333        0.32        49.84
  1.50     0.337688       0.333333   0.500000        1.31        32.46

theta を1/2から外すと右の傾きは2に付き、1/2でだけ3に落ちます。

左の青い領域が $e^{-s/\lambda}$ の面積で、桃色の破線が $\theta$ 則の代わりに使う高さです。$\theta = 1/2$ だけが $e^{-z}$ の2次係数 $1/2$ と噛み合います。表の半ステップ誤差0.01%は、その中点が合っているときの残りです。

$\tau - 1/2$ 側の誤差は0.01〜1.31%です。格子の離散化が残す高次誤差の水準です。$\tau$ 側は32〜98%です。誤差と呼べる大きさではありません。

初期振幅を0.01と低くしたのは圧縮性の影響を消すためです。マッハ数が0.02程度なら、$f^{\rm eq}$ の打ち切られた3次項が減衰率に混ざりません。最初の100ステップを捨てるのも同じ理由です。初期条件を平衡分布で与えると $f^{(1)}$ がまだ0なので、非平衡成分が落ち着くまで数十ステップかかります。その区間を入れて傾きを測ると粘性が実際より小さく出ます。

$\tau$ が大きくなるほど $\tau - 1/2$ 側の誤差が少しずつ育つのも読みどころです。$\tau = 1.5$ で1.31%になります。展開で捨てた $O(\Delta t^2)$ の項が $\tau$ とともに大きくなるからです。半ステップの補正は完結ではなく、2次までの話です。

表を $\tau/(\tau - 1/2)$ で読み直すと、51倍は一行です。$\tau = 0.51$ なら $0.51/0.01 = 51$ です。$0.170000/0.003333$ も51.0です。$\tau = 0.55$ は $0.55/0.05 = 11$ で、$0.183333/0.016667$ も11です。$\tau = 0.60$ は6、$\tau = 1.0$ は2、$\tau = 1.5$ は1.5です。ナイーブ誤差の98.04%、90.91%、83.33%、62.46%、49.84%、32.46%は、この比が1へ近づきながら縮んだ跡です。

半ステップ側はもっと静かです。$\tau = 0.51$ では測定も予測も0.003333で、誤差は0.01%です。$\tau = 0.80$ は測定0.100115、予測0.100000、差0.000115で0.12%です。$\tau = 1.50$ は0.337688対0.333333、差0.004355で1.31%です。捨てた $O(\Delta t^2)$ が $\tau$ とともに大きくなる、という話がこの三つのマスに入っています。初期振幅0.01、先頭100ステップ、格子32は、その残りを圧縮性から切り離すために選んだ値です。マッハ0.02なら、切られた3次項は減衰率に混ざりません。$\tau = 0.55$ のナイーブ誤差90.91%は $1 - 0.016667/0.183333$ と同じです。その比は1/11なので、連続式 $\nu = c_s^2\tau$ で読むと十一倍の粘性をコードに書くことになります。$\tau = 1.5$ でも比は1.5で、ナイーブ誤差32.46%を小さい誤差とは呼べません。

境界条件と外力項が同じ半ステップを再度請求する

この半ステップは一度しか現れないわけではありません。

bounce-back境界が壁をノード上ではなく二つのノードの中間に置く理由も同じです。分布関数が反射する時点がステップの中間だからです。$\tau$ が0.5から離れるとこの実効的な壁位置が $\tau$ に依存して動き、それが有名な2次精度の崩れです。

Guo形式の外力項に $\left(1 - \frac{1}{2\tau}\right)$ が付く理由も同じです。外力も同じ特性曲線上で積分され、同じ台形則を受けます。$\bar f$ 変換が $\Omega_i$ にしか作用しないので係数が残ります。この括弧を落とすと、静止しているべき流体がゆっくり流れたり、重力下の静水圧分布がわずかに傾いたりします。

格子幅を半分にするとき $\tau$ を単純に2倍してはいけない理由もここにあります。$\nu$ を保つには $\tau_{\rm fine} - 1/2 = 2(\tau_{\rm coarse} - 1/2)$ を解く必要があります。細分化ループの中でこの一行をどこに置くかが、誤差を分けます。

$\tau$ をノードごとに変える場合、下限はまさにこの半ステップです。$\tau \to 0.5$ で $\nu \to 0$ なので、粘性が局所的に0へ近づく問題は格子そのものの壁に当たります。

0.5は下限ではなく座標原点である#

$\tau > 0.5$ を安定条件として暗記すると、原因を落とします。物理的な緩和時間は $\tau - 1/2$ です。$\tau = 0.51$ では $c_s^2\tau = 0.170000$、$c_s^2(\tau - 1/2) = 0.003333$ なので比は51です。負の緩和時間は物理ではありません。

コードに書くべきは $\tau$ ではなく $\nu$ です。$\tau = 1/2 + \nu/c_s^2$ を入力の一箇所に置きます。外力項に付く $\left(1 - \frac{1}{2\tau}\right)$ も、同じ半ステップの請求書です。その括弧を力の3次モーメントと並べて置く限界は力項が守るべき線に続きます。減衰が51倍遅くなれば、どこかで $\tau$ を $\nu$ として読んでいるということです。