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

気液界面が音波の99.9%を反射する理由と、にじんだ界面が生む偽の音速

潜っている人にはプールサイドからの呼び声が届きません。水面を通り抜ける音響エネルギーは0.1%ほどだからです。残りはすべて空へ戻ります。

同じ数字が圧縮性二相流ソルバの中にも現れます。空気中に水滴を一つ浮かべて音波を当てると、計算格子の上でも波の99.9%は界面で引き返します。このきわめて大きな対比を数値的に扱うことが、液体を含む圧縮性解析の本当の難所です。

論文: 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$ を厚さゼロの面として扱い、その上で満たすべき条件を直接課すという意味です。論文が使うのは次の三つです。

[\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$ )。

注目すべきは密度がこのリストにないことです。密度は各相の状態方程式がそれぞれ決めます。界面付近で密度を平均しないため、どちらの物質にも属さない中間状態が生まれません。次に見る偽の音速の問題が、ここで断たれます。

閉じるための状態方程式には3次のEoSを使います。

\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$ — 界面で圧力はほぼ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.

観察のポイントは二つです。medium 2をwaterにすると、黄色の圧力波はほぼそのまま戻り、界面では振幅が2倍に跳ね上がります。一方、水色の速度波は符号が反転し、水側の振幅は $10^{-4}$ 程度まで落ちます。圧力は2倍、速度はゼロ — 界面が壁のように振る舞っています。次に $\rho_2$ を下げながら $c_2$ を同じ比率で上げてみると、絵はまったく変わりません。積しか見ていないからです。

にじんだ界面は存在しない音速を作る

今度は界面を格子数セルにわたってにじませたらどうなるかを見ます。拡散界面(diffuse interface)系の手法が必然的に置かれる状況です。

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

\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です。問題に登場するどの物質も、その速さで信号を伝えません。この数字は物理ではなく平均化の副産物です。

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セルに上げてみてください。8セルは、どの捕捉法でも数百ステップで到達する厚みです。水色の曲線がどちらの物質にもない谷を掘り、下のレースではにじんだ信号がその区間でつまずきます。 $\rho_\ell$ を下げると谷は埋まります — この罠を作るのは音速の対比ではなく密度の対比です。

副作用は時間刻みにも及びます。陽解法なら、その谷のセルで $\Delta t \le \Delta x / c_m$ を守らなければなりません。にじみがCFLを締めつけます。

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

論文の定式化で議論を呼ぶのは、保存変数 $(\rho, \rho\vec{u}, \rho E)$ ではなく原始変数 $(p, \vec{u}, T)$ で解くという判断です。エネルギー方程式を圧力方程式に書き換えます。

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

失うものも明確です。全エネルギー保存が悪化します。そのためこのソルバは適用範囲を亜音速に自ら限定しています。衝撃波を捉えるつもりなら、別の定式化を使うべきです。

圧力の正値性(positivity)も保証しません。著者らはこれを欠陥ではなく余地と見ています。キャビテーション時の準安定液体は実際に負の圧力を持ちます。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

$\alpha = 0.01$ の行を見てください。水に気泡が1%混ざるだけで、音速が1500から120 m/sまで落ちます。実際の気泡水(bubbly water)で観測される値であり、同時に、界面が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です。重力はゼロにします。

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

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

液滴はこの振動速度場の真ん中に置かれます。圧力は液体の中も伝播し、初期の円形液滴の内外圧力差はLaplace圧 $2\sigma/r = 116$ Paで一定です。この二つの数字の大きさの差が結果を支配します。音響擾乱は2026 Pa、表面張力が守ろうとする圧力差は116 Paです。一周期ごとに液滴は軸方向に伸びては潰れ、振幅が大きくなれば分裂に至ります。

この計算がなぜ先の定式化を必要とするのかに注目してください。液体を非圧縮性にすると液滴内部の圧力波が消えます。界面をにじませると、界面近傍のセルで24 m/sの偽の音速が定在波の節の位置を揺らします。sharp界面かつ両相とも圧縮性 — この組み合わせでこそ、この絵が出てきます。

覚えておくこと三行

界面での音波の運命は $Z = \rho c$ 一つで決まります。空気-水は $Z$ 比が3600なので、エネルギーの99.9%が反射し、界面圧力は2倍に跳ね、液体側の速度は事実上ゼロです。
界面を格子数セルにわたってにじませるとWood音速が24 m/sまで下がります。どちらの物質にもない値であり、波の到達時刻と陽的CFLを同時に壊します。sharp界面はこれを根元から断ちます。
移流は陽的、音響は射影法で陰的 — この分離のおかげで時間刻みから $c$ が抜けます。代償は原始変数定式化での全エネルギー保存の低下であり、だからこそ適用範囲は亜音速までです。