Por qué un grito sobre el agua nunca llega abajo — Acústica a través de una interfaz difusa

Un coeficiente de transmisión de presión mayor que 1 no aumenta la energía. El 99,9% queda cerca del cuadrado del salto de impedancia.

Un grito desde el borde de una piscina apenas llega a quien nada debajo. Cerca del 99,9% del sonido en el aire rebota en la superficie. La misma física separa la acústica oceánica, el ultrasonido médico y el ruido submarino. Siguiendo a Ballout et al. (2025), estas notas rehacen la reflexión, la transmisión y la refracción con un solver acústico 1D.

Al abrir $z_2/z_1$ con las razones de velocidad y de densidad, el pulso reflejado en la interfaz central crece y el transmitido se encoge.

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

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

Con ambas impedancias en 1:1 el pulso cruza y la reflexión desaparece. Un salto aire–agua deja solo el pulso reflejado a la izquierda.

Abbas Ballout et al., "Acoustic Propagation/Refraction Through Diffuse Interface Models", arXiv:2504.01727 (2025). Ajustan el término de compresibilidad débil (weak compressibility) de un sistema incompresible de Navier–Stokes/Cahn–Hilliard para que dos fases tengan velocidades del sonido distintas, y la onda viaja sobre ellas.

Impedancia acústica — la resistencia que un medio opone al sonido#

La resistencia que un medio ofrece a una onda sonora que lo atraviesa es su impedancia acústica (densidad por velocidad del sonido).

z = \rho\, c_s

Aquí $\rho$ es la densidad y $c_s$ la velocidad del sonido. El aire da $z_\text{air} = 1.2 \times 343 \approx 410$ ; el agua da $z_\text{water} = 1000 \times 1481 \approx 1.5\times10^6$ . Un factor de 3600 de diferencia.

La impedancia mide con qué facilidad el sonido se acopla a un medio. Cuando dos impedancias son parecidas, la onda cruza con suavidad. Cuando difieren enormemente, la mayor parte rebota. Ese salto en la superficie del agua es la razón de la reflexión.

Ese salto vuelve a $R=(z_2-z_1)/(z_1+z_2)$ . Con $z_1\approx 410$ y $z_2\approx 1.5\times 10^6$ la razón es cerca de 3600, $R$ se pega en 0.999 y el cuadrado es el 99,9% del inicio. $T=2z_2/(z_1+z_2)$ queda justo bajo 2.

Reflexión y transmisión — lo decide el salto de impedancia#

Para una onda plana con incidencia normal, el coeficiente de reflexión $R$ y el de transmisión $T$ dependen solo de la impedancia (ec. 37 del artículo).

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

$z_1$ y $z_2$ son las impedancias del lado de incidencia y del lado de transmisión. Para aire-a-agua, $z_2 \gg z_1$ , así que $R \to 1$ : reflexión casi total. El coeficiente de transmisión de presión $T$ puede superar 1. Una amplitud de presión mayor no viola la conservación de energía: la velocidad de partícula transmitida disminuye, de modo que el flujo de energía sigue equilibrado.

El tubo de arriba ya empujó ese salto. Al crecer las razones de velocidad y de densidad, el pulso transmitido se encoge y el reflejado crece. En 1:1 la reflexión desaparece. Otro sitio donde la misma onda no debe volver desde una salida es una salida NSCBC.

El muro de los 13 grados — la ley de Snell y el ángulo crítico#

Con incidencia oblicua, la onda transmitida se desvía. La misma ley de Snell (la razón de los senos de los ángulos de incidencia y transmisión es igual a la razón de las velocidades del sonido) que rige la luz vale también para el sonido (ec. 41 del artículo).

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

$\theta_i$ es el ángulo de incidencia y $\theta_t$ el de transmisión. Cuando el medio de transmisión es más rápido ( $c_{s2} > c_{s1}$ , como del aire al agua), $\sin\theta_t$ supera 1 en cuanto el ángulo de incidencia pasa cierto umbral. No puede existir onda transmitida. Ese umbral es el ángulo crítico (critical angle).

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

$343/1481$ vale cerca de 0.232. El arcoseno de eso es 13,4 grados. Si una fuente en el aire incide sobre la superficie con más oblicuidad que ese ángulo, el sonido no entra al agua. Reflexión total.

Al mover el ángulo de incidencia y la razón de velocidades, el rayo transmitido se tumba sobre la superficie o vuelve a levantarse.

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

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

Al llevar el ángulo de incidencia hacia los 13 grados, el rayo transmitido casi se tumba sobre la superficie antes de desaparecer. En ese instante solo queda el rayo reflejado. Al reducir la razón de velocidades, el ángulo crítico crece y admite transmisión incluso con incidencia más pronunciada.

Cargar el sonido sobre un modelo de compresibilidad débil#

Un solver incompresible usa la presión solo como herramienta para imponer la restricción de divergencia. Ese campo de presión no puede propagar sonido. El artículo introduce una ecuación de compresibilidad débil (un término de derivada temporal añadido a la presión que otorga una velocidad del sonido finita).

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

$p$ es la presión, $\mathbf{u}$ la velocidad y $c_s$ la velocidad del sonido local. El sitio editado es la interpolación de $c_s$ por fase. La densidad y la velocidad del sonido se mezclan mediante la concentración $c$ .

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

Con este cambio, el mismo solver transporta el sonido a $c_{s1}$ en el medio 1 y a $c_{s2}$ en el medio 2. Los autores lo discretizan con un método de Galerkin discontinuo de elementos espectrales (DGSEM) estable en entropía y muestran convergencia espectral.

El error de modelado que introduce una interfaz difusa#

A diferencia de VOF o level set, que recortan la interfaz a una sola celda, el modelo de Cahn–Hilliard la difumina en un ancho $\varepsilon$ .

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

Esa difuminación introduce un error de modelado. El artículo confirma numéricamente que el error decae como el cuadrado del ancho de interfaz, exactamente como $\mathcal{O}\!\left((\varepsilon/\lambda)^2\right)$ una vez normalizado por la longitud de onda. Cuando $\varepsilon \to 0$ , converge a la solución analítica de interfaz nítida. El control "interface width" del tubo 1D es este $\varepsilon$ . Al subirlo, los pulsos transmitido y reflejado se difuminan un poco.

Python — reproducir los coeficientes en un tubo 1D#

Las ecuaciones acústicas lineales se integran con un esquema FDTD en malla escalonada. Se dispara un pulso gaussiano y se comparan las amplitudes reflejada y transmitida tras la interfaz con los valores analíticos.

import numpy as np
 
def impedance(rho, c):
    """Impedancia acustica z = rho * c (densidad x velocidad del sonido)"""
    return rho * c
 
def analytic_coeffs(z1, z2):
    """Coeficientes de reflexion / transmision, onda plana incidencia normal (ec. 37)"""
    R = (z2 - z1) / (z1 + z2)
    T = 2.0 * z2 / (z1 + z2)
    return R, T
 
def tanh_interface(x, eps):
    """Campo de concentracion que difumina la interfaz en un ancho eps (ec. 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):
    """Integra el sistema acustico lineal 1D y mide las amplitudes de los pulsos."""
    dx = L / N
    x = -L / 2 + (np.arange(N) + 0.5) * dx
    conc = tanh_interface(x, eps)            # 1: medio 1, 0: medio 2
    rho = rho1 * conc + rho2 * (1.0 - conc)  # interpolacion de densidad
    cs  = c1   * conc + c2   * (1.0 - conc)  # interpolacion de velocidad del sonido
    K   = rho * cs**2                        # modulo de compresibilidad rho c^2
 
    p = np.exp(-((x + 0.55) / 0.12)**2)      # pulso gaussiano hacia la derecha
    u = p / (rho1 * c1)                      # caracteristica del medio 1: u = p/(rho c)
 
    dt = 0.4 * dx / cs.max()                 # paso temporal limitado por CFL
    left, right = N // 4, 3 * N // 4         # sondas izquierda / derecha
    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:                   # despues de que el pulso cruza
            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"analitico  |R|={abs(R):.3f}  T={T:.3f}")
    print(f"medido     |R|={refl / p_inc:.3f}  T={tran / p_inc:.3f}")

Con $z_1=1.0$ y $z_2=4.5$ , lo medido sigue a $|R|=0.636$ y $T=1.636$ . Los valores por defecto $\rho_2=2.5$ , $c_2=1.8$ dan $z_2=4.5$ , y $R=3.5/5.5=0.636$ , $T=9/5.5=1.636$ . El cuadrado de 0.636 queda cerca de 0.405.

La malla es $N=2000$ , $L=4$ , $dx=0.002$ . $\varepsilon=0.04$ ocupa 20 celdas, por encima de los 10 puntos del caso 3D. Con $c_2=1.8$ , $dt=0.4\,dx/1.8$ . Los picos se registran tras el paso 1300 de 2600.

Los valores por defecto también fijan $\rho_1=1$ y $c_1=1$ , así que $z_1=1$ . El pulso sale en $x=-0.55$ con ancho 0.12, y las sondas quedan en $N/4$ y $3N/4$ . Medir solo después del paso 1300 separa la onda que ya se fue de la que vuelve, aunque el solape temprano sigue ensuciando el pico. Una malla más fina y un $\varepsilon$ menor acercan 0.636 y 1.636 al analítico.

Lo que generó dudas al reproducirlo#

La primera trampa fue separar la onda reflejada. La sonda izquierda ve la suma de la presión incidente y la reflejada. El artículo corre una simulación monofásica (single-phase) aparte para restar la parte incidente. El código de arriba mide solo la señal residual a la izquierda tras la marcha del pulso, pero se desvía mientras las ondas incidente y reflejada se solapan.

La compresibilidad débil no es compresibilidad real. No sirve para ondas de choque ni números de Mach altos. Los propios autores la restringen a acústica lineal de bajo Mach. Su propósito difiere de resolver la acústica directamente con un solver compresible.

El caso 3D aire–agua consumió 67,5 millones de grados de libertad, porque la interfaz difusa necesita más de diez puntos comprimidos dentro. El costo se dispara al adelgazar la interfaz. Que este método reemplace a VOF depende de la escala del problema.

Un artículo para leer a continuación#

Al reducir $\varepsilon$ , el error baja como $(\varepsilon/\lambda)^2$ , pero en 3D los 10 puntos dentro de la interfaz ya costaron 67,5 millones de grados de libertad. El defecto 1D es $\varepsilon/dx=20$ . Si 0.636 y 1.636 no coinciden, falla el instante de la sonda, no el afilado por antidifusión.

Manzanero et al. (2020) es la discretización DGSEM sobre la que se monta el término. Antes de adelgazar la interfaz, la incidencia normal con razón 3600 tiene que dejar $R$ en 0.999.