[論文レビュー] 水と空気を半分ずつ混ぜたら音は24 m/sで這った

体積分率の式に付くK∇·u項こそが5方程式モデルにWood音速で音を伝えさせる項であり、これを消すと空気50%で音速が31倍速くなります。

1930年の教科書にある一行の式#

1930年、英国海軍研究所のA. B. Woodは音響学の教科書を出版しました。その中に、気泡を含む水の音速を与える一行の式があります。混合物の圧縮率は各相の圧縮率の体積加重和になる、という式です。これに従うと、空気と水を半分ずつ混ぜた泡の音速は約24 m/sです。水(1625 m/s)よりも、空気(342 m/s)よりもはるかに遅い値です。

本記事ではMurrone–GuillardのINRIA報告書RR-4778(2003)を読みます。扱うのは圧縮性二相流の5方程式縮約モデル(reduced model、7方程式モデルから緩和の極限で導くモデル)です。問いは一つです。体積分率の式に付く $K\nabla\cdot u$ 項は何をしているのか。その答えをPythonの1次元計算で確かめます。同じ管、同じパルスに $K$ を入れた場合と外した場合で、音の速さを測ります。

速度を揃えたあと圧力緩和で式を減らす道はSaurel(2009)の圧力緩和にあります。今回は、その緩和を最後まで進めたとき体積分率の式に残る $K$ を測ります。

αスライダーを動かすと、管AのWoodパルスと管Bの速いパルスが同じ画面で分かれます。

air volume fraction α = 0.500slow motion ×20

K = 0.4999 (≈ α when water is stiff)physical time 0.0 msTube A laps: 0Tube B laps: 0 (×31.5 faster)

αスライダーを0.5にすると、管Aのパルスは這うように進み、管Bのパルスは何十周もします。αを1%付近に動かすと、青い曲線がすでに水の音速の1/14まで落ちているのが分かります。

2003年の報告書が問い直したこと#

Woodの式自体は古い結果です。問題は、数値モデルがその式を自ら再現するかどうかです。7方程式のBaer–Nunziatoモデルは両相に別々の速度と圧力を持たせます。物理的には豊かですが計算コストが大きいモデルです。両相がすばやく同じ速度・同じ圧力に揃うなら、変数を減らせます。

Kapilaら(2001)はこの極限で5方程式モデルを導きました。MurroneとGuillardは同じモデルを漸近展開(小さなパラメータで解を級数展開する手法)で導き直しました。さらに双曲性(実数の有限な波速を持つ性質)と波の構造を整理し、二つの数値スキームを試しています。式をさらに減らして界面平衡だけを残す選択は4方程式モデルが取った側です。今回は5本の式を保ったまま、 $K$ の一行だけを入れたり外したりします。

報告書の一枚の図が今回の出発点です。空気–水混合物の音速を空気体積分率に対して描いた図です。曲線は二本あります。"5 Eq. Reduc"は縮約モデル、"5 Eq. Trans"は体積分率を輸送するだけのモデルです。二本の曲線は大きく離れています。

緩和を無限に速くすると残る5本の式#

7方程式モデルの速度緩和と圧力緩和の時間スケールを0にすると、次が残ります。

\begin{aligned} &\partial_t(\alpha_1\rho_1) + \nabla\cdot(\alpha_1\rho_1 u) = 0, \qquad \partial_t(\alpha_2\rho_2) + \nabla\cdot(\alpha_2\rho_2 u) = 0,\\ &\partial_t(\rho u) + \nabla\cdot(\rho u\otimes u + p I) = 0, \qquad \partial_t(\rho E) + \nabla\cdot\big((\rho E + p)u\big) = 0,\\ &\partial_t\alpha_1 + u\cdot\nabla\alpha_1 = K\,\nabla\cdot u . \end{aligned}

$\alpha_k$ は相 $k$ の体積分率、 $\rho_k$ は相密度、 $u$ と $p$ は共通の速度と圧力、 $E$ は混合物の全エネルギーです。最初の4本は保存形です。最後の式だけが非保存形で、右辺の $K\nabla\cdot u$ が今回の主役です。

K∇·uはどこから来るのか#

圧力平衡なら、両相は同じ圧力変化を受けます。衝撃波のない区間では相 $k$ は等エントロピー的に圧縮されるので、 $D\ln\rho_k/Dt = (\rho_k c_k^2)^{-1} Dp/Dt$ です。これに相の質量保存 $D(\alpha_k\rho_k)/Dt = -\alpha_k\rho_k\nabla\cdot u$ を組み合わせます。

\frac{D\alpha_k}{Dt} = -\alpha_k\nabla\cdot u - \frac{\alpha_k}{\rho_k c_k^2}\frac{Dp}{Dt}

二相の式を足すと $\alpha_1+\alpha_2=1$ より左辺は0になります。すると圧力の式が出てきます。

\frac{Dp}{Dt} = -\rho c_W^2\,\nabla\cdot u, \qquad \frac{1}{\rho c_W^2} = \frac{\alpha_1}{\rho_1 c_1^2} + \frac{\alpha_2}{\rho_2 c_2^2}

$c_W$ がWood音速です。これを最初の式に戻すと $K$ が決まります。

K = \frac{\alpha_1\alpha_2\,(\rho_2 c_2^2 - \rho_1 c_1^2)}{\alpha_1\rho_2 c_2^2 + \alpha_2\rho_1 c_1^2}

$\rho_k c_k^2$ は相 $k$ の体積弾性率(bulk modulus、圧縮への抵抗)です。水は空気より約1.9万倍硬いです。そのため空気を相1とすると $K\approx\alpha_1$ になります。混合物が圧縮されると、その体積変化をほぼすべて空気が引き受けるという意味です。

Wood音速:圧縮率は体積で、慣性は密度で#

Wood音速が遅いのは、二つの性質が別々の相から来るからです。圧縮率は柔らかい空気が決めます。慣性(密度)は重い水が決めます。空気50%なら圧縮率は純粋な空気の約半分です。密度は空気の約400倍です。音速は $\sqrt{1/(\text{密度}\times\text{圧縮率})}$ なので、空気の音速の $1/\sqrt{200}$ 程度まで下がります。

$K=0$ にすると話が変わります。体積分率は流れに乗って運ばれるだけです。このとき混合物は一つのstiffened gas(水のような硬い液体を近似する状態方程式)のように振る舞います。音速は $\rho\xi c^2 = \sum_k \alpha_k\xi_k\rho_k c_k^2$ 、 $\xi_k = 1/(\gamma_k-1)$ で与えられます。水の大きな体積弾性率がエネルギー加重で混ざり込みます。そのため空気が半分でも音は速いままです。

コード:同じ管に二つのモデルを載せる

状態方程式は空気 $\gamma=1.4,\ \pi=0$ 、水 $\gamma=4.4,\ \pi=6\times10^8$ Paです。1 mの管に空気50%の混合物を満たします。中央に振幅0.1%のガウス型圧力パルスを置きます。右向きパルスが0.3 m進む時間だけ計算し、ピーク位置から速度を測ります。保存変数にはRusanovフラックス、体積分率には中心差分と同じ大きさの人工拡散を使います。

import math
 
# stiffened gas: 相 1 = 空気, 相 2 = 水
G = (1.4, 4.4)
PI = (0.0, 6.0e8)
P0 = 1.0e5
RHO0 = (1.2, 1000.0)
 
def wood_speed(a1):
    rho = a1 * RHO0[0] + (1 - a1) * RHO0[1]
    comp = a1 / (G[0] * (P0 + PI[0])) + (1 - a1) / (G[1] * (P0 + PI[1]))
    return 1.0 / math.sqrt(rho * comp)
 
def isobaric_speed(a1):
    """alpha を輸送するだけの場合 (K = 0): 混合物は一つの stiffened gas"""
    rho = a1 * RHO0[0] + (1 - a1) * RHO0[1]
    xi = a1 / (G[0] - 1) + (1 - a1) / (G[1] - 1)
    gpi = a1 * G[0] * PI[0] / (G[0] - 1) + (1 - a1) * G[1] * PI[1] / (G[1] - 1)
    gam = 1 + 1 / xi
    return math.sqrt(gam * (P0 + gpi / xi / gam) / rho)
 
def mix_pressure(a1, rhoe):
    a2 = 1 - a1
    num = rhoe - a1 * G[0] * PI[0] / (G[0] - 1) - a2 * G[1] * PI[1] / (G[1] - 1)
    return num / (a1 / (G[0] - 1) + a2 / (G[1] - 1))
 
def kapila_K(a1, p):
    k1, k2 = G[0] * (p + PI[0]), G[1] * (p + PI[1])
    a2 = 1 - a1
    return a1 * a2 * (k2 - k1) / (a1 * k2 + a2 * k1)
 
def pulse_speed(a_bg, use_K, n=400, travel=0.3, amp=1e-3):
    dx = 1.0 / n
    x = [(i + 0.5) * dx for i in range(n)]
    a1 = [a_bg] * n
    p = [P0 * (1 + amp * math.exp(-((xi - 0.5) / 0.03) ** 2)) for xi in x]
    m1 = [a_bg * RHO0[0]] * n
    m2 = [(1 - a_bg) * RHO0[1]] * n
    mom = [0.0] * n
    E = [a_bg * (pi + G[0] * PI[0]) / (G[0] - 1) + (1 - a_bg) * (pi + G[1] * PI[1]) / (G[1] - 1) for pi in p]
    c = wood_speed(a_bg) if use_K else isobaric_speed(a_bg)
    s = 1.05 * c                        # Rusanov 速度 = そのモデル自身の音速
    t_end = travel / c
    steps = math.ceil(t_end / (0.4 * dx / s))
    dt = t_end / steps
    g = lambda q, i: q[min(max(i, 0), n - 1)]   # 勾配ゼロのゴーストセル
    for _ in range(steps):
        u = [mom[i] / (m1[i] + m2[i]) for i in range(n)]
        pr = [mix_pressure(a1[i], E[i] - 0.5 * mom[i] * u[i]) for i in range(n)]
        U = (m1, m2, mom, E)
        F = ([m1[i] * u[i] for i in range(n)], [m2[i] * u[i] for i in range(n)],
             [mom[i] * u[i] + pr[i] for i in range(n)], [(E[i] + pr[i]) * u[i] for i in range(n)])
        new = []
        for q, f in zip(U, F):
            fl = [0.5 * (g(f, i - 1) + g(f, i)) - 0.5 * s * (g(q, i) - g(q, i - 1)) for i in range(n + 1)]
            new.append([q[i] - dt / dx * (fl[i + 1] - fl[i]) for i in range(n)])
        na = []
        for i in range(n):
            K = kapila_K(a1[i], pr[i]) if use_K else 0.0
            dudx = (g(u, i + 1) - g(u, i - 1)) / (2 * dx)
            dadx = (g(a1, i + 1) - g(a1, i - 1)) / (2 * dx)
            lap = g(a1, i + 1) - 2 * a1[i] + g(a1, i - 1)
            na.append(a1[i] + dt * (-u[i] * dadx + K * dudx + s * lap / (2 * dx)))
        m1, m2, mom, E = new
        a1 = na
    pr = [mix_pressure(a1[i], E[i] - 0.5 * mom[i] ** 2 / (m1[i] + m2[i])) for i in range(n)]
    ipk = max((i for i in range(n) if x[i] > 0.6), key=lambda i: pr[i])
    return c, (x[ipk] - 0.5) / t_end, steps
 
print("alpha_air  wood[m/s]  K=0[m/s]")
for a in [0.0, 1e-4, 1e-3, 0.01, 0.1, 0.5, 0.9, 0.99, 1.0]:
    print(f"{a:9.4f}  {wood_speed(a):9.1f}  {isobaric_speed(a):8.1f}")
amin = min((k / 1000 for k in range(1, 1000)), key=wood_speed)
print(f"Wood minimum: {wood_speed(amin):.2f} m/s at alpha_air = {amin:.3f}")
print(f"K at alpha_air=0.5: {kapila_K(0.5, P0):.6f}")
print("\n1D pulse, alpha_air = 0.5, N = 400")
for useK in (True, False):
    c, meas, st = pulse_speed(0.5, useK)
    tag = "with K" if useK else "K = 0 "
    print(f"{tag}  theory {c:7.2f}  measured {meas:7.2f}  ratio {meas / c:.3f}  steps {st}")

alpha_air  wood[m/s]  K=0[m/s]
   0.0000     1624.9    1624.9
   0.0001      956.6    1624.3
   0.0010      364.8    1618.9
   0.0100      118.6    1567.2
   0.1000       39.4    1228.3
   0.5000       23.6     745.3
   0.9000       39.2     581.7
   0.9900      112.4     540.9
   1.0000      341.6     341.6
Wood minimum: 23.65 m/s at alpha_air = 0.501
K at alpha_air=0.5: 0.499947
 
1D pulse, alpha_air = 0.5, N = 400
with K  theory   23.65  measured   23.55  ratio 0.996  steps 316
K = 0   theory  745.30  measured  742.19  ratio 0.996  steps 315

上のコードをブラウザで直接動かしてみましょう。同じRusanov計算が二本の管で1ステップずつ進みます。

cells N = 400air fraction α = 0.50steps per frame = 2

model	theory c	measured	ratio	step
with K	23.65	—	—	0 / 0
K = 0	745.30	—	—	0 / 0

measured speed = (peak position − 0.5 m) / t, shown after τ > 0.25

既定値(N = 400, α = 0.5)でPlayを押すと、表のmeasuredは23.55と742.19で終わります。Nを100に下げると二つの比は一緒に下がりますが、モデル間の31倍の差は変わりません。αを動かせば、出力表にない混合比のパルス速度を自分で測れます。

23.55と742.19が語ること#

二つの計算は同じ格子、同じフラックス、同じ初期パルスを使います。違うのは $K$ の一行だけです。それなのにパルス速度は23.55 m/sと742.19 m/sで31.5倍も分かれます。どちらもそれぞれのモデルの理論音速と0.4%以内で一致します。残る0.4%は1次Rusanovの拡散がピークを後ろに引きずる分です。

つまり数値スキームが間違っているのではありません。二つのモデルが別の物理を解いているのです。 $K=0$ モデルはこのパルスを正確に、しかし間違った速度で運びます。格子をいくら細かくしても742が23には下がりません。 $K$ を消して界面付近の水を過剰に圧縮するのは、同じ項の別の顔です。今回の数字は、一様な混合領域にその項が残す痕跡です。

表の二列を $\alpha = 0.5$ で割ると $745.3/23.6 \approx 31.6$ です。パルスの測定は $742.19/23.55 \approx 31.5$ で、ほぼ同じです。各測定は理論の 0.996 倍です。23.65 から 23.55 を引くと 0.10 m/s、745.30 から 742.19 を引くと 3.11 m/s です。絶対誤差は速い側が大きいですが、比率は同じです。1次 Rusanov がピークを同じ比率だけ後ろへ引いています。ステップ数は 316 と 315 です。各実行が自分の音速で CFL を取るので、遅い音も速い音も同じ歩数を使います。違うのは所要時間で、その比は音速比と同じです。

$K$ も出力にあります。 $\alpha = 0.5$ で $K = 0.499947$ です。二つの体積弾性率が無限に離れると、この式の $K$ は 0.5 になります。印刷値は 0.5 より $5.3\times10^{-5}$ 小さいです。水が空気より約 1.9 万倍硬いという文と合います。 $\alpha = 0.01$ の Wood 音速は 118.6 m/s で、純水 1624.9 m/s の約 1/14 です。 $\alpha = 0.0001$ では 956.6 m/s になり、1624.9 からすでに 668 m/s 落ちています。同じ位置の $K=0$ 列は 1624.3 m/s です。気泡が音速を削るのは $K$ があるときだけです。

最小の 23.65 m/s は $\alpha = 0.501$ にあります。0.5 の 23.6 m/s と実質同じです。体積弾性率の比が大きければ最小は 0.5 に付きます。 $\alpha = 0.99$ で $K=0$ は 540.9 m/s で、純空気 341.6 m/s より 199 m/s 速いです。水 1% の $\pi$ が空気側へ漏れた結果です。同じ $\alpha$ の Wood 値は 112.4 m/s で、純空気よりむしろ遅いです。

空気1%から99%まで、二本の曲線が分かれる場所#

左の列は空気がひとつまみでも敏感です。0.01% あるだけで 1625 m/s が 957 m/s になり、1% なら 119 m/s です。最小が $\alpha\approx0.5$ に付くのも、体積弾性率の差が大きいためです。

右の列は空気 99% で純空気より速くなります。水 1% の $\pi$ が空気側へ漏れた結果で、報告書の "5 Eq. Trans" が縮約モデルから外れる場所もここです。水滴が音波をほとんど跳ね返す話は音響インピーダンスと sharp 界面で別に見ました。今回の二本の曲線は、界面での反射の前に、一様な混合物がどの音速を持つかを問います。

報告書が実験と照らし合わせたもの

報告書は音速曲線で終わりません。空気–水衝撃波管で、二つのモデルの相密度と内部エネルギーを並べて描いています。さらにエポキシ–スピネル粉末混合物の衝撃Hugoniot、つまり衝撃波速度 $U_s$ と粒子速度 $U_p$ の関係を描いています。図には5方程式縮約モデル、7方程式モデル、実験点が一緒に載っています。

この比較が重要なのには理由があります。 $K\nabla\cdot u$ は非保存項なので、衝撃波の跳び条件はモデルだけでは一意に決まりません。数値粘性の形が衝撃波後の状態を変えることがあります。だからこそ、縮約モデルが7方程式モデルや実験にどれだけ近いかを図で確認しているのです。今回の計算のような弱い波では差は消えます。強い衝撃波では別途検証が必要です。

23 m/s付近でなければ、そのコードはKを捨てている#

一様な混合物 $\alpha=0.5$ に小さな圧力パルスを置き、速度を測ります。23 m/s 付近なら $K$ は生きています。745 m/s 付近、つまり表の $K=0$ の理論値なら、体積分率は運ばれているだけです。CFL に入れる音速も同じ検査に入ります。Wood 音速のコードに凍結音速を入れると時間刻みは数十倍保守的になり、 $K=0$ モデルに 23.6 m/s を入れると CFL を超えます。その剛性を音響と輸送のどちらへ置くかは、Peluchon の音響-輸送 IMEXが選んだ場所と重なります。

空気 0.1% だけでも Wood 音速は 365 m/s で、水 1625 m/s の 1/4 を下回ります。キャビテーションの初期や、気泡を吸い込むポンプでは、この数字が圧力波の到達時間を決めます。そのセルが 1600 m/s と出ているなら、モデルは $K$ のない側に寄っています。0.01% でもすでに 957 m/s まで落ちるので、純水の音速をそのまま使うと到達時間がずれます。

[論文レビュー] 水と空気を半分ずつ混ぜたら音は24 m/sで這った — Murrone–Guillard 5方程式モデルのK∇·u項