音速1625 m/sの水が圧力部分系では18 m/sだった — Toro–Vázquez分割が剛性を移す場所

音速1625 m/sの水が18 m/sと計算された#

水の音速はおよそ1625 m/sです。ところが今回の論文の圧力部分系では、1 barの水の波速が18.44 m/sになります。計算ミスではありません。一つの式からそのまま出てくる値です。

読むのはSchropff、Petitpas、Danielの2026年JCP論文です。題名は"Semi-implicit method for two-phase flows in pressure disequilibrium"です。圧力の異なる二相を一つの速度で解くIMEX(陰的・陽的の混合)手法です。物質をWaterやSpinelに切り替えるとオレンジの線がほぼ垂直になり、圧力とマッハ数のスライダーが $A$ と時間刻みの比を一緒に動かします。

灰色の破線は元の系の特性線 $u\pm c$ です。オレンジの実線は圧力部分系の $\pm A$、青い実線は移流速度 $u$ です。右側の棒は、低マッハで時間刻みが何倍になるかを対数目盛で示します。水でオレンジの線が垂直に近いほど、$A/c$ は0.01を下回っています。

1993年、圧力を対流から切り離した分割#

出発点はZha–Bilgen分割(flux splitting、流束を性質ごとに分けて別々に扱う手法)です。オイラー方程式の流束を二つに分けます。一つは流速で運ばれる対流部分です。もう一つは圧力が押す部分です。

ToroとVázquezは2012年にこの分割をGodunov型に書き直しました。各部分系に別々にリーマン問題を立てられるようにしたのです。Schropffらはこの考えを圧力非平衡の二相モデルへ持ち込みました。圧力の効果を含む流束と非保存項を第一の部分系に集めます。残る輸送項は第二の部分系に入ります。

2026年、二つの相を二つの部分系に分ける#

原始変数は $W=(\alpha_1,\alpha_2,\rho_1,\rho_2,p_1,p_2,u)^T$ です。圧力部分系は次のとおりです。

$\alpha_k$ は相 $k$ の体積分率、$\rho_k$ は相の密度、$e_k$ は比内部エネルギーです。$p=\sum_k\alpha_k p_k$ は混合圧力です。

移流部分系には圧力がありません。体積分率、質量、内部エネルギー、運動量がすべて $u$ で運ばれます。そのため7個の固有値はすべて $u$ です。二つの部分系を足すと元のモデルに戻ります。緩和項は解析から外しています。

圧力段階では密度が凍結している

圧力部分系の式に戻ります。$\alpha_k$ と $\alpha_k\rho_k$ は時間に対して一定です。したがって $\rho_k$ も動きません。圧力を変える経路は内部エネルギーだけです。原始変数で書くとこうなります。

$\Gamma_k$ はグリュナイゼン係数(Grüneisen parameter、密度一定でエネルギーを加えたときの圧力上昇率)です。$e_{p,k}$ は密度一定での $\partial e_k/\partial p_k$ です。この式と運動量式を組み合わせると、波速 $\pm A$ が出てきます。

ここで $\rho=\sum_k\alpha_k\rho_k$ は混合密度です。$A$ には $u$ が含まれません。そのため論文は、この波の構造が常に亜音速だと述べています。

硬化気体(stiffened gas)の状態方程式なら $\Gamma_k=\gamma_k-1$ です。すると論文の式(14)になります。

$c_f$ は凍結音速(frozen sound speed、相間緩和がないときの混合音速)、$p_{\infty,k}$ は硬化気体の剛性圧力です。理想気体なら $A^2=c^2-p/\rho$ で、空気は $A/c=\sqrt{(\gamma-1)/\gamma}=0.535$ になります。水は事情が違います。$\gamma p_\infty$ の項が $c^2$ のほぼすべてを占め、その分がまるごと抜け落ちます。

4物質のAとcを一つの表に並べた#

論文の表1の硬化気体パラメータ(空気、水、エポキシ、スピネル)で $A$ と $c$ を計算しました。同じコードで空気の圧力部分系だけを取り出し、パルス実験も行います。陽的Lax–Friedrichsと陰的後退オイラーを同じCourant数で比べます。最後は低マッハでの時間刻みの比です。

import math
 
# Stiffened-gas parameters (gamma, p_inf [Pa]) from the paper's Table 1
MATS = {"Air": (1.4, 0.0), "Water": (4.4, 6e8), "Epoxy": (2.43, 5.3e9), "Spinel": (1.62, 141e9)}
RHO = {"Air": 1.2, "Water": 1000.0, "Epoxy": 1185.0, "Spinel": 3622.0}
 
def split_speeds(gam, pinf, p, rho):
    c = math.sqrt(gam * (p + pinf) / rho)     # full sound speed
    A = math.sqrt((gam - 1.0) * p / rho)       # pressure-subsystem speed, Eq. (14)
    return c, A
 
def thomas_solve(rhs, beta):
    # (I - beta*L) x = rhs, compact Laplacian, zero-gradient ends
    n = len(rhs)
    cp, dp = [0.0] * n, [0.0] * n
    b0 = 1 + beta
    cp[0], dp[0] = -beta / b0, rhs[0] / b0
    for i in range(1, n):
        b = 1 + beta if i == n - 1 else 1 + 2 * beta
        m = b + beta * cp[i - 1]
        cp[i] = -beta / m
        dp[i] = (rhs[i] + beta * dp[i - 1]) / m
    x = [0.0] * n
    x[-1] = dp[-1]
    for i in range(n - 2, -1, -1):
        x[i] = dp[i] - cp[i] * x[i + 1]
    return x
 
def ddx(f, dx):
    n = len(f)
    return [(f[min(i + 1, n - 1)] - f[max(i - 1, 0)]) / (2 * dx) for i in range(n)]
 
def pulse_run(courant, implicit, n=200, travel=0.3):
    # linearised pressure subsystem (air, 1 bar): p_t + K u_x = 0, u_t + p_x / rho = 0
    rho, K = 1.2, 0.4 * 1e5
    A = math.sqrt(K / rho)
    dx = 1.0 / n
    x = [(i + 0.5) * dx for i in range(n)]
    p = [100.0 * math.exp(-((xi - 0.5) / 0.04) ** 2) for xi in x]
    u = [0.0] * n
    dt = courant * dx / A
    steps = math.ceil(travel / A / dt)
    for _ in range(steps):
        if implicit:
            ux = ddx(u, dx)
            p = thomas_solve([p[i] - dt * K * ux[i] for i in range(n)], dt * dt * K / rho / dx ** 2)
            px = ddx(p, dx)
            u = [u[i] - dt / rho * px[i] for i in range(n)]
        else:  # Lax-Friedrichs
            P = [p[0]] + p + [p[-1]]; U = [u[0]] + u + [u[-1]]
            p = [0.5 * (P[i] + P[i + 2]) - dt * K * (U[i + 2] - U[i]) / (2 * dx) for i in range(n)]
            u = [0.5 * (U[i] + U[i + 2]) - dt / rho * (P[i + 2] - P[i]) / (2 * dx) for i in range(n)]
        if max(abs(v) for v in p) > 1e9:
            return steps, float("inf"), float("nan")
    right = max(range(n // 2, n), key=lambda i: p[i])
    speed = (x[right] - 0.5) / (steps * dt)
    return steps, max(abs(v) for v in p[n // 2:]), speed
 
print("material  p[bar]     c[m/s]     A[m/s]      A/c")
for name, (g, pinf) in MATS.items():
    for pbar in (1.0, 100.0):
        rho = RHO[name] * (pbar if name == "Air" else 1.0)  # isothermal air
        c, A = split_speeds(g, pinf, pbar * 1e5, rho)
        print(f"{name:7s} {pbar:7.0f} {c:10.1f} {A:10.2f} {A / c:9.5f}")
 
print("\nwater/air mixture, 1 bar, Eq. (13)")
for aw in (0.01, 0.5, 0.99):
    rho = aw * 1000.0 + (1 - aw) * 1.2
    A = math.sqrt((aw * 1e5 * 3.4 + (1 - aw) * 1e5 * 0.4) / rho)
    print(f"alpha_w={aw:4.2f}  rho={rho:7.1f}  A={A:7.2f} m/s")
 
A_air = math.sqrt(0.4e5 / 1.2)
print(f"\npressure pulse, air 1 bar, A = {A_air:.1f} m/s, c = {math.sqrt(1.4e5 / 1.2):.1f} m/s")
print("Courant_A steps  exp_peak  exp_speed  imp_peak  imp_speed")
for cr in (0.5, 0.9, 1.5, 5.0, 20.0):
    s, pe, ve = pulse_run(cr, False)
    _, pi, vi = pulse_run(cr, True)
    print(f"{cr:7.1f} {s:6d} {pe:9.2f} {ve:9.1f} {pi:9.2f} {vi:9.1f}")
 
print("\nlow-Mach air: dt ratio at equal Courant number = (|u|+c)/|u|")
for M in (0.3, 0.1, 0.01, 0.001):
    print(f"M={M:6.3f}  ratio={(1 + M) / M:9.1f}")

material  p[bar]     c[m/s]     A[m/s]      A/c
Air           1      341.6     182.57   0.53452
Air         100      341.6     182.57   0.53452
Water         1     1624.9      18.44   0.01135
Water       100     1638.3     184.39   0.11255
Epoxy         1     3296.7      10.99   0.00333
Epoxy       100     3299.8     109.85   0.03329
Spinel        1     7941.3       4.14   0.00052
Spinel      100     7941.6      41.37   0.00521
 
water/air mixture, 1 bar, Eq. (13)
alpha_w=0.01  rho=   11.2  A=  62.00 m/s
alpha_w=0.50  rho=  500.6  A=  19.48 m/s
alpha_w=0.99  rho=  990.0  A=  18.45 m/s
 
pressure pulse, air 1 bar, A = 182.6 m/s, c = 341.6 m/s
Courant_A steps  exp_peak  exp_speed  imp_peak  imp_speed
    0.5    120     25.59     184.1     35.90     181.1
    0.9     67     42.37     183.2     30.52     180.2
    1.5     40    240.14     129.3     25.76     178.0
    5.0     12 102438.55      13.7     16.02     168.9
   20.0      3   3336.76      32.0      9.82     129.3
 
low-Mach air: dt ratio at equal Courant number = (|u|+c)/|u|
M= 0.300  ratio=      4.3
M= 0.100  ratio=     11.0
M= 0.010  ratio=    101.0
M= 0.001  ratio=   1001.0

Courantのスライダーを0.5から20へ動かし、陽解法と陰解法を切り替えると、表の exp_peak と imp_peak が画面で計算し直されます。

0.9まではどちらも波速が180 m/s付近にあります。1.5を超えると陽解法の山は240 Paまで跳ね、陰解法は20でも9.82 Paまで平たくなるだけで、数値は飛びません。

1 barの水で $A/c$ は0.01135です。エポキシは0.00333、スピネルは0.00052まで下がります。圧力を100 barに上げると、水の $A$ はちょうど10倍になります。$A^2$ が $p$ に比例するからです。一方 $c$ は1624.9から1638.3へ、0.8%しか動きません。空気は圧力に関係なく0.53452で固定です。

水と空気の混合物も押さえておきます。水の体積分率0.5で $A$ は19.48 m/sです。0.99のときの18.45 m/sとほぼ同じです。混合密度も分子 $\sum\alpha_k p_k\Gamma_k$ も水が決めるからです。空気が99%になって初めて62.00 m/sに跳ね上がります。

1 barから100 barへ上げると、水の $A$ は18.44から184.39になります。ちょうど10倍です。$A^2$ が $p$ に比例するので、圧力が100倍なら波速は10倍です。同じ区間で $c$ は1624.9から1638.3へ、0.8%しか動きません。エポキシの $A$ も10.99から109.85へ10倍です。$c$ は3296.7から3299.8に留まります。スピネルは $A$ が4.14から41.37、$c$ は7941.3と7941.6です。$A/c$ は0.00052から0.00521へ増えても0.01を超えません。

空気は1 barと100 barが同じ行です。$c=341.6$、$A=182.57$、$A/c=0.53452$。理想気体なので密度を圧力に比例させ、温度が固定されると $c$ も固定されます。$\sqrt{(\gamma-1)/\gamma}=\sqrt{0.4/1.4}$ が0.535と書かれる理由と同じです。

水が少しでもあると $A$ は水の値に張り付きます。$\alpha_w=0.01$ では $\rho=11.2$、$A=62.00$ です。$\alpha_w=0.50$ では $\rho=500.6$、$A=19.48$。$\alpha_w=0.99$ では $\rho=990.0$、$A=18.45$ で、純水の18.44と0.01 m/s差です。分母の混合密度と分子の $\sum \alpha_k p_k \Gamma_k$ を水が一緒に持っていくためです。パルス表の陽的な山は Courant 0.9 で42.37 Pa、1.5で240.14 Pa、5.0で102438.55 Paです。陰的な側は同じ三つの欄で30.52、25.76、16.02 Paとむしろ減ります。Courant 20の陰的な山9.82 Paは、右向き波の基準50 Paの1/5で、測定速度は129.3 m/sまで下がります。

剛性は消えず、段階の間へ移る

では水の1625 m/sはどこへ行ったのでしょうか。二つの部分系のどちらにも $u\pm c$ はありません。圧力部分系は $\pm A$、移流部分系は $u$ だけです。

本来の音速は、二つの段階を交互に回すときに戻ってきます。移流段階で $\rho_k$ が変わります。状態方程式がその変化を $\rho c^2$ の比率で圧力変化に変えます。その圧力が次の圧力段階の入力になります。水の剛性は部分系の固有値から抜け、段階間の状態方程式の再計算へ移ったのです。音響を一つの段階、輸送を次の段階へ分けて剛性を移す解き方は、音響と輸送を段階で分けたIMEXにもあります。

ここから先は論文の記述ではなく、筆者の解釈です。圧力段階だけを切り出すと、$\Delta t\le\Delta x/A$ で陽的に扱えそうに見えます。しかし合成されたスキームの安定性をこの不等式は保証しません。$A$ が小さい水でも論文が圧力段階を常に陰的に解くのは、この文脈では自然です。

圧力段階を陰的に解くとパルスは生き残り、平たくなる

出力の二つ目のブロックがパルス実験です。陽的Lax–FriedrichsはCourant 0.9までは持ちこたえ、1.5で崩れます。1.5では最大圧力が240 Paと初期値を超えました。5.0では10万Paを超えました。陰的スキームは20でも発散しません。その代わり最大値は9.82 Paと、厳密解50 Paの1/5まで平たくなりました。

この減衰がIMEXの契約条件です。低マッハ問題では音響波の振幅は関心の対象ではありません。圧力場が速度場にすぐ追従すれば十分です。音響波を正確に追う必要がある問題なら、この取引は損になります。

論文は陰的圧力段階で相ごとの圧力方程式を立てます。$\Gamma_k$ が圧力の関数なので、式は非線形です。これをPicard反復(一部の係数を直前の反復値に固定して線形化する反復法)で解きます。空間離散化はコンパクトなステンシルなので、行列は帯幅の狭い帯行列になります。上のコードの thomas_solve は、1Dでその帯行列を一つ解いていることになります。

低マッハでΔtが1001倍になる計算#

陽的スキームの時間刻みは $\Delta x/(|u|+c)$ に縛られます。IMEXでは移流段階だけが陽的なので、$\Delta x/|u|$ まで緩みます。論文は $\Delta t_{IMEX}=0.5\,\Delta x/|\sigma|$ を使います。同じCourant数を仮定すると、比は次のとおりです。

$M$ はマッハ数です。$M=0.3$ なら4.3倍、$M=0.001$ なら1001倍です。ただしこれは1ステップのコストを同じとみなした上限です。陰的段階はPicard反復のたびに線形系を解きます。論文の結論も同じことを述べています。低マッハ問題では利得が大きく、衝撃波管問題では小さいか、ないこともあります。

非保存項が残したエネルギーを戻す最後の段階

分割には後始末が一つ付いてきます。相ごとの内部エネルギー式には $\alpha_k p_k\nabla\cdot u$ のような非保存項があります。これを離散化すると、全エネルギーを厳密には保存できません。

そこで論文は補正演算子を第三の段階として置きます。質量、運動量、混合全エネルギー、体積分率は予測値をそのまま使います。相ごとの内部エネルギーにだけ $\alpha_k^n\,\delta\varepsilon$ を加えます。混合全エネルギーと一致するように差を配分するのです。圧力緩和はこのIMEX段階が終わった後に別途行います。

水の $A$ だけを見て圧力段階を陽的に置くと、合成スキームは持たない#

スピネルの $A/c$ は0.00052、エポキシは0.00333、1 barの水は0.01135です。三つとも0.01付近か、その下です。圧力段階だけを見ると $\Delta t\le\Delta x/A$ で陽的に扱えそうです。その不等式は部分系の話です。合成されたスキームの剛性は $A$ ではなく $c$ にあります。1 barの水では18.44と1624.9は別の数字です。$M=0.3$ なら時間刻みの比は4.3倍、$M=0.001$ なら1001倍です。1Dで $1+1/M$ がPicard反復のコストを超えるときだけ、この比は利得として残ります。

Courant 20の後退オイラーは、この出力でパルスの山を50 Paの1/5である9.82 Paまで削ります。音響振幅が結果に入るなら、その数字を書く必要があります。陰的な圧力段階が $L_\infty$ の山を平たくする場所は、陰的圧力段階が $L_\infty$ の山を平たくする場合と同じです。振幅を解析に使うなら、時間刻みを音響に合わせるか、後退オイラーより高次の陰的積分へ移ります。