水上で叫んだ声が水中に届かない理由 — 拡散界面を越える音響伝播

インピーダンスの跳びが音波を反射・屈折させる仕組みを弱圧縮性モデルで再現

プールの縁から友人を大声で呼んでも、水に潜った人にはほとんど聞こえません。空気中の音の約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)シミュレーションを別に走らせて入射波を引きます。上のコードはパルスが十分右へ去った後の左側領域の残留信号だけを測ることでこれを回避しましたが、入射と反射が重なる区間では不正確です。

第二に、弱圧縮性は本当の圧縮性ではありません。衝撃波や大きなマッハ数には適用できません。著者ら自身も低マッハ・線形音響に限定しています。OpenFOAMの圧縮性ソルバーで直接音響を解くのとは目的が異なります。

第三に、3D空気-水の例は6750万自由度を使いました。拡散界面に10点以上を詰め込まないと精度が出ないためです。界面が薄くなるほどコストが爆発します。実務でこの手法がVOFを置き換えるかは問題規模次第です。

次に読む論文#

要点だけ残します。音波の運命はインピーダンスの跳び $z_2/z_1$ が分けます。斜め入射はSnellの法則で屈折し、 $\theta_c$ を超えると全反射します。弱圧縮性項を一つ相ごとに補間すれば、非圧縮ソルバーでも音波を正しい速度で運べます。

この流れをさらに掘るなら、同じグループのエントロピー安定DGSEM iNS/CH原論文(Manzanero 他、2020)や、Ffowcs-Williams–Hawkings音響アナロジーを多相に拡張した研究を続けて読むことをお勧めします。