[论文评述] 液滴几乎把声波全部弹回 — 声阻抗与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。仅仅一个液体网格就把全局时间步压紧四倍以上。

而液体内部物质的运动速度只有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$ )。

值得注意的是密度不在这份清单上。密度由各相自己的状态方程决定。界面附近不对密度做平均，因此不会产生任何一种物质都不具备的中间状态。下面要看的虚假声速问题，就在这里被切断。

封闭方程用的是三次状态方程。

\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$ — 界面处压力几乎翻倍。

在下面的模拟中亲自调节看看。

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，黄色的压力波几乎原样返回，而在界面处振幅跃升一倍。青色的速度波则符号翻转，水侧振幅跌到 $10^{-4}$ 量级。压力翻倍、速度归零 — 界面表现得像一堵墙。接着把 $\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条件收紧了。

用原始变量求解要付的代价

论文的表述中最具争议的部分，是选择用原始变量 $(p, \vec{u}, T)$ 而非守恒变量 $(\rho, \rho\vec{u}, \rho E)$ 求解，并把能量方程改写成压力方程。

收获有两点。第一，状态方程给出的压力与半隐式压力修正算得的压力在每一步都一致。在守恒变量表述中，这两者会相互偏离是已知问题。第二，导热项可以隐式处理。对自然对流这类以热为主角的低马赫数问题，这是很大的差别。

失去的同样明确：总能量守恒变差。因此该求解器自行把适用范围限定在亚声速。若打算捕捉激波，应当换一种表述。

压力的正定性(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)中实测到的数值，同时也说明界面只要弥散一个网格，网格看到的会是什么。

驻波中的液滴 — 论文最后展示的算例

论文的最后一个算例，是在高 $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%的能量被反射，界面压力翻倍，液体侧速度实际上为零。
把界面弥散到若干网格上，Wood声速会降到24 m/s。这个值不属于任何一种物质，同时破坏波的到达时刻与显式CFL限制。sharp界面从源头上杜绝了它。
对流显式、声学用投影法隐式 — 这一分裂让 $c$ 从时间步中消失。代价是原始变量表述下总能量守恒的下降，因此适用范围止于亚声速。