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

圧力透過係数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 を超えてもエネルギーが増えたわけではありません。水中の粒子速度はインピーダンスに反比例して小さくなります。

反射と透過 — インピーダンスの跳びが分ける

垂直入射する平面波の反射係数 $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 万自由度を使いました。拡散界面に 10 点以上を詰め込まないと精度が出ないためです。界面が薄くなるほどコストが爆発します。この手法が 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 に付くかを先に確認します。