[论文评述] 液滴几乎把声波全部弹回 — 声阻抗与sharp界面两相求解器

腔体里的2026 Pa是Laplace压力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。一个液体单元就把整体时间步收紧四倍以上。

液体里面物质移动的速度大约是 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%。显式时间步被那个单元的声速拴住，所以只糊开一格，就会收紧十二倍。

驻波是2026 Pa对Laplace的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倍的扰动就进不了液滴。界面只混开一格，120.5 m/s的单元就会把21.7 kHz的节推走。

透射的0.111%再拆开看，压力透射系数接近2，水侧速度却被阻抗比压到近乎为零，能量因此进不去。显式步长认的是这一格里的最低声速，不是两侧材料自己的声速。

再看 $\alpha = 0.01$ 这一行。水里只混进 1% 的气体，声速就从 1500 掉到 120 m/s。这是含气水里测到的值，也是界面只糊开一格时格子会看见的值。

驻波里的液滴 — 论文最后给出的图

最后的试验把半径 1 mm 的液滴放进高度 $L = 8$ 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}

液滴待在这片振荡速度场的正中间。压力在液体内部也会传播，最初圆形液滴内外的压差是 Laplace 压力 $2\sigma/r = 116$ Pa，保持不变。声学扰动 2026 Pa 和表面张力守着的 116 Pa 之间的差距，决定结果。每个周期液滴沿轴向拉长再压扁，振幅够大就走向破碎。

这张图需要前面那套写法。把液体当成不可压，液滴内部的压力波就没了。把界面糊开，附近单元里 24 m/s 的假声速会晃动节的位置。sharp 界面，而且两相都可压，这张图才会出来。

只糊开一格就会顶走压力节的选择

$\alpha = 0.5$ 的 24 m/s，是界面一直糊着时求解器必然踩到的数。把这个声速留在方程里的五方程模型，是 Murrone–Guillard 的 Wood 声速，它故意带着同一个 24 m/s 走。

在液滴的驻波里，留下这 24 m/s，节还站得住吗？还是必须抹掉那一格 120.5 m/s，21.7 kHz 才会回到原位？