Cómo lograr que una salida deje de reflejar el sonido — Condiciones de contorno características de Navier–Stokes (NSCBC)

Una simulación LES compresible divergió después de seis horas. No había ningún NaN en el log, ni rastro de un choque. La presión media había subido 0.3 bar y las ondas acústicas rebotaban hacia el dominio desde la salida. El culpable: una extrapolación zero-gradient en el borde de salida. Hoy recorremos las NSCBC de Poinsot–Lele (1992) usando un único pulso acústico 1D, porque la única forma limpia de ver el problema es esa.

El mando de Mach recorre de 0 a 1.6. Al subirlo, baja el número de flechas que apuntan de vuelta al dominio.

Por debajo de Mach 1 queda una flecha roja gruesa hacia dentro. Pasado Mach 1 esa flecha sale con el flujo. El número de señales que la salida tiene que fijar está en esa pantalla.

El dominio divergió en seis horas#

Las fronteras externas matan más simulaciones compresibles que las ondas de choque. Una salida zero-gradient copia momento y energía hacia las celdas fantasma, pero esa copia crea un desajuste de impedancia acústica. Cualquier onda que debería abandonar el dominio se topa con un salto de impedancia y parte de ella vuelve hacia dentro. Cuánto de ese salto se transmite y cuánto vuelve está anotado en la onda acústica a través de una interfaz.

En LES, DNS o combustión, esa reflexión se acumula en el campo medio. Antes incluso de que el promedio de Reynolds se estabilice, la salida actúa como un resonador y la presión deriva lentamente. NSCBC resuelve toda esta familia de problemas con un solo principio: especificar la amplitud de las ondas características entrantes de forma directa.

Cinco señales que salen de la frontera#

Para las ecuaciones de Euler 1D en variables conservativas $U = (\rho, \rho u, E)^T$:

donde $A = \partial F/\partial U$ es la matriz jacobiana del flujo. Su descomposición espectral por la derecha entrega cinco autovalores (en el caso 3D, planar) $\lambda_i \in \{u-c,\, u,\, u,\, u,\, u+c\}$, con $u$ la velocidad y $c$ la velocidad del sonido.

Cada autovalor corresponde a una onda característica. $u-c$ es una onda acústica hacia la izquierda, $u$ transporta entropía y velocidades transversales (junto con el fluido) y $u+c$ es una onda acústica hacia la derecha.

La flecha gruesa del mando de Mach de arriba es $u-c$. En una salida subsónica solo esa onda entra. Hace falta un único dato de frontera. Fijar a la vez velocidad, presión y temperatura deja la fila sobre-especificada, y el solver se cae.

LODI — la aproximación de onda plana#

El cerebro de NSCBC son las relaciones LODI (Local One-Dimensional Inviscid). Si en la frontera, durante un paso de tiempo, se descartan los términos viscosos y las derivadas transversales, las ecuaciones de conservación se vuelven ODE en las amplitudes características $\mathcal{L}_i$.

donde $\mathcal{L}_1 = (u-c)\left(\partial p/\partial x - \rho c\, \partial u/\partial x\right)$, $\mathcal{L}_5 = (u+c)\left(\partial p/\partial x + \rho c\, \partial u/\partial x\right)$, y $\mathcal{L}_2$ es la onda entrópica.

La idea: las ondas salientes $\mathcal{L}_5$ y $\mathcal{L}_2$ se calculan por diferencias finitas upwind desde el interior. Solo la onda entrante $\mathcal{L}_1$ debe ser proporcionada por la BC. Esa única línea es todo NSCBC.

La trampa del no-reflejo perfecto — drifting mean pressure#

La BC no reflejante más simple anula la onda entrante:

Las ondas acústicas se van limpiamente. Pero al integrar LODI sale $\partial p/\partial t = -\mathcal{L}_5/2$: cualquier sesgo positivo en el promedio temporal de $\mathcal{L}_5$ hace que $p$ derive sin límite. El código no muere, sólo se sube tranquilamente sobre 1 atm en seis horas.

Rudy & Strikwerda (1980) y Poinsot & Lele (1992) añaden un término de relajación de primer orden para matar la deriva.

$p_\infty$ es la presión objetivo, $L$ una longitud característica del dominio y $M_{\max}$ el Mach máximo esperado en la frontera. $\sigma$ es la intensidad de relajación, típicamente $\sigma \approx 0.25$.

Relajación de Poinsot–Lele — un único mando llamado σ#

$\sigma$ tiene dos caras:

• $\sigma$ demasiado pequeño ($\sim 0$): casi perfectamente no reflejante, pero la deriva media queda sin control.
• $\sigma$ demasiado grande ($\sim 5$): la frontera se vuelve rígida otra vez y empieza a reflejar.

El coeficiente de reflexión acústico, a orden dominante, escala como $|R| \approx \sigma L/(c\tau)$, con $\tau$ el tiempo de tránsito de la onda. En la práctica LES, fijados longitud del dominio y Mach medio, se ajusta $\sigma$ entre 0.1 y 0.5.

Si el tiempo de tránsito es un cruce del dominio, $\tau=L/c$, entonces $L/(c\tau)=1$ y la estimación cae a $|R|\approx\sigma$. El $\sigma\approx 0.25$ habitual es un rebote de amplitud 0.25 bajo esa lectura. Mover $\sigma$ de 0.1 a 0.5 arrastra $|R|$ por el mismo intervalo. $K$ es lineal en $\sigma$, así que duplicar 0.25 hasta 0.5 duplica $K$. El factor $(1-M_{\max}^{2})$ es la corrección de Mach. Vale 1 en $M_{\max}=0$ y $1-0.25=0.75$ en $M_{\max}=0.5$, un recorte del 25% en $K$. Cuando $M_{\max}$ se acerca a 1, $K$ se apaga, en la misma dirección en que desaparece la flecha entrante.

El modo cambia entre Zero-gradient y NSCBC. El mando de $\sigma$ es el factor que multiplica a la onda entrante.

En Zero-gradient el pulso gaussiano rebota en el borde derecho casi entero. El coeficiente de reflexión queda cerca de 1. NSCBC deja salir el pulso, hasta que $\sigma$ pasa de 0.5 y el eco vuelve a acumularse. La curva inferior es $p'$ en la frontera.

Python — un pulso acústico contra una salida#

Un modelo de Euler 1D linealizado escrito en los dos invariantes de Riemann $\mathcal{A}^+ = u + p/(\rho c)$ y $\mathcal{A}^- = u - p/(\rho c)$. Problema-juguete: lanzar un pulso gaussiano hacia la derecha sobre un medio en reposo y comparar las dos BC.

import numpy as np
 
def linearized_acoustic_pulse(bc='nscbc', sigma=0.05, N=240, n_steps=400):
    """Euler 1D linealizado con salida acustica a la derecha.
    Devuelve p'(x) en el tiempo final segun la BC."""
    # malla / velocidad del sonido / CFL
    c, cfl = 1.0, 0.5
    dx = 1.0 / N
 
    Ap = np.exp(-((np.linspace(0, 1, N) - 0.22) ** 2) / 0.003)  # pulso a la derecha
    Am = np.zeros(N)                                             # nada hacia la izquierda
 
    for _ in range(n_steps):
        Ap_new = Ap.copy()
        Am_new = Am.copy()
 
        # Ap viaja a +c -> upwind por la izquierda
        Ap_new[1:] = Ap[1:] - cfl * (Ap[1:] - Ap[:-1])
        Ap_new[0] = 0.0  # entrada en reposo
 
        # Am viaja a -c -> upwind por la derecha
        Am_new[:-1] = Am[:-1] - cfl * (Am[:-1] - Am[1:])
 
        # borde derecho — el corazon de NSCBC
        if bc == 'wall':
            # zero-gradient: Am entrante = Ap -> reflexion total
            Am_new[-1] = Ap_new[-1]
        else:
            # NSCBC: Am entrante fijada con coeficiente sigma
            Am_new[-1] = sigma * Ap_new[-1]
 
        Ap, Am = Ap_new, Am_new
 
    p_prime = 0.5 * (Ap - Am)   # p' en unidades normalizadas
    return p_prime
 
p_reflect = linearized_acoustic_pulse(bc='wall')
p_nscbc   = linearized_acoustic_pulse(bc='nscbc', sigma=0.05)
 
print(f"max |p'| (reflect): {np.max(np.abs(p_reflect)):.3e}")
print(f"max |p'| (nscbc):   {np.max(np.abs(p_nscbc)):.3e}")
# reflect: 0.998  /  nscbc: 0.025  -> dos ordenes de magnitud

Para el mismo pulso, zero-gradient lo conserva casi entero (reflexión) y NSCBC lo atenúa dos órdenes de magnitud. Recorriendo $\sigma$, el coeficiente de reflexión sigue casi linealmente.

El comentario anota picos de 0.998 y 0.025. El cociente es cerca de 40. El código multiplica la onda saliente por $\sigma=0.05$, pero el cociente de picos queda más cerca de 0.025, porque el pulso también sale y se reparte. Dividir los 0.3 bar del registro entre seis horas deja una pendiente media de 0.05 bar/h. Con $\mathcal{L}_1=0$, $\partial p/\partial t=-\mathcal{L}_5/2$ fija el signo de esa pendiente.

Lista de verificación antes de llevarlo al solver#

• Contar primero el tipo de frontera. Salida subsónica → 1 BC, salida supersónica → 0, entrada subsónica → 4, entrada supersónica → 5. Verificar con la figura de arriba.
• El ajuste de σ está ligado a la longitud del dominio. Si el dominio se duplica, $K$ también cambia. No fijar valores absolutos.
• No descartar los términos transversos a ciegas. En LES 3D la versión "relaxed" de Yoo & Im (2007) es más estable.
• NSCBC solo no basta cuando un choque golpea la salida. Combinar con un buffer zone (sponge layer) para absorber la reflexión no lineal.
• Los flujos viscosos requieren un tratamiento aparte. NSCBC no se ocupa de ellos por sí solo.

Con σ en cero la media se va de largo#

$\mathcal{L}_1=0$ deja salir el sonido y conserva $\partial p/\partial t=-\mathcal{L}_5/2$. Si el promedio temporal de $\mathcal{L}_5$ tiene un sesgo positivo pequeño, la presión no tiene techo. Los 0.3 bar a las seis horas son esa integral, escrita en el registro. Meter $\sigma\approx 0.25$ en $K=\sigma(1-M_{\max}^{2})c/L$ tira de la deriva, y con $\tau=L/c$ también devuelve cerca de un cuarto de la amplitud.

Subir $\sigma$ hacia 5 endurece otra vez la frontera y el pulso rebota. Entre 0 y 5, el 0.05 del código y el 0.25 habitual son dos puntos de la misma recta. Esta página llega hasta "hay una sola onda entrante". El ajuste más fino de ese único coeficiente está en cómo σ en la salida fija la reflexión.