速度を上げたら拡散が27%減った — LBM移流拡散モデルの余分なフラックス

格子ボルツマン法では拡散係数が緩和時間ひとつで決まります。$D = c_s^2(\tau - 1/2)\delta t$。短い式です。速度は入っていません。流れをオンにしても同じ $D$ が出るでしょうか。一様な移流速度 $u$ のもとでは $D(1 - u^2/c_s^2)$ まで落ちます。$u = 0.3$ なら27%が消えます。Chapman–Enskog 展開のどこでその項が落ちるのか、源項ひとつがどう戻すのかを、展開の順に追います。

痛いのは phase-field 多相流です。Cahn–Hilliard や Allen–Cahn を格子ボルツマンで解くと、界面厚さは移動度に縛られます。移動度が27%ずれると厚さがずれ、厚さがずれると表面張力係数がついてずれます。

拡散係数は合わせたのに界面が薄くなる

まず目で確かめましょう。以下は1次元の移流拡散方程式

を D1Q3 格子ボルツマンで解いたものです。初期条件はガウス分布ひとつ、厳密解もガウス分布で、幅は $\sigma^2(t) = \sigma_0^2 + 2Dt$ で育ちます。u を0から0.30へ押すと山が破線より細くなり、tau を回してもその比率は動きません。

u = 0.30 では実線（計算）が破線（厳密解）より鋭いです。下段の $\sigma^2$ は破線の傾き $2D$ に届きません。tau を0.6に置いても1.5に置いても、欠損の比率は同じマスに残ります。

D1Q3 が守るべきモーメントは3つだけ#

Navier–Stokes を解く LBM と違い、スカラー輸送用の分布関数 $g_i$ が満たすべきモーメントは少なくて済みます。Guo が2009年の非線形移流拡散方程式モデルで使った平衡分布はこれです。

$w_i$ は格子重み、$\mathbf{c}_i$ は格子速度、$c_s^2 = 1/3$ は格子音速の二乗です。この分布が満たすモーメントは3つあります。

流れ用の平衡分布と分かれるのが3番目です。ここには $\phi\mathbf{u}\mathbf{u}$ 項がありません。速度の2次項を最初から入れていません。スカラー方程式には応力テンソルがないので、2次モーメントに等方成分さえあれば足ります。格子も D2Q9 ではなく D2Q5 で済みます。拡散係数の前の $\tau - 1/2$ に入っている半ステップは、離散 LBE の2次展開が粘性の式に残すものと同じ項です。

十分に見えますし、$u = 0$ なら実際に完璧です。問題は、切り落とした分がただではないことです。

Chapman-Enskog が置いていく4番目の項#

BGK 格子ボルツマン方程式に源項 $S_i$ をつけて書きます。

$g_i = g_i^{(0)} + \varepsilon g_i^{(1)} + \varepsilon^2 g_i^{(2)}$ と展開し、時間微分を $\partial_t = \varepsilon\partial_{t_1} + \varepsilon^2\partial_{t_2}$ に分けます。$\varepsilon$ 次の0次モーメントは目標方程式の移流部分をそのまま返します。

同じ $\varepsilon$ 次の式の1次モーメントが核心です。ここで $g^{(1)}$ が運ぶフラックスが決まります。

括弧内の2番目の項 $c_s^2\nabla\phi$ が求めていた拡散フラックスです。ところが1番目の項が一緒に入っています。流れの LBM では平衡分布の2次モーメントに $\rho\mathbf{u}\mathbf{u}$ があってこの項がほぼ相殺されますが、スカラーモデルにはその項がありません。だから残ります。

$\varepsilon^2$ 次の0次モーメントまで合わせると最終形になります。

$D = c_s^2(\tau - 1/2)\delta t$ は予想どおりです。右辺の2番目の項が、誰も注文していない余分なフラックスです。定常状態で $u$ が時間的に一定、$\phi$ も変化しなければ消えます。しかし $\phi$ が動いている限り、つまり計算が進んでいる限り $\partial_t(\phi\mathbf{u})$ はゼロではありません。

一様な $u$ ではその項の正体は負の拡散です#

$u$ が空間・時間について定数という最も単純な場合を見ます。すると $\partial_t(\phi u) = u\,\partial_t\phi$ であり、先行次数で $\partial_t\phi \simeq -u\,\partial_x\phi$ です。代入すると

拡散項と符号が逆の拡散項です。物理拡散項と合わせると

欠損の形は二乗ひとつです。格子音速の二乗は3分の1なので、速度0.3の二乗0.09を3倍すると0.27になります。緩和時間が1なら、設定した拡散は3分の1掛ける2分の1、すなわち0.166667です。残る値は0.166667に0.73を掛けた0.121667です。差の0.045は、設定拡散に0.27を掛けたものと一致します。出力の $u = 0.30$ の行がこの引き算です。低速では二乗が小さくなります。速度0.05の二乗0.0025を3倍すると0.0075で、欠損は0.75%にとどまります。山は厳密解のすぐ隣に張り付きます。目では見逃しやすいです。

比率が緩和時間に縛られない理由も同じ行にあります。緩和時間から2分の1を引いた因子が、物理拡散と余分な項の両方に付きます。二つを割ると因子は消えます。だから比率は緩和時間と無関係です。緩和時間を0.6から1.5へ上げると、設定拡散は0.033333から0.333333へ10倍になります。速度0.25では0.0625を3倍した0.1875を1から引き、0.8125が残ります。その値が4行すべてに固定されます。絶対欠損だけが0.006付近から0.062付近へ一緒に育ちます。拡散を大きくして誤差を埋める場所はありません。

速度が格子音速、つまり3分の1の平方根の逆数である0.577付近に届くと、括弧の中が0になります。その先で有効拡散は負です。山は広がる代わりに自ら鋭くなり、分布は発散します。

また現れる $1 - 1/(2\tau)$#

源項 $S_i$ の係数を決めます。上の最終方程式で余分な項と源項が互いに消えるには

この条件を満たしつつ $\sum_i S_i = 0$ となる最も簡単な形はひとつだけです。

$1 - 1/(2\tau)$ がまた出てきました。この係数は力項が3次モーメントで越えてはいけない線と同じ場所にあります。離散時間の格子では源項が2回、一度は $g^{(1)}$ を通して、もう一度は Taylor 展開の2次項を通して効きます。係数 $1/2$ の差はそこで生まれます。

係数を落として $\sum_i \mathbf{c}_i S_i = \partial_t(\phi\mathbf{u})$ とするとどうなるでしょうか。$\tau = 1$ でちょうど2倍が入り、27%不足だった拡散が27%過剰になります。誤差の絶対値は変わらないので、対数スケールで収束次数を測っても気づきにくいのです。

$\partial_t(\phi\mathbf{u})$ はコード上、前の時間ステップの $\phi u$ を保存しておいて後退差分で求めれば済みます。配列がひとつ増えるだけがコストの全部です。

Python 60行で測った $D_{\text{eff}}$#

言葉で終わらせず測りましょう。ガウス分布を流しながら2次モーメント $\sigma^2$ の増加傾きを最小二乗で取り出せば、それが $2D_{\text{eff}}$ です。

import numpy as np
 
CS2 = 1.0 / 3.0
C = np.array([0, 1, -1])
W = np.array([2 / 3, 1 / 6, 1 / 6])
 
 
def d1q3_equilibrium(phi, u):
    """g_i^eq = w_i phi (1 + c_i u / cs^2) — モーメント3つだけ合わせた平衡分布"""
    return np.stack([W[i] * phi * (1.0 + C[i] * u / CS2) for i in range(3)])
 
 
def gaussian_moments(x, phi):
    m0 = phi.sum()
    mean = (x * phi).sum() / m0
    return mean, (((x - mean) ** 2) * phi).sum() / m0
 
 
def run_cde_lbm(L, steps, tau, u, sigma0, x0, corrected):
    x = np.arange(L, dtype=float)
    phi = np.exp(-((x - x0) ** 2) / (2 * sigma0**2))
    g = d1q3_equilibrium(phi, u)
    phi_old = phi.copy()
    hist = []
    for n in range(steps + 1):
        if n % 100 == 0:
            hist.append((n, gaussian_moments(x, phi)[1]))
        src = np.zeros_like(g)
        if corrected and n > 0:
            # S_i = w_i (1 - 1/(2 tau)) c_i d_t(phi u) / cs^2,  dt = 1
            dt_phiu = (1.0 - 1.0 / (2 * tau)) * u * (phi - phi_old)
            for i in range(3):
                src[i] = W[i] * C[i] * dt_phiu / CS2
        geq = d1q3_equilibrium(phi, u)
        g = g - (g - geq) / tau + src          # 衝突
        for i in range(3):
            g[i] = np.roll(g[i], C[i])         # 並進
        phi_old = phi
        phi = g.sum(axis=0)
    return np.array(hist)
 
 
def fit_diffusivity(hist):
    """sigma^2 = sigma0^2 + 2 D_eff t の傾きから D_eff を読む"""
    return np.polyfit(hist[:, 0], hist[:, 1], 1)[0] / 2.0
 
 
L, STEPS, SIG0, X0 = 800, 1600, 10.0, 80.0
 
tau = 1.0
D = CS2 * (tau - 0.5)
print(f"tau = {tau},  D = cs^2 (tau-1/2) = {D:.6f},  cs^2 = {CS2:.6f}")
print(f"{'u':>6} {'u^2/cs^2':>9} | {'D_eff (no src)':>14} {'ratio':>7} {'1-u^2/cs^2':>11} |"
      f" {'D_eff (src)':>12} {'ratio':>7}")
for u in [0.05, 0.10, 0.20, 0.30]:
    d_raw = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, False))
    d_fix = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, True))
    print(f"{u:>6.2f} {u * u / CS2:>9.4f} | {d_raw:>14.6f} {d_raw / D:>7.4f} {1 - u * u / CS2:>11.4f} |"
          f" {d_fix:>12.6f} {d_fix / D:>7.4f}")
 
u = 0.25
print(f"\nu = {u} fixed, tau sweep   (theory: ratio = 1 - u^2/cs^2 = {1 - u * u / CS2:.4f}, tau-independent)")
print(f"{'tau':>6} {'D':>10} | {'D_eff (no src)':>14} {'ratio':>7} | {'D_eff (src)':>12} {'ratio':>7}")
for tau in [0.6, 0.8, 1.0, 1.5]:
    D = CS2 * (tau - 0.5)
    d_raw = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, False))
    d_fix = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, True))
    print(f"{tau:>6.1f} {D:>10.6f} | {d_raw:>14.6f} {d_raw / D:>7.4f} | {d_fix:>12.6f} {d_fix / D:>7.4f}")

出力はこうなります。

tau = 1.0,  D = cs^2 (tau-1/2) = 0.166667,  cs^2 = 0.333333
     u  u^2/cs^2 | D_eff (no src)   ratio  1-u^2/cs^2 |  D_eff (src)   ratio
  0.05    0.0075 |       0.165417  0.9925      0.9925 |     0.166666  1.0000
  0.10    0.0300 |       0.161667  0.9700      0.9700 |     0.166666  1.0000
  0.20    0.1200 |       0.146667  0.8800      0.8800 |     0.166663  1.0000
  0.30    0.2700 |       0.121667  0.7300      0.7300 |     0.166658  0.9999
 
u = 0.25 fixed, tau sweep   (theory: ratio = 1 - u^2/cs^2 = 0.8125, tau-independent)
   tau          D | D_eff (no src)   ratio |  D_eff (src)   ratio
   0.6   0.033333 |       0.027096  0.8129 |     0.033345  1.0004
   0.8   0.100000 |       0.081258  0.8126 |     0.100006  1.0001
   1.0   0.166667 |       0.135417  0.8125 |     0.166661  1.0000
   1.5   0.333333 |       0.270794  0.8124 |     0.333275  0.9998

u を上げると赤い余分フラックスだけが $u^2$ で育ち、alpha を0から1へ引くと源項がその幅を消します。

alpha = 0 では赤いローブが青い物理フラックスを上下反転して敷かれます。alpha が1に届くと緑が赤の上に重なり、黄色い合計は青い線に戻ります。源項をオンにした列の ratio が1.0000に貼り付く場面がこれです。

1つ目の表の ratio 列と 1-u^2/cs^2 列が小数第4位まで一致しています。予測が近似ではなく厳密な先行次数だという意味です。源項をオンにすると4つの速度すべてが1.0000に戻ります。

2つ目の表のほうが手厳しいです。$\tau$ を0.6から1.5まで、$D$ を10倍動かしても欠損比率は0.8129から0.8124までしか動きません。拡散を大きく取って誤差を埋めようとする試みが通用しないということです。$D$ を10倍にすれば消える量も10倍になります。

格子速度の予算はすでに使われていました

$u^2/c_s^2$ は格子マッハ数の二乗です。$u = 0.1$ なら $0.01 \times 3 = 0.03$ で、拡散の3%が消えます。$u = 0.2$ なら $0.04 \times 3 = 0.12$、12%です。0.05以下では0.75%で、格子誤差に埋もれます。圧縮性誤差のために $u$ を0.1付近へ縛る慣行は、スカラー側の請求書をまだ見ていません。界面厚さや Sherwood 数を数字で出すなら、$u$ が0.1を超えた瞬間に源項をオンにするのが合います。その下では配列をひとつ省いても、表の0.9925は残ります。

phase-field では $u^2/c_s^2$ だけを信じてはいけません。余分な項は $\partial_t(\phi\mathbf{u})$ 全体で、二乗の形は $u$ が一様で定常なときだけ成り立ちます。界面付近で $u$ が跳ね、$\partial_t u$ もゼロでなければ $\phi\,\partial_t u$ まで生き残ります。$\partial_t(\phi\mathbf{u})$ は前ステップの $\phi u$ を後退差分で入れます。$1 - 1/(2\tau)$ の $\tau$ は1次モーメントの緩和時間です。速度フレームをずらしたとき粘性がどう減るかは、平衡分布の3次モーメント欠損と同じ種類の請求書です。移動度の算数が合っているのに界面だけ薄くなるなら、$\tau$ が違うのではありません。流れが拡散の予算を先に持っていったのです。