水面上的呼喊为何传不进水下 — 跨越扩散界面的声传播

压力透射系数T超过1，能量也不会变多。99.9%更接近阻抗比的平方，而不是比值本身。

在泳池边呼喊朋友，潜在水里的人几乎听不见。空气中约 99.9% 的声音从水面弹回。同样的物理把海洋声学、医学超声和水下噪声分开。顺着 Ballout 等人（2025），这里用一个 1D 声学求解器把界面上的反射、透射和折射再写一遍。

用声速比和密度比把 $z_2/z_1$ 拉开，中央界面上的反射脉冲变大，透射脉冲变小。

sound speed ratio c2/c1: 1.80density ratio ρ2/ρ1: 2.50interface width ε: 0.060

z1 = 1.00z2 = 4.50R = 0.636T = 1.636

把两个阻抗调成 1:1，脉冲就穿过界面，反射消失。像空气–水那样的跳变，左边只剩下反射。

Abbas Ballout 等，"Acoustic Propagation/Refraction Through Diffuse Interface Models"，arXiv:2504.01727（2025）。他们改了不可压 Navier–Stokes/Cahn–Hilliard 的弱可压缩（weak compressibility）项，让两相带上不同的声速，声波就沿着这些声速走。

声阻抗 — 介质对声音施加的阻力

声音穿过介质时感受到的阻力称为声阻抗（密度×声速）。

z = \rho\, c_s

其中 $\rho$ 是密度， $c_s$ 是声速。空气为 $z_\text{air} = 1.2 \times 343 \approx 410$ ，水为 $z_\text{water} = 1000 \times 1481 \approx 1.5\times10^6$ 。相差约 3600 倍。

阻抗衡量声音能多容易地载入某介质。两种介质的阻抗相近时，声波平滑越过；相差悬殊时，大部分被弹回。声音在水面几乎被全反射，正是因为这个跳变。

把这个跳变放回反射系数。上面写出的舍入值 $z_1\approx 410$ 、 $z_2\approx 1.5\times 10^6$ ，比值约 3600。 $R=(z_2-z_1)/(z_1+z_2)$ 贴在 0.999。它的平方也在 0.999 附近，和前面说的 99.9% 对得上。压力透射 $T=2z_2/(z_1+z_2)$ 就在 2 下面。超过 1 并不表示能量变多。水里的质点速度按 $1/z$ 变小。

反射与透射 — 由阻抗跳变决定

对于垂直入射的平面波，反射系数 $R$ 和透射系数 $T$ 仅由阻抗决定（论文式 37）。

R = \frac{z_2 - z_1}{z_1 + z_2}, \qquad T = \frac{2 z_2}{z_1 + z_2}

$z_1$ 、 $z_2$ 分别是入射侧与透射侧的阻抗。空气→水时 $z_2 \gg z_1$ ，故 $R \to 1$ ，几乎完全反射。压力透射系数 $T$ 可以超过 1。压力幅值增大并不违反能量守恒：透射波的质点速度变小，能量通量依然守恒。

上面的声管已经把这个跳变压过。声速比和密度比越大，透射脉冲越小，反射脉冲越大。1:1 时反射消失。同一列声波不从出口弹回来的另一个位置，是 NSCBC 出口边界。

13度之墙 — Snell 定律与临界角#

斜入射时透射波的方向会偏折。支配光折射的同一条 Snell 定律（入射角与透射角的正弦之比等于声速之比） 对声波同样成立（论文式 41）。

\frac{\sin\theta_i}{c_{s1}} = \frac{\sin\theta_t}{c_{s2}}

$\theta_i$ 是入射角， $\theta_t$ 是透射角。当透射侧声速更快时（ $c_{s2} > c_{s1}$ ，如空气到水），一旦入射角越过某阈值， $\sin\theta_t$ 就超过 1，透射波无法存在。这个阈值就是临界角（critical angle）。

\theta_c = \arcsin\!\left(\frac{c_{s1}}{c_{s2}}\right) = \arcsin\!\left(\frac{343}{1481}\right) \approx 13.4^\circ

$343/1481$ 约等于 0.232。它的反正弦是 13.4 度。空气中的声源若以比这个角度更斜的方向触及水面，声音进不了水。全反射。

移动入射角和声速比，透射射线会贴到水面，或者再立起来。

incidence angle θi: 10°sound speed ratio c2/c1: 4.32

critical angle θc = 13.38°wave transmits into medium 2

把入射角推向 13 度，透射射线几乎贴着水面后消失。那一刻只剩反射射线。降低声速比则临界角增大，更斜的入射也允许透射。

用弱可压缩把声波载上

不可压求解器仅把压力当作满足散度约束的工具。那样的压力场无法传播声波。论文引入了**弱可压缩（给压力加上时间导数项以赋予有限声速）**方程。

\partial_t p + \rho\, c_s^2\, \nabla\cdot\mathbf{u} = 0

$p$ 是压力， $\mathbf{u}$ 是速度， $c_s$ 是该点的声速。改动的位置是按相插值 $c_s$ 。用浓度 $c$ 混合密度与声速。

(\,\cdot\,) = (\,\cdot\,)_1\, c + (\,\cdot\,)_2\,(1 - c)

凭这一处修改，同一求解器在介质 1 中以 $c_{s1}$ 、在介质 2 中以 $c_{s2}$ 传播声波。作者用熵稳定的间断 Galerkin 谱元法（DGSEM）离散，并展示了谱收敛。

扩散界面带来的建模误差

与把界面切到单个网格的 VOF 或水平集不同，Cahn–Hilliard 模型将界面抹开到宽度 $\varepsilon$ 。

c_0 = 1 - 0.5\left(1 + \tanh\frac{2x}{\varepsilon}\right)

这种抹开会带来建模误差。论文用数值实验确认，该误差按界面宽度的平方衰减，精确地说是以波长归一化后的 $\mathcal{O}\!\left((\varepsilon/\lambda)^2\right)$ 。当 $\varepsilon \to 0$ 时，它收敛到锐界面的解析解。前面 1D 声管里的 "interface width" 滑块就是这个 $\varepsilon$ 。调大它，透射与反射脉冲会略微变模糊。

Python — 在 1D 管中再现系数#

线性声学方程用交错网格 FDTD 积分。发射一个高斯脉冲，把越过界面后的反射、透射幅值与解析值比较。

import numpy as np
 
def impedance(rho, c):
    """声阻抗 z = rho * c (密度 × 声速)"""
    return rho * c
 
def analytic_coeffs(z1, z2):
    """垂直入射平面波的反射 / 透射系数 (论文式37)"""
    R = (z2 - z1) / (z1 + z2)
    T = 2.0 * z2 / (z1 + z2)
    return R, T
 
def tanh_interface(x, eps):
    """把界面抹开到宽度 eps 的浓度场 (论文式32)"""
    return 1.0 - 0.5 * (1.0 + np.tanh(2.0 * x / eps))
 
def run_acoustic_tube(rho1=1.0, c1=1.0, rho2=2.5, c2=1.8,
                      N=2000, L=4.0, eps=0.04, steps=2600):
    """积分1D线性声学系统并测量反射、透射脉冲幅值。"""
    dx = L / N
    x = -L / 2 + (np.arange(N) + 0.5) * dx
    conc = tanh_interface(x, eps)            # 1: 介质1, 0: 介质2
    rho = rho1 * conc + rho2 * (1.0 - conc)  # 密度插值
    cs  = c1   * conc + c2   * (1.0 - conc)  # 声速插值
    K   = rho * cs**2                        # 体积模量 rho c^2
 
    p = np.exp(-((x + 0.55) / 0.12)**2)      # 右行高斯脉冲
    u = p / (rho1 * c1)                      # 介质1特征: u = p/(rho c)
 
    dt = 0.4 * dx / cs.max()                 # CFL限制的时间步长
    left, right = N // 4, 3 * N // 4         # 左 / 右探针
    p_inc = p.max()
    refl_peak = tran_peak = 0.0
    for n in range(steps):
        u[:-1] += -dt / (0.5 * (rho[:-1] + rho[1:]) * dx) * (p[1:] - p[:-1])
        p[1:]  += -dt * K[1:] / dx * (u[1:] - u[:-1])
        if n > steps // 2:                   # 脉冲越过界面之后
            refl_peak = max(refl_peak, np.abs(p[:left]).max())
            tran_peak = max(tran_peak, np.abs(p[right:]).max())
    return p_inc, refl_peak, tran_peak
 
if __name__ == "__main__":
    rho1, c1, rho2, c2 = 1.0, 1.0, 2.5, 1.8
    z1, z2 = impedance(rho1, c1), impedance(rho2, c2)
    R, T = analytic_coeffs(z1, z2)
    p_inc, refl, tran = run_acoustic_tube(rho1, c1, rho2, c2)
    print(f"z1={z1:.3f}  z2={z2:.3f}")
    print(f"解析解   |R|={abs(R):.3f}  T={T:.3f}")
    print(f"测量值   |R|={refl / p_inc:.3f}  T={tran / p_inc:.3f}")

在 $z_1=1.0$ 、 $z_2=4.5$ 时，测量值贴近解析的 $|R|=0.636$ 、 $T=1.636$ 。默认值是 $\rho_1=1$ 、 $c_1=1$ 、 $\rho_2=2.5$ 、 $c_2=1.8$ ，所以 $z_2=\rho_2 c_2=4.5$ 。 $R=(4.5-1)/(4.5+1)=3.5/5.5=0.636$ ， $T=2\times 4.5/5.5=1.636$ 。 $0.636$ 的平方在 0.405 附近。即便 $T=1.636$ ，这个平方也不会把反射能量推过 1。质点速度按阻抗比减小，通量仍然对齐。

网格是 $N=2000$ 、 $L=4$ ，所以 $dx=0.002$ 。 $\varepsilon=0.04$ 时，宽度参数是 20 个格子。这和 3D 空气–水要求界面内多于 10 个点是同一处。声速最大值是 $c_2=1.8$ ，时间步为 $dt=0.4\,dx/1.8$ 。CFL 0.4 盯着快的一侧。2600 步的一半，也就是 1300 步之后，左右探针才开始记峰值。加密网格并减小 $\varepsilon$ ，和解析解的差距会进一步缩小。

再现过程中令人生疑之处

首先撞上的陷阱是分离反射波。左探针看到的是入射波与反射波之和。论文另跑一次单相（single-phase）模拟来减去入射部分。上面的代码只测量脉冲离开后左侧区域的残余信号，从而绕开这点，但在入射与反射重叠的区间会对不齐。

弱可压缩并非真正的可压缩。它无法处理激波或大马赫数。作者本人也将其限定于低马赫、线性声学。它的目的与用可压缩求解器直接求解声学不同。

3D 空气–水算例用了 6750 万自由度，因为扩散界面里必须塞进十个以上的点才有精度。界面越薄，代价越往上冲。这种方法能否取代 VOF，取决于问题规模。

下一篇值得读的论文

把 $\varepsilon$ 减小，误差按 $(\varepsilon/\lambda)^2$ 下降。在 3D 里，界面内 10 个点的要求把自由度推到了 6750 万。1D 默认的 $\varepsilon/dx=20$ ，是那 10 个点的两倍厚。再切薄之前，先看 0.636 和 1.636 是否已经贴上。贴不上的话，问题在探针时刻，不在反扩散界面锐化。

同一团队的熵稳定 DGSEM iNS/CH 原始论文（Manzanero 等，2020）是这个弱可压缩项坐上去的离散。出口把声波吞掉的边界，和上面的 NSCBC 接在一起。在强行把界面切薄之前，先确认阻抗比为 3600 的垂直入射里， $R$ 是否贴在 0.999。