Kを消したら水が1.86倍圧縮された — 格子を細かくしても減らない界面の誤差

界面セルではψ+Kがほぼゼロになります。Kはψ∇·uを消すために入る係数で、その項を外すと液体が気体の圧縮率で押されます。

圧縮性VOFの色関数の移流式には項がもう一つ付きます。 $-(\psi+K)\nabla\cdot u$ です。その $K$ をコードに移して値を出力してみました。水と空気の界面セルで $\psi=0.5$ のとき $K=-0.499947$ 、したがって $\psi+K=5.3\times10^{-5}$ でした。式を正しくするために入れた係数が、その項をほぼゼロに消してしまったわけです。

本記事ではDennerの圧力ベース完全連成圧縮性多相有限体積法の論文から、 $K$ を等圧クロージャで直接導きます。そのうえで $K$ をゼロにした場合と符号を反転した場合に、1次元圧縮問題で液体体積がそれぞれ1.86倍・4.23倍多く減ることを再現します。最後に、その誤差が格子を8倍細かくしても7%しか減らない理由を見ます。

圧縮性VOFに $-(\psi+K)\nabla\cdot u$ が付いている理由#

非圧縮性VOFでは色関数はただ流れていきます。 $\nabla\cdot u=0$ なので移流だけで済みます。圧縮性になると流体が膨張し収縮します。色関数は体積率なので、その膨張・収縮に追従しなければなりません。

論文の式はこうです。

\frac{\partial \psi}{\partial t} + \frac{\partial u_i \psi}{\partial x_i} - (\psi + K)\frac{\partial u_i}{\partial x_i} = 0

$\psi$ は液体の体積率、 $u_i$ は速度、 $K$ は物質によって変わる圧縮率因子(compressibility factor)です。第三の項を移項すると物質微分の形になります。

\frac{\mathrm{D}\psi}{\mathrm{D}t} = K\,\nabla\cdot u

$K=0$ なら体積率は物質に沿って変化しません。もっともらしく見えますが誤りです。界面セルが見る $\nabla\cdot u$ は混合物全体の値であり、混合物の圧縮率は柔らかい側が支配します。

数字で見ましょう。1気圧の水の体積弾性率は $\rho_l a_l^2 = 2.64\times10^9$ Pa、空気は $\rho_g a_g^2 = 1.40\times10^5$ Pa です。比は18,857倍。混合物の圧縮率は体積率で重み付けした和なので、 $\psi=0.5$ でも空気の項がほぼ全部を占めます。その結果、混合物の音速は23.65 m/sまで落ちます。水より69倍、空気より14倍遅い値です。

下のシミュレーションで体積率を直接動かしてみてください。

ψ (liquid volume fraction) = 0.500pressure = 1 bar

a_mix = 23.65 m/s  |  ρ_mix = 500.60 kg/m³  |  a_mix / a_gas = 0.069
crossing 2.0 m: liquid 1.23 ms  |  gas 5.85 ms  |  mixture 84.57 ms

$\psi$ スライダーをゼロの少し上まで上げるだけで、音速が341 m/sから100 m/s以下へ垂直に落ちます。下のトラックは、同じ2 mの区間を音響パルスが渡るのにかかる時間です。液体が体積比で0.1%混ざるだけで、混合物のレーンが気体のレーンから目に見えて遅れます。

$K$ はどこから出てくる数字か#

$K$ は仮定から出てきます。等圧クロージャ(isobaric closure) — 界面セル内の2相が同じ圧力を共有し、それぞれ自分の体積弾性率で応答するという仮定です。論文がACIDで使う局所単相仮定と同じ内容です。

ラグランジュ的なセルを追いながら、相 $k$ の質量が保存すると書くと

\frac{1}{V_k}\frac{\mathrm{D}V_k}{\mathrm{D}t} = -\frac{1}{\rho_k a_k^2}\frac{\mathrm{D}p}{\mathrm{D}t}

になります。 $V_k$ はその相が占める体積、 $\rho_k a_k^2$ は体積弾性率です。両相を足すと混合物の発散が出ます。

\nabla\cdot u = -\left[\frac{\psi}{\rho_l a_l^2} + \frac{1-\psi}{\rho_g a_g^2}\right]\frac{\mathrm{D}p}{\mathrm{D}t}

角括弧の中が混合物の圧縮率です。ここからWood音速が自動的に出てきます。次に $\psi = V_l/V$ を微分して上の2式を代入すると

\frac{\mathrm{D}\psi}{\mathrm{D}t} = \psi(1-\psi)\left[\frac{1}{\rho_g a_g^2} - \frac{1}{\rho_l a_l^2}\right]\frac{\mathrm{D}p}{\mathrm{D}t}

が残ります。掛かっているのは2つの圧縮率の差です。両相の体積弾性率が同じなら体積率は動きません。 $\mathrm{D}\psi/\mathrm{D}t = K\nabla\cdot u$ と比べて $\mathrm{D}p/\mathrm{D}t$ を消去すると $K$ が落ちてきます。

K = \frac{\rho_g a_g^2 - \rho_l a_l^2}{\dfrac{\rho_g a_g^2}{1-\psi} + \dfrac{\rho_l a_l^2}{\psi}}

分母に $1/\psi$ と $1/(1-\psi)$ があります。純粋セルで発散するので、コードには約分した形で書くべきです。

K = \frac{(\rho_g a_g^2 - \rho_l a_l^2)\,\psi(1-\psi)}{\rho_g a_g^2\,\psi + \rho_l a_l^2\,(1-\psi)}

この形は $\psi=0$ と $\psi=1$ でそのままゼロになります。分子の符号は $\psi$ がどちらの相を指すかに完全に依存します。論文の式をそのまま写すとラベルが入れ替わって符号が反転しかねません。その代償は後で測ります。

$\psi+K$ が $5\times10^{-5}$ なら両方の項を消してよいか#

$\rho_l a_l^2 \gg \rho_g a_g^2$ の極限を入れてみます。分子は $-\rho_l a_l^2$ 、分母は $\rho_l a_l^2(1-\psi)$ に収束するので

K \to -\psi, \qquad \psi + K \to 0

です。物理的には当然です。水はほとんど圧縮されないので、水の体積は移流されるだけであるべきです。 $(\psi+K)\nabla\cdot u$ をゼロにすることが、まさにその挙動にあたります。

では両方の項を使わなければよいのでは。純粋な液体セルは別の話です。 $\psi=1$ かつ $K=0$ なので $(\psi+K)\nabla\cdot u = \nabla\cdot u$ です。液体も自分の圧縮率では圧縮されるべきで、その項がその役目を果たします。相殺が完了するかどうかも体積弾性率の比に依存します。比が $10$ 前後なら(超臨界流体同士、気体と気体)、 $\psi=0.5$ で $\psi+K=0.0909$ です。ゼロではありません。

1 barから10 barまで押すと、どれだけずれるか#

オイラー格子で直接測ります。区間 $[0,1]$ の左半分を水、右半分を空気で埋めます。左は壁( $u=0$ )、右は空気が流入する開放境界です。圧力は空間的に一様なまま1 barから10 barへ線形に上昇させます。Stiffened-gas状態方程式( $\gamma_l=4.4$ 、 $\Pi_l=6\times10^8$ Pa)を使い、面速度は上で導いた混合物の圧縮率から積分します。

水は右境界に達しないので液体質量が保存します。したがって厳密な液体体積は $\psi_0\,\rho_l(p_0)/\rho_l(p)$ ですぐ出ます。比較の基準が無料で手に入ります。

import math
 
P0, PEND = 1.0e5, 1.0e6                    # 1 bar -> 10 bar
GAM_L, PI_L, RHO_L0 = 4.4, 6.0e8, 1000.0   # 水 (stiffened gas)
GAM_G, PI_G, RHO_G0 = 1.4, 0.0,   1.2      # 空気
 
def stiffened_rho(p, gam, pi, rho0):
    return rho0 * ((p + pi) / (P0 + pi)) ** (1.0 / gam)
 
def stiffened_bulk(p, gam, pi):
    return gam * (p + pi)                  # rho*a^2
 
def wood_speed(psi, p):
    rho_m = psi * stiffened_rho(p, GAM_L, PI_L, RHO_L0) \
          + (1 - psi) * stiffened_rho(p, GAM_G, PI_G, RHO_G0)
    comp  = psi / stiffened_bulk(p, GAM_L, PI_L) \
          + (1 - psi) / stiffened_bulk(p, GAM_G, PI_G)
    return math.sqrt(1.0 / (rho_m * comp)), rho_m
 
def compressibility_factor(psi, p, mode):
    bl = stiffened_bulk(p, GAM_L, PI_L)
    bg = stiffened_bulk(p, GAM_G, PI_G)
    k = (bg - bl) * psi * (1 - psi) / (bg * psi + bl * (1 - psi))
    return 0.0 if mode == 'zero' else (-k if mode == 'flip' else k)
 
def squeeze_column(mode, ncell, nstep, psi0=0.5):
    """左壁 u=0、右は純粋な気体が流入、p(t) は 1->10 bar の線形"""
    dx, dt = 1.0 / ncell, 1.0 / nstep
    psi = [1.0 if (i + 0.5) * dx < psi0 else 0.0 for i in range(ncell)]
    for n in range(nstep):
        p    = P0 + (PEND - P0) * (n * dt)
        pdot = PEND - P0
        bl   = stiffened_bulk(p, GAM_L, PI_L)
        bg   = stiffened_bulk(p, GAM_G, PI_G)
        u = [0.0] * (ncell + 1)            # 混合物の圧縮率 -> du/dx
        for i in range(ncell):
            u[i + 1] = u[i] - dx * pdot * (psi[i] / bl + (1 - psi[i]) / bg)
        nxt = psi[:]
        for i in range(ncell):
            ur, ul = u[i + 1], u[i]
            right  = psi[i + 1] if i + 1 < ncell else 0.0
            fr = ur * (psi[i] if ur > 0 else right)
            fl = (ul * (psi[i - 1] if ul > 0 else psi[i])) if i > 0 else 0.0
            src = (psi[i] + compressibility_factor(psi[i], p, mode)) * (ur - ul) / dx
            nxt[i] = min(1.0, max(0.0, psi[i] - dt * (fr - fl) / dx + dt * src))
        psi = nxt
    exact = psi0 * stiffened_rho(P0, GAM_L, PI_L, RHO_L0) \
                 / stiffened_rho(PEND, GAM_L, PI_L, RHO_L0)
    return sum(psi) * dx, exact
 
print("[1] stiffened-gas properties")
for p in (P0, PEND):
    rl = stiffened_rho(p, GAM_L, PI_L, RHO_L0)
    rg = stiffened_rho(p, GAM_G, PI_G, RHO_G0)
    print(f"  p={p/1e5:4.1f} bar  rho_l={rl:8.3f}  rho_g={rg:6.3f}"
          f"  rho_l*a_l^2={stiffened_bulk(p,GAM_L,PI_L):.3e}"
          f"  rho_g*a_g^2={stiffened_bulk(p,GAM_G,PI_G):.3e}")
 
print("\n[2] K and (psi+K) at 1 bar")
for psi in (0.10, 0.25, 0.50, 0.75, 0.90):
    k = compressibility_factor(psi, P0, 'acid')
    a, _ = wood_speed(psi, P0)
    print(f"  psi={psi:4.2f}  K={k:+9.6f}  psi+K={psi+k:+.3e}  a_wood={a:7.2f} m/s")
 
print("\n[3] liquid-volume change after 1->10 bar compression")
print("  N     exact        ACID         K=0          K flipped    | x exact")
for ncell in (50, 100, 200, 400):
    out = {}
    for mode in ('acid', 'zero', 'flip'):
        v, exact = squeeze_column(mode, ncell, 10 * ncell)
        out[mode] = v - 0.5
    d = exact - 0.5
    print(f"  {ncell:4d}  {d:.4e}  {out['acid']:.4e}  {out['zero']:.4e}  "
          f"{out['flip']:.4e} |  {out['zero']/d:5.3f}  {out['flip']/d:5.3f}")

[1] stiffened-gas properties
  p= 1.0 bar  rho_l=1000.000  rho_g= 1.200  rho_l*a_l^2=2.640e+09  rho_g*a_g^2=1.400e+05
  p=10.0 bar  rho_l=1000.341  rho_g= 6.215  rho_l*a_l^2=2.644e+09  rho_g*a_g^2=1.400e+06
 
[2] K and (psi+K) at 1 bar
  psi=0.10  K=-0.099994  psi+K=+5.891e-06  a_wood=  39.23 m/s
  psi=0.25  K=-0.249982  psi+K=+1.767e-05  a_wood=  27.28 m/s
  psi=0.50  K=-0.499947  psi+K=+5.302e-05  a_wood=  23.65 m/s
  psi=0.75  K=-0.749841  psi+K=+1.590e-04  a_wood=  27.32 m/s
  psi=0.90  K=-0.899523  psi+K=+4.770e-04  a_wood=  39.43 m/s
 
[3] liquid-volume change after 1->10 bar compression
  N     exact        ACID         K=0          K flipped    | x exact
    50  -1.7027e-04  -1.7027e-04  -3.1695e-04  -7.1967e-04 |  1.861  4.227
   100  -1.7027e-04  -1.7027e-04  -3.1559e-04  -7.0606e-04 |  1.853  4.147
   200  -1.7027e-04  -1.7027e-04  -3.1256e-04  -6.7737e-04 |  1.836  3.978
   400  -1.7027e-04  -1.7027e-04  -3.0670e-04  -6.2727e-04 |  1.801  3.684

K = 0、符号反転、体積弾性率比のスライダーが、液体柱の減る速さを変えます。

ψ₀ (initial liquid fraction) = 0.50final p/p₀ = 10log₁₀(ρ_l a_l² / ρ_g a_g²) = 4.3

p/p₀ = 1.00  |  ψ = 0.50000  |  K = 0.00000  |  ψ+K = 5.00e-1
V_liquid: scheme 0.50000 vs exact 0.50000  →  +0.00%

K = 0 を選んで再生すると、青いブロックが厳密解より速く減ります。比を $1.0$ まで下げると三本の柱が重なり、 $\psi+K$ の表示は $0$ から離れます。

圧力を10倍にする間に水が実際に失う体積は $1.70\times10^{-4}$ 、初期体積の0.034%です。ACIDクロージャはこの値を小数第5位まで合わせます。 $K=0$ は代わりに $3.13\times10^{-4}$ を削りました。物理的な圧縮量の1.84倍です。

0.034%と0.062%の差は小さく見えます。しかし比較対象は全体積ではなく物理的な圧縮量です。その量の84%をスキームが上乗せしたことになります。気泡崩壊やキャビテーションのように体積変化そのものが答えになる問題では、この比率が答えを分けます。

表[2]はほぼ対称です。 $\psi=0.10$ では $K=-0.099994$ 、 $\psi+K=5.891\times10^{-6}$ 、音速 $39.23$ m/sです。 $\psi=0.90$ では $K=-0.899523$ 、 $\psi+K=4.770\times10^{-4}$ 、音速 $39.43$ m/sです。両端の音速差は $0.20$ m/sなのに、 $\psi+K$ は約 $81$ 倍開きます。底は $\psi=0.50$ の $23.65$ m/sです。 $39.23$ より $15.58$ m/s、 $40\%$ 近く遅い値です。体積弾性率の比 $2.64\times10^{9}/1.40\times10^{5}$ は $18857$ です。その比が $\psi=0.5$ の $K$ を $-0.5$ のすぐ隣、表の $-0.499947$ まで寄せます。 $\psi+K$ は $+5.302\times10^{-5}$ です。表[3]で $K=0$ の倍率は $N=50$ の $1.861$ から $N=400$ の $1.801$ へ、 $0.060$ 減ります。符号を反転した倍率は $4.227$ から $3.684$ へ、 $0.543$ 減ります。同じ $8$ 倍の細分化なのに、符号側の減少は $0.060$ の約 $9$ 倍です。 $N=400$ でも $1.801$ と $3.684$ は $1$ から遠いです。 $N=50$ で $3.1695\times10^{-4}$ と $1.7027\times10^{-4}$ の差は $1.4668\times10^{-4}$ で、本文の $1.467\times10^{-4}$ に丸められます。 $N=400$ は $3.0670\times10^{-4}$ から $1.7027\times10^{-4}$ を引いて $1.3643\times10^{-4}$ です。 $10$ barで液体密度は $1000.000$ から $1000.341$ へ $0.341$ 上がり、気体は $1.200$ から $6.215$ へ跳ねます。無い $K$ は、その小さな液体圧縮の $1.86$ 倍を界面に請求します。

格子を8倍細かくすれば誤差は消えるか#

表[3]の最後の2列がその答えです。 $N$ を50から400へ8倍にする間、 $K=0$ の比は1.861から1.801へ動きました。誤差の大きさでは $1.467\times10^{-4}$ から $1.364\times10^{-4}$ 、7%の減少です。収束とは呼べません。

理由は $K=0$ が離散化誤差ではないからです。 $K=0$ とすると解いている方程式自体が $\mathrm{D}\psi/\mathrm{D}t=0$ に変わります。別の方程式の解に収束しているのです。 $\Delta x\to0$ で差が消える理由はありません。

残る7%は、なまされた界面の帯の形が格子で少し変わるためです。誤差の大きさはその帯の中で $\psi\,\nabla\cdot u$ を積分した量で、1次風上差分の数値拡散は格子を細かくしても帯をセル数で広げます。帯の物理的な厚みがほぼ変わらないので、積分値もほぼ変わりません。その帯を1ステップでどれだけ押せるかは、VOF界面移流のCFL上限が別に決めます。

この性質は悪い知らせであり、よい知らせでもあります。格子収束性の試験ではこのバグは見つかりません。逆に、格子を2つ回して比が一定のままなら、それ自体がクロージャが誤っているという信号です。

符号を反転するとどんな形で誤るか

3列目が符号を反転した結果です。 $K\to-K$ なら係数は $\psi+K \to 2\psi$ になり、圧縮量は物理値の3.68〜4.23倍に跳ねます。向きも反転します。圧縮しているのに液体の体積率が下がるので、画面上は気泡が成長するように見えます。

厄介なのは、こちらもそれらしい絵を出すところです。界面はなめらかなまま、体積率は $[0,1]$ の内側にとどまり、発散もしません。単相極限を確認しても捕まりません。 $\psi\to0$ と $\psi\to1$ では $K$ がどちらの符号でもゼロだからです。

ですから論文の式をコードに移すときは、 $\psi$ がどちらの相の体積率かをコメントで固定し、上の1次元圧縮テスト一つで符号を確定させるのが早いです。50セルなら1秒で回ります。

同じ仮定がエンタルピーでもう一度使われる場所

論文のACID(Acoustically-Conservative Interface Discretisation)は、この等圧クロージャを密度とエンタルピーにそれぞれ適用します。密度は簡単な側です。部分密度を線形補間すると

\rho = \rho_g + \psi(\rho_l - \rho_g)

となり、これが等圧クロージャと数学的に同一です。質量と運動量の保存がそのまま従います。

エンタルピーは違います。 $h$ は線形系の中で直接解く1次変数なので、密度のように補間式を差し替えられません。そこでACID仮定を適用したターゲット面エンタルピーを別に計算します。

\hat{h}_f = \frac{1}{\tilde{\rho}_f}\left[\rho_U h_U + \xi_f L_f (\rho_D h_D - \rho_U h_U)\right]

$U$ と $D$ は上流・下流セル、 $L_f$ は補間重み、 $\xi_f$ は制限関数です。そして現在値との差を遅延補正項としてエネルギー方程式の右辺に加えます。

\delta h_f = \hat{h}_f - \tilde{h}_f

$K$ と $\delta h_f$ は同じ仮定の二つの顔です。体積率の式では代数的に閉じて係数一つに落ち、1次変数であるエンタルピーでは閉じずに陽的なソース項へ先送りされます。先送りされた側は反復の中で収束しなければならないので、非線形反復の回数に影響します。等圧で一つの圧力に貼り付ける仕上げは、二つの圧力を解いてから一つの圧力へ緩和する経路の終点と同じです。違いは、こちらがその緩和を $K$ 一つに折り畳む点です。

Dennerがリーマンソルバーを捨てられた理由#

密度ベースの解法は界面をリーマン問題として解き、接触不連続面を保ちます。圧力ベース完全連成の解法は $\rho$ を解析変数から外し、 $u$ 、 $p$ 、 $h$ を一つの線形系に束ねます。同一位置格子で生じる圧力-速度のデカップリングはMWIで抑えます。極端な密度比では調和平均した面密度を重みに使います。接触面を再構成で守る道は、より小さい振幅を選ぶMUSCL-THINC-BVDです。

リーマンソルバーを外すと、界面の物理を自動で面倒みてくれる装置も一緒に外れます。その穴はクロージャ仮定で手作業で埋めるしかありません。 $K$ と $\delta h_f$ がその埋め物です。

混合物音速23.65 m/sを陰的に飛び越えながら接触不連続面をなまさないためには、そのクロージャが正確である必要があります。音響を陰的にすると、CFLの請求書は輸送側へ移ります。音響と輸送を分けてCFLを解く分離がその場所です。1.86倍は、そのクロージャが果たしている仕事の大きさを測った数字です。

Kを消したら水が1.86倍圧縮された — 格子を細かくしても減らない界面の誤差

圧縮性VOFに −(ψ+K)∇⋅u-(\psi+K)\nabla\cdot u−(ψ+K)∇⋅u が付いている理由#

KKK はどこから出てくる数字か#

ψ+K\psi+Kψ+K が 5×10−55\times10^{-5}5×10−5 なら両方の項を消してよいか#