[論文レビュー] 液滴は音波をほとんど跳ね返す — 音響インピーダンスとsharp界面二相ソルバ

空洞の2026 Paはラプラス116 Paの17.5倍なので、液滴は一周期のあいだ円を保てません。

ダイバーは水面の上から呼ばれても、ほとんど聞こえません。水へ入る音響エネルギーは0.1%ほどで、残りは空へ戻ります。空気中の液滴に音波を当てると、格子の上でも99.9%が界面で引き返します。液体を含む圧縮性計算の難しさは、この対比です。

medium 2をwaterにすると、黄色い圧力波はほぼ戻り、界面で振幅が二倍になります。 $\rho_2$ を下げる分だけ $c_2$ を上げると、絵はそのままです。

medium 1

medium 2 →

ρ₂998 kg/m³c₂1500 m/sinterface

What to watch: set medium 1 = air, medium 2 = water. The yellow pressure pulse comes back almost intact and nearly doubles at the interface (T = 1.998), while the cyan velocity pulse flips sign and the liquid barely moves. Now press medium 2 = air — R collapses to 0 and the pulse walks straight through. Only the ratio Z₂/Z₁ matters: slide ρ₂ down and c₂ up so that the product stays put, and the picture does not change at all.

水色の速度波は符号が反転し、水側の振幅は $10^{-4}$ 付近で死にます。圧力は二倍、速度はほぼ0です。界面が壁のように振る舞う場面です。

論文: A. Urbano, M. Bibal, S. Tanguy, A semi implicit compressible solver for two-phase flows of real fluids, Journal of Computational Physics (2022). ISAE-SUPAERO / IMFT。Level set + ghost fluidのsharp界面、原始変数の半陰的投影法、van der Waals状態方程式。

液体が入ると音響の時間刻みが先に死ぬ

圧縮性二相流をすべて陽的に解くと、時間刻みは音響CFLに縛られます。

\Delta t \le \frac{\Delta x}{\max(|u| + c)}

$c$ は音速、 $\Delta x$ は格子間隔です。 $c$ は相ごとに違います。空気は347 m/s、水は1500 m/sです。液体セル一つが、全体の時間刻みを4倍以上締めます。

液体の中で物質が動く速さは1 m/sほどです。マッハ数（流速/音速）では $10^{-3}$ です。見たいのは界面の変形なのに、時間刻みを決めるのは誰も追っていない音波です。

以前の回り道は、液体を非圧縮にする混合でした。時間刻みは解けますが、液体の中の音波が消えます。液滴が音にどう応えるかを見ようとした目的も一緒に消えます。

この論文は方程式を対流と音響に分けます。対流は陽的、音響は投影法で陰的です。残る条件は次だけです。

\frac{1}{\Delta t} = \frac{1}{\Delta t_{\mathrm{conv}}} + \frac{1}{\Delta t_{\mu}} + \frac{1}{\Delta t_{\sigma}}

$\Delta t_{\mathrm{conv}} = \Delta x / \max\lVert \vec{u} \rVert$ は対流、 $\Delta t_{\mu} = \rho \Delta x^{2} / 2\mu$ は粘性、 $\Delta t_{\sigma} = \tfrac{1}{2}\sqrt{\max(\rho)\Delta x^{3}/\sigma}$ は表面張力です。 $c$ はどこにもありません。

界面で何がつながり、何が切れるか

sharp界面は $\Gamma$ を厚さ0の面として、その上で守る条件を直接課します。論文の三つの条件は次です。

[\vec{u}]_{\Gamma} = 0, \qquad [p]_{\Gamma} = \sigma \kappa + 2\left[\mu \frac{\partial u_n}{\partial n}\right]_{\Gamma} - \frac{2}{3}\left[\mu \nabla \cdot \vec{u}\right]_{\Gamma}, \qquad [-k \nabla T \cdot \vec{n}]_{\Gamma} = 0

$[\,\cdot\,]_{\Gamma}$ は界面をまたぐジャンプ、 $\sigma$ は表面張力、 $\kappa$ は曲率、 $\mu$ は粘性係数、 $k$ は熱伝導率です。速度はつながり、圧力は表面張力の分だけ切れ、熱流束はつながります。質量移動はないと置くので、温度も連続です（ $[T]_{\Gamma}=0$ ）。

リストに密度がありません。密度は各相の状態方程式が別に決めます。界面付近で密度を平均しないので、二つの相のあいだにない中間状態は生まれません。あとの偽の音速はここで止まります。

状態方程式は三次です。

\rho^{3} + a_{1}\rho^{2} + a_{2}\rho + a_{3} = 0

係数 $a_1, a_2, a_3$ は圧力・温度・物性で決まります。van der Waals、Peng-Robinson、Redlich-Soave-Kwongがすべてこの形です。気体と液体を一つの式が一緒に持ちます。完全気体は液体を扱えず、Tait式は気体を扱えません。

インピーダンス一つが反射率を全部決める

界面で波が分かれる仕方は、音響インピーダンス一つで終わります。

Z = \rho c

媒質1から媒質2へ入るとき、圧力の反射・透過係数は次です。

R = \frac{Z_2 - Z_1}{Z_2 + Z_1}, \qquad T = 1 + R = \frac{2 Z_2}{Z_2 + Z_1}

$R$ は反射波の圧力振幅比、 $T$ は透過波の圧力振幅比です。密度だけでも音速だけでもありません。積だけを見ます。

空気と水を入れると $Z_{\text{air}} = 416$ 、 $Z_{\text{water}} = 1.50 \times 10^{6}$ 、比は3600です。 $R = 0.9994$ 、エネルギーでは99.89%が戻ります。 $T = 1.999$ なので界面の圧力はほぼ二倍です。前の図で $\rho_2$ と $c_2$ を逆に動かしても変わらなかった理由が、この積です。

界面が送り返した波は、計算領域の出口でもう一度会います。出口がそれを飲み込むか、また跳ね返すかはNSCBCの $\sigma$ が決めます。

ぼかした界面は存在しない音速を作る

界面を格子の数セルにまたがってぼかすと、拡散界面の手法がいつも出会う状況になります。

一つのセルに気体と液体が混ざると、密度と圧縮率はこう平均されます。

\rho_m = \alpha \rho_g + (1-\alpha)\rho_\ell, \qquad \frac{1}{\rho_m c_m^{2}} = \frac{\alpha}{\rho_g c_g^{2}} + \frac{1-\alpha}{\rho_\ell c_\ell^{2}}

$\alpha$ は気体の体積分率、下付きの $g$ / $\ell$ は気体/液体です。Woodの関係式です。混合物は気体のように柔らかく、液体のように重い。

$\alpha = 0.5$ なら $c_m = 24$ m/sです。空気の1/14、水の1/60です。問題の中のどの物質も、その速さで信号を送りません。平均が作った数字です。ぼかし幅を手で広げる競走は、出力表の直後に置きます。

原始変数で解く選択の代償

論文で論争になる部分は、保存変数 $(\rho, \rho\vec{u}, \rho E)$ ではなく原始変数 $(p, \vec{u}, T)$ で解く決定です。エネルギー方程式を圧力方程式に書き換えています。

得るものは、二つの圧力の一致と、熱伝導の陰的な扱いです。状態方程式が返した圧力と、半陰的補正が計算した圧力が毎ステップ一致します。保存変数で解くと、この二つがずれることが知られています。熱伝導項も陰的に渡せるので、自然対流のように熱が主役の低マッハ数では差が大きいです。

失うものは全エネルギーの保存です。だからソルバは亜音速へ範囲を自分で切ります。衝撃波を捕らえるつもりなら、別の定式化が要ります。

圧力の正値性も保証しません。著者たちはこれを欠陥より余地と見ます。キャビテーションの準安定な液体は実際に負の圧力を持ち、van der Waalsのような状態方程式は負の圧力に正の温度を対応させられます。

Pythonで確かめる反射率とWoodの音速#

前の図に入った数字を打ち出します。

# 空気-水界面の反射・透過、そしてにじんだ界面のWood音速
from math import sqrt
 
AIR   = (1.2,   347.0)    # (rho, c)
WATER = (998.0, 1500.0)
 
def impedance_split(m1, m2):
    z1, z2 = m1[0] * m1[1], m2[0] * m2[1]
    R = (z2 - z1) / (z2 + z1)
    return z1, z2, R, 1.0 + R, R * R
 
def wood_speed(alpha, gas=AIR, liq=WATER):
    rho_m = alpha * gas[0] + (1.0 - alpha) * liq[0]
    inv   = alpha / (gas[0] * gas[1] ** 2) + (1.0 - alpha) / (liq[0] * liq[1] ** 2)
    return 1.0 / sqrt(rho_m * inv)
 
z1, z2, R, T, E = impedance_split(AIR, WATER)
print(f"Z_air = {z1:8.1f}   Z_water = {z2:10.1f}   Z2/Z1 = {z2/z1:6.0f}")
print(f"R = {R:.5f}   T = 1+R = {T:.5f}   reflected energy = {100*E:.3f} %")
print(f"velocity amp in water / incident = {(z1/z2)*T:.3e}")
 
print("\nalpha   c_wood [m/s]")
for a in (0.0, 0.01, 0.1, 0.5, 0.9, 0.99, 1.0):
    print(f"{a:5.2f}   {wood_speed(a):8.1f}")

Z_air =    416.4   Z_water =  1497000.0   Z2/Z1 =   3595
R = 0.99944   T = 1+R = 1.99944   reflected energy = 99.889 %
velocity amp in water / incident = 5.562e-04
 
alpha   c_wood [m/s]
 0.00     1500.0
 0.01      120.5
 0.10       40.1
 0.50       24.0
 0.90       39.9
 0.99      114.3
 1.00      347.0

smear n8 cellsρ_liquid998 kg/m³

What to watch: at n = 0 the two lanes are the same run. Pull n up to 8 cells — a width any capturing scheme reaches within a few hundred steps — and the cyan curve digs a hole down to ~24 m/s that neither material has, so the smeared signal falls behind and arrives late. Drop ρ_liquid and the hole fills in: the dip is driven by the density ratio, not by the sound speeds. The yellow line is the bill an explicit solver pays for that hole.

nスライダーを0から8セルへ上げると、水色の曲線がどちらの物質にもない谷を掘ります。 $\rho_\ell$ を下げると谷は埋まります。出力の $\alpha = 0.01$ の行、120.5 m/sが、その谷の入口です。密度比を削ると、音速比をそのままでも谷は消えます。

反射エネルギーは $R^{2}$ です。出力の $R = 0.99944$ を二乗すると99.889%で、印刷された百分率と一致します。通り抜ける側は0.111%で、導入の0.1%と同じ場所です。音速比1500/347は4.32なのに、インピーダンス比は3595です。密度比998/1.2に音速比を掛けると、同じ3595になります。その壁を立てているのは音速ではなく密度です。

$\alpha = 0.01$ で、すでに1500が120.5へ落ちます。12.4倍です。 $\alpha = 0.5$ の24 m/sまでは、そこからさらに約5倍です。被害の大部分は最初の1%にあります。陽的な時間刻みはそのセルの音速に縛られるので、1セルぼかすだけで十二倍締まります。

定在波は2026 Pa対ラプラスの116 Pa、比は17.5です。速度振幅4.9 m/sは空気の音速347 m/sの1.4%で、マッハ0.014です。先に書いた液体の対流1 m/sの五倍近くです。 $f = 347/(2 \times 0.008) = 21688$ Hzで、書かれた21.7 kHzと一致します。液体を非圧縮にすると、この17.5倍の擾乱は液滴の中へ入りません。界面を1セルだけ混ぜると、120.5 m/sのセルが21.7 kHzの節の位置を押します。

$\alpha = 0.01$ の行を見直すと、水に気体が1%混ざるだけで音速は1500から120 m/sへ落ちます。気泡水で観測される値であり、界面を1セルぼかしただけで格子が見る値でもあります。

定在波の中の液滴 — 論文が最後に見せたもの

最後の試験は、高さ $L = 8$ mmの円筒空洞に半径1 mmの液滴を置き、上の境界に圧力振動を与えます。

p(t) = p_0 + \Delta p \sin(\omega t)

$\omega = 2\pi f$ で、 $f = c_0/(2L) \approx 21.7$ kHzは空洞の最初の音響モードです。 $\Delta p$ は $p_0$ の2%、2026 Paです。重力は0です。

定在波が立つと $z = L/2$ に圧力の節（node）ができます。圧力の節は速度の腹（antinode）です。その速度振幅はインピーダンスから見積もられます。

\Delta v = \frac{\Delta p}{\rho_0 c_0} \approx \frac{2026}{416} \approx 4.9\ \text{m/s}

液滴はこの振動速度場の真ん中にあります。圧力は液体の中でも伝わり、最初の円形液滴の内外の圧力差はラプラス圧 $2\sigma/r = 116$ Paで一定です。音響の擾乱2026 Paと、表面張力が守る116 Paの差が結果を決めます。周期ごとに液滴は軸方向へ伸びて押し戻され、振幅が大きくなると分裂へ向かいます。

この絵が前の定式化を必要とした理由です。液体を非圧縮にすると、液滴内部の圧力波が消えます。界面をぼかすと、近くのセルの24 m/sという偽の音速が節の位置を揺らします。sharp界面で両方の相が圧縮性であって、この絵が出ます。

1セルのぼかしが節を押し出す選択#

$\alpha = 0.5$ の24 m/sは、界面をぼかしたままにするとソルバが踏む数字です。その音速を式の中に残す5方程式の側は、Murrone–GuillardのWood音速が同じ24 m/sをわざと連れていきます。

液滴の定在波で、その24 m/sを残すと節は持つか。それとも1セルの120.5 m/sを消して、21.7 kHzが元の位置に立つか。