存在しない流体で界面を埋める — Ghost Fluid Method

水と空気が出会うセルに、本当はそこに無い偽の空気を流し込みます。このゴースト流体(界面を越えて外挿した仮想流体)は格子のどこにも存在しません。それでもそのゴーストこそが、衝撃波が界面を横切るときに計算を安定させます。コピーする量と外挿する量は別の仕事です。オリジナルGFMが震える場所は、その分け方を誤ったところです。

密度ジャンプを上げると階段は急になり、圧力は赤い破線をそのまま通ります。「show ghost」は界面の向こうの延長線を付けたり消したりします。

破線を付けると、各流体が自分のエントロピーを持ったまま反対側へ伸びます。界面の位置を動かしても圧力曲線は一本のままです。切れるのは密度です。

一つのセルで二つの物質がぶつかるとき

気体と液体、火薬と金属のように、物性が大きく異なる二つの物質が同じ格子の上で出会います。単純に平均物性を使うと界面が数値的ににじみます。密度比1000倍の水-空気界面では、このにじみがすぐに非物理的な振動へと広がります。

GFMの発想は違います。界面をぼかして混ぜません。各流体を自分の領域の外へ偽として延ばし、片方の物質しか無いかのように単一物質ソルバを二度走らせます。界面はレベルセット(符号で領域を分ける距離関数)$\phi$ の零点として追跡します。

ここで $\phi$ は符号付き距離関数、その零点 $\phi=0$ が界面です。

界面を横切って切れるものと続くもの

ゴーストを正しく埋めるには、界面を横切って何が連続かを正確に知る必要があります。圧力 $p$ と法線速度 $u_n$ は界面で連続です。一方、エントロピー $s$ と接線速度 $u_t$ は不連続です。

角括弧 $[\cdot]$ は界面を横切る飛びを表します。この四行がGFMのすべてです。連続量はコピーし、不連続量は外挿します。上の図で圧力だけが滑らかで密度だけが階段なのは、この分け方の結果です。

一つのセルの中で二相を平均すると、別の故障が起きます。6方程式 two-fluid モデルは密度比が大きい界面で固有値が複素になり、格子を細かくすると発散が早く来ます。GFMはその平均を取りません。セルはいつも片方の物質だけを持ちます。

水を気体のように — Stiffened EOS#

水はほぼ非圧縮なので、理想気体の状態方程式では扱えません。GFMは二つの物質を同じ形式のMie–Grüneisen系にまとめるためにStiffened Gas EOSを使います。

$\rho$ は密度、$e$ は単位質量あたりの内部エネルギー、$\gamma$ は比熱比、$p_\infty$ は物質剛性を表す定数です。理想気体は $p_\infty = 0$、水はおよそ $\gamma \approx 4.4$、$p_\infty \approx 6 \times 10^8\,\mathrm{Pa}$ に合わせます。音速は次のようになります。

$p_\infty$ のおかげで、同じ圧力でも水の音速は空気よりずっと大きく出ます。EOSを一行に統一すれば、二つの流体を同じソルバに入れられます。

ゴーストセルを埋める二つの手つき

ここが核心の手順です。流体Aの領域を解くとき、Bが占めるセルをAのゴーストで埋めます。手つきは二つに分かれます。

1. コピー：連続量である圧力と法線速度は、界面の向こうの実際のBセルの値をそのまま取ります。
2. 外挿：不連続量であるエントロピー(つまり密度)と接線速度は、最も近い実際のAセルから界面を越えて等エントロピーで外挿します。

等エントロピー外挿はエントロピー $s = p / \rho^\gamma$ を保ったまま密度を解き直します。

$s_A$ は実際のA流体のエントロピー、$p_\text{copy}$ はコピーしてきた圧力です。外挿は通常、界面の法線方向へ速く広げるfast marching法で行います。一度のスイープでAのゴースト場が完成したら、Bの存在を忘れてAだけの単一物質Euler方程式を解きます。Bについても同じ手順を繰り返します。

コードに書いてある初期値だけでゴースト密度を計算すると、外挿が実密度をどうずらすかが見えます。左側は密度1、圧力1、比熱比1.4です。エントロピー $s = p / \rho^\gamma$ は1です。右側の圧力0.1をコピーして左側の流体を延ばすと、ゴースト密度は $0.1^{1/1.4}$ で約0.193になります。実際の右側密度0.125ではありません。単一物質ソルバが見る密度比は8倍ではなく約5.2倍です。

逆向きも同じ式です。右側の比熱比は1.67です。エントロピーは0.1を $0.125^{1.67}$ で割った約3.22です。左側の圧力1をこのエントロピーに載せると、ゴースト密度は約0.496です。実際の左側密度1の半分近くです。圧力は渡し、各流体のエントロピーはそのゴーストに残します。

音速も状態方程式からそのまま出ます。剛性定数が0なら、音速は比熱比と圧力と密度から作る平方根です。左側は約1.183、右側は約1.156です。点400個で長さ1なら間隔は $1/399$、約0.00251です。CFL 0.4を左側の音速に掛けると時間刻みは約0.00085です。終了時刻0.14を割ると165ステップ前後です。音波が長さ1を渡る時間は約0.85なので、0.14はその5分の1に届きません。水-空気の密度比1000に比べれば、この衝撃管の8倍は穏やかです。それでもゴースト密度は0.193と0.496で、すでに実値からずれています。

オリジナルGFMが振動する場所#

オリジナルGFMは、界面の片側が非常に剛い(stiff、水のように $p_\infty$ が大きい)とき問題が起きます。強い衝撃波やデトネーション波が界面を横切ると、単純なコピー・外挿で作ったゴースト状態が本当のRiemann解とずれます。そのずれが界面付近で圧力振動として現れ、振動はしばしば計算の発散につながります。

下のシミュレーションで圧力パルスを剛い流体(ピンク領域)に通してみましょう。stiffness ratioを上げるほど、界面を通過した直後のさざ波が大きくなります。

stiffness ratioを9まで上げると、界面直後に高周波振動がはっきり見えます。この振動は格子を細かくしても消えません。原因が格子ではなく境界条件だからです。

修正GFM — 界面Riemann問題で値を決める#

解法は、ゴースト状態を推測せず解いて求めることです。修正GFM(Modified GFM、MGFMまたはIGFM)は界面両側の実際の状態を左右として、Riemann問題を解きます。

$W_A, W_B$ は界面両側の実際の状態、$\mathcal{R}$ はRiemann解法、$p^*, u^*$ は接触不連続の圧力・速度です。この $p^*, u^*$ を両側のゴーストに同じく与えると、圧力・速度の連続条件が正確に満たされます。上のシミュレーションで「use modified GFM」をオンにすると、同じ剛性でも振動が消えます。

Python — 気体-気体GFM衝撃管#

二つの理想気体($\gamma_A = 1.4$、$\gamma_B = 1.67$)が出会う衝撃管を、等エントロピー外挿GFMで解きます。HLLフラックスを使い、各ステップで二度のゴーストスイープの後にレベルセットで併合します。

import numpy as np
 
GAMMA = (1.4, 1.67)   # 左(空気)・右(単原子気体)の比熱比
 
def primitive(U, g):
    rho = U[0]; u = U[1] / rho
    e = U[2] / rho - 0.5 * u * u
    p = (g - 1.0) * rho * e          # Stiffened EOS, p_inf = 0 (理想気体)
    return rho, u, p
 
def conserved(rho, u, p, g):
    E = p / (g - 1.0) + 0.5 * rho * u * u
    return np.array([rho, rho * u, E])
 
def hll(UL, UR, g):
    rL, uL, pL = primitive(UL, g); rR, uR, pR = primitive(UR, g)
    aL = np.sqrt(g * pL / rL); aR = np.sqrt(g * pR / rR)
    sL = min(uL - aL, uR - aR); sR = max(uL + aL, uR + aR)
    FL = np.array([rL * uL, rL * uL * uL + pL, uL * (UL[2] + pL)])
    FR = np.array([rR * uR, rR * uR * uR + pR, uR * (UR[2] + pR)])
    if sL >= 0: return FL
    if sR <= 0: return FR
    return (sR * FL - sL * FR + sL * sR * (UR - UL)) / (sR - sL)
 
def extrapolate_ghost(rho, u, p, mat, side):
    # side 流体で全領域を埋める: 圧力・速度はコピー, 密度は等エントロピー外挿
    rg, ug, pg = rho.copy(), u.copy(), p.copy()
    g = GAMMA[side]
    idx = np.where(mat == side)[0]
    lo, hi = idx[0], idx[-1]
    ref = lo if side == 1 else hi      # 界面に最も近い実セル
    s_ref = p[ref] / rho[ref] ** g     # 保存するエントロピー
    for i in range(len(rho)):
        if mat[i] != side:
            j = min(max(i, lo), hi)            # 最も近い実セル
            ug[i] = u[j]; pg[i] = p[j]         # コピー
            rg[i] = (pg[i] / s_ref) ** (1.0 / g)   # 等エントロピー外挿
    return rg, ug, pg
 
def run_ghost_fluid_tube(N=400, t_end=0.14, cfl=0.4):
    x = np.linspace(0, 1, N); dx = x[1] - x[0]
    x0 = 0.5                                   # 初期界面位置
    mat = np.where(x < x0, 0, 1)
    rho = np.where(x < x0, 1.0, 0.125)
    u   = np.zeros(N)
    p   = np.where(x < x0, 1.0, 0.1)
    t = 0.0
    while t < t_end:
        a = np.sqrt(np.where(mat == 0, GAMMA[0], GAMMA[1]) * p / rho)
        dt = min(cfl * dx / np.max(np.abs(u) + a), t_end - t)
        new = {}
        for side in (0, 1):                    # 二度のゴーストスイープ
            rg, ug, pg = extrapolate_ghost(rho, u, p, mat, side)
            g = GAMMA[side]
            U = np.array([conserved(rg[i], ug[i], pg[i], g) for i in range(N)]).T
            F = np.zeros((3, N + 1))
            for i in range(1, N):
                F[:, i] = hll(U[:, i - 1], U[:, i], g)
            F[:, 0] = F[:, 1]; F[:, N] = F[:, N - 1]
            new[side] = U - dt / dx * (F[:, 1:] - F[:, :-1])
        U = np.where(mat[None, :] == 0, new[0], new[1])   # レベルセット併合
        for i in range(N):
            rho[i], u[i], p[i] = primitive(U[:, i], GAMMA[mat[i]])
        x0 += dt * u[np.argmin(np.abs(x - x0))]           # 界面の移動
        mat = np.where(x < x0, 0, 1)
        t += dt
    return x, rho, u, p
 
if __name__ == "__main__":
    x, rho, u, p = run_ghost_fluid_tube()
    print(f"t=0.14:  max p = {p.max():.3f}, "
          f"接触面付近の密度 {rho[180]:.3f} -> {rho[220]:.3f}")

出力は、界面を横切る密度の飛びが生き残り、圧力が滑らかにつながったことを示します。二度の単一物質スイープがGFMの背骨です。

コードを書くとき外してはいけない落とし穴

外挿を格子軸に流すと、斜めの界面でエントロピーが混ざります。広げる方向は界面法線だけです。剛性比が大きい水-空気では、オリジナルのコピーと外挿が圧力のさざ波をほぼ必ず作ります。その場合は最初から界面Riemann解でゴーストを埋めます。レベルセットを符号付き距離へ戻さないと、零点は残るのに距離の値が崩れ、界面速度がずれます。

気体同士で通った外挿を水へそのまま置くとき

$p_\infty$ が0の理想気体二つ($\gamma = 1.4$ と $1.67$)は、等エントロピー外挿だけでも接触面の密度段差が残ります。水は $\gamma \approx 4.4$、$p_\infty \approx 6 \times 10^8\,\mathrm{Pa}$ です。同じ圧力で音速の分子にこの定数が乗ると、デモで stiffness ratio を9まで上げて見たさざ波は、もっと大きいずれの縮小版になります。格子を細かくしても、そのさざ波は境界条件の中に残ります。

ゴーストに入れる圧力と速度をコピーのままにするか、界面Riemann解の $p^*, u^*$ に変えるかは、密度ジャンプを数える前に決めます。外挿密度が負へ崩れる場所は、リミッタの問題が別です。正値保存 flux リミッタが、質量は合うのに密度が負になるセルをどう止めるか。剛い界面の次に見る価値があります。