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

在泳池边大声呼喊朋友,潜在水里的人却几乎听不见。空气中约99.9%的声音直接从水面弹回。同样的物理支配着海洋声学、医学超声和水下噪声预测。本文跟随 Ballout 等人(2025)的论文,用一个从零写起的1D声学求解器再现声波跨越界面时的反射、透射与折射。读完之后,仅凭一个阻抗跳变就能预测这一切。

先看论文信息。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倍。

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

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

对于垂直入射的平面波,反射系数 $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。压力幅值增大并不违反能量守恒:透射波的质点速度变小,能量通量依然守恒。

在下面的模拟中亲手调节参数吧。一个高斯脉冲从左侧出发,撞击中央的界面。

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

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

提高声速比和密度比以拉大 $z_2/z_1$ ,透射脉冲会变小,反射脉冲会变大。反之,把两个阻抗调成1:1,脉冲就直接穿过界面。反射消失了。

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

空气中的声源若以比13.4度更斜的角度触及水面,声音一丝也进不去水中。全反射。

在下面改变入射角和声速比。

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

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

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

用弱可压缩把声波载上#

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

\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$ 。加密网格并减小 $\varepsilon$ ,两者差距会进一步缩小。

再现过程中令人生疑之处#

首先撞上的陷阱是分离反射波。左探针看到的是入射波与反射波之和。论文另跑一次单相(single-phase)模拟来减去入射部分。上面的代码通过只测量脉冲离开后左侧区域的残余信号绕开了这点,但在入射与反射重叠的区间并不准确。

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

第三,3D空气-水算例用了6750万自由度,因为扩散界面里必须塞进十个以上的点才有精度。界面越薄,代价越爆炸。这种方法在实际中能否取代VOF,取决于问题规模。

下一篇值得读的论文#

只留要点。声波的命运取决于阻抗跳变 $z_2/z_1$ 。斜入射按Snell定律折射,超过 $\theta_c$ 则全反射。每相插值一个弱可压缩项,即便是不可压求解器也能以正确速度传播声波。

若想深入,推荐继续阅读同一团队的熵稳定DGSEM iNS/CH 原始论文(Manzanero 等,2020),或把 Ffowcs-Williams–Hawkings 声学类比扩展到多相流的研究。