La salida está abierta y la onda regresa — NSCBC y el precio de σ

La salida está abierta. Aun así, la onda regresa. Basta con tratar la frontera de salida de un código compresible mediante una extrapolación descuidada y calcular una llama o una zona turbulenta: en pleno centro del dominio crece una oscilación sin causa física. Al medir su periodo, casi siempre coincide con la longitud del dominio dividida por la velocidad del sonido. La caja de cálculo se ha convertido en un tubo resonante. En la frontera que fabrica esa resonancia se separa lo que se calcula de lo que hay que inventar. El coeficiente que mueve ese invento tiene un precio, y el precio está en la tabla.

Al llevar el deslizador de Mach desde un valor subsónico positivo más allá de lo sónico, las flechas que entraban salen. Con el signo negativo, el número de condiciones de esa misma cara salta.

Una salida subsónica lleva una condición. Pasado lo sónico, la cuenta es cero. Al volver la cara una entrada, pasa a cuatro. La onda es un número inventado solo mientras la flecha apunta hacia el interior.

El borde de la caja no es física#

Un dominio físico no tiene borde. Más allá de la salida de la cámara de combustión el espacio continúa. La malla, en cambio, tiene que terminar en algún punto, y la última celda no tiene vecina. Sin vecina no hay derivada, y sin derivada no hay ecuación de gobierno que avanzar.

En los códigos RANS el problema estuvo oculto mucho tiempo. La viscosidad turbulenta y la artificial son grandes, así que una onda mal fabricada en la frontera muere a las pocas celdas. Con LES y DNS las cuentas cambian. La viscosidad artificial es casi nula y la turbulenta está en su mínimo. El error creado por la frontera no se disipa: cruza el dominio y vuelve.

La receta que Poinsot y Lele articularon en 1992 invierte el enfoque. En lugar de extrapolar variables en la frontera, se cuentan las ondas que la atraviesan y se fija la amplitud de cada una. Extiende las condiciones de contorno características de Euler (ECBC) a Navier–Stokes con sus términos viscosos, de ahí el nombre NSCBC (Navier–Stokes Characteristic Boundary Conditions). No se usa ni una línea de extrapolación.

Lo que una frontera puede contar y lo que debe inventar#

La frontera queda en $x_1 = L$. Al reagrupar los términos en dirección $x_1$ como ondas, la ecuación de continuidad queda:

Aquí $\rho$ es la densidad, $u_i$ las componentes de velocidad y $d_1$ reúne las contribuciones normales a la frontera. Ese vector $d$ es el producto del análisis característico, y dentro viven las amplitudes de onda $\mathcal{L}_i$.

$c$ es la velocidad local del sonido ($c^2 = \gamma p / \rho$) y $p$ la presión. $\mathcal{L}_1$ es la variación de amplitud de la onda acústica que corre en el sentido negativo de $x_1$, y $\mathcal{L}_5$ la que corre en el positivo. Las tres restantes viajan con el fluido: $\mathcal{L}_2$ transporta entropía, $\mathcal{L}_3$ y $\mathcal{L}_4$ las velocidades transversales $u_2$ y $u_3$, todas a velocidad $u_1$.

Lo que separa los casos es el signo de esas velocidades. Si una onda sale del dominio, su amplitud se calcula desde los puntos interiores. Si entra, esa información no está en ninguna parte de la solución: hay que inventarla. El número de ondas entrantes es exactamente el número de condiciones de contorno físicas que admite esa frontera.

Cómo las características parten el sonido en todo el dominio está en por qué un flujo supersónico no conoce su aguas arriba. Aquí solo se cuenta cuántas de esas líneas corta una cara de salida.

Al subir $M$ de $+0.3$ a $+1.4$, la flecha que cambió de sentido arriba es $\mathcal{L}_1$. La cuenta pasa de uno a cero. Con signo negativo la cara es una entrada y la cuenta es cuatro. Un código que corría con salida subsónica y diverge al volverse supersónica suele imponer una condición más de las que este registro admite.

LODI: la regla que fabrica una amplitud#

Inventar una amplitud entrante exige una justificación. En cada punto de la frontera, NSCBC construye un sistema local unidimensional no viscoso, eliminando todos los términos transversales, viscosos y de reacción. Son las relaciones LODI (Local One Dimensional Inviscid).

Las relaciones LODI no son condiciones físicas ni son las ecuaciones que se resuelven. Existen solo para estimar las $\mathcal{L}_i$ entrantes. El procedimiento tiene tres pasos: eliminar del sistema cada ecuación de conservación cuya variable esté impuesta físicamente, usar la relación LODI correspondiente para escribir la $\mathcal{L}_i$ desconocida en función de las conocidas, y avanzar el resto de variables con las ecuaciones que quedan.

Una salida perfectamente no reflectante toma aquí la decisión más simple posible: declara que ninguna onda acústica llega desde fuera.

Onda entrante nula, reflexión nula. Parece limpio. Pero en esa afirmación no aparece por ningún lado la presión exterior $p_\infty$.

Lo que cuesta σ = 0: la presión que nunca vuelve a p∞#

Sin $p_\infty$, la frontera no tiene idea de cuál debería ser la presión. El calor liberado dentro del dominio la sube, y no hay en ninguna parte una fuerza restauradora que devuelva ese desplazamiento. El problema deja de estar bien planteado.

La corrección de Rudy y Strikwerda consiste en no anular la onda entrante, sino atarla a la diferencia de presión:

$L$ es una longitud característica del dominio, $M$ el número de Mach máximo y $\sigma$ el único parámetro libre de toda la receta. Con $\sigma = 0$ se vuelve al caso perfectamente no reflectante. Al subir $\sigma$, la frontera empieza a arrastrar la presión hacia $p_\infty$.

El conducto cerrado, con sigma y source q, queda justo después de la tabla impresa.

La ventana estrecha entre dos fracasos#

$\sigma$ vive entre dos fracasos que tiran en sentidos opuestos.

Con $\sigma$ pequeño la presión media se va a la deriva; con $\sigma$ grande la frontera se endurece hasta devolver energía acústica al interior. Los valores que Poinsot y Lele emplearon fueron $\sigma \approx 0.25$, y $\sigma = 0.58$ para el coeficiente equivalente en el B1 basado en extrapolación. Ninguno de los dos sale de la teoría: ambos se eligieron entre los dos fracasos.

La dependencia con la frecuencia explica por qué es un compromiso. Para una onda acústica de frecuencia angular $\omega$, esta frontera refleja con

Cuanto más baja la frecuencia, mayor la reflexión. Es decir, $\sigma$ funciona como un filtro que agarra las frecuencias bajas y deja pasar las altas. La presión media es la componente $\omega \to 0$, así que queda agarrada; la acústica que se quiere evacuar pasa. La ventana en la que esa separación funciona es estrecha.

Medir reflexión y desplazamiento con código#

En acústica lineal unidimensional el estado se parte exactamente en dos amplitudes características: $A^\pm = p' \pm \rho c u'$, que viajan a $\pm c$. Con $\Delta t = \Delta x / c$ cada una se desplaza exactamente una celda por paso, de modo que toda ondulación en pantalla proviene de la condición de contorno y no del esquema.

import numpy as np
 
N, C, L = 240, 1.0, 1.0
DX = L / N
DT = DX / C                       # desplazamiento exacto de una celda: sin difusion del esquema
 
 
def outlet_relax_k(sigma, mach=0.0):
    """Coeficiente de relajacion NSCBC K = sigma (1 - M^2) c / L"""
    return sigma * (1.0 - mach ** 2) * C / L
 
 
def duct_step(ap, am, am_b, k_relax, q):
    """A+ una celda a la derecha, A- una a la izquierda. En la salida hay que inventar A-."""
    ap[1:] = ap[:-1].copy()
    ap[0] = 0.0
    am[:-1] = am[1:].copy()
 
    p_b = 0.5 * (ap[-1] + am[-1])
    am_b -= DT * k_relax * p_b     # L1 = K (p - p_inf)
    am[-1] = am_b
 
    ap[0] = am[0]                  # extremo cerrado en x = 0 (u = 0)
    ap += q * DT                   # liberacion de calor debil y uniforme
    am += q * DT
    return am_b
 
 
def measure_outlet(sigma, q, steps, pulse):
    ap, am, am_b = np.zeros(N), np.zeros(N), 0.0
    x = (np.arange(N) + 0.5) * DX
    if pulse:
        ap += np.exp(-((x - 0.30) / 0.09) ** 2)
    k = outlet_relax_k(sigma)
    refl = 0.0
    for n in range(steps):
        am_b = duct_step(ap, am, am_b, k, q)
        if pulse and n > 0.85 * N:
            refl = max(refl, np.abs(am[:-8]).max())
    return refl, float(np.mean(0.5 * (ap + am)))
 
 
for sigma in (0.0, 0.25, 1.0, 4.0, 10.0):
    r, _ = measure_outlet(sigma, q=0.0, steps=650, pulse=True)
    _, p = measure_outlet(sigma, q=0.3, steps=6000, pulse=False)
    print(f"sigma={sigma:5.2f}   reflejado={r * 100:5.1f}%   p media - p_inf={p:+.4f}")

Lo que imprime:

sigma= 0.00   reflejado=  0.0%   p media - p_inf=+0.3000
sigma= 0.25   reflejado=  1.9%   p media - p_inf=+0.0007
sigma= 1.00   reflejado=  7.3%   p media - p_inf=+0.0006
sigma= 4.00   reflejado= 24.5%   p media - p_inf=+0.0227
sigma=10.00   reflejado= 46.8%   p media - p_inf=-0.1142

Cuando el pulso de fire pulse llega a la salida, una curva roja $A^-$ que sube es la reflexión. Con source q alto y $\sigma = 0$, la curva blanca se queda por encima de $p_\infty$. La reflexión no está. La presión tampoco queda sujeta.

Con $\sigma = 0$ la reflexión desaparece por completo, y el desplazamiento de 0.3 creado por la liberación de calor se queda intacto. Con $\sigma = 0.25$ la reflexión baja del 2% y la presión se mantiene a menos de 0.0007 de $p_\infty$. Hasta aquí, lo esperado.

La sorpresa está en las dos últimas filas. Llevar $\sigma$ a 4 y a 10 sube la reflexión al 25% y al 47%, lo cual no extraña; pero el desplazamiento de presión vuelve a empeorar. Una frontera tan rígida empieza a resonar por su cuenta, y esa resonancia mueve la media. Subir $\sigma$ no es un canje en el que se compre control de presión pagando con reflexión. Pasado cierto punto se pierden las dos cosas.

La columna de reflexión lee 0,0, 1,9, 7,3, 24,5 y 46,8%. En el código $C = 1$, $L = 1$ y el Mach es 0, así que $K = \sigma \left( 1 - M^2 \right) \frac{c}{L}$ se reduce a $K = \sigma$. Al pasar $\sigma$ de 0,25 a 10, un factor 40, la reflexión pasa de 1,9 a 46,8. Unas 25 veces. No es lineal. $|R| = \frac{K}{\sqrt{K^2 + 4\omega^2}}$ se satura hacia 1 cuando K crece. El pulso no es una sola frecuencia: es una gaussiana de ancho 0,09. No hace falta invertir un $\omega$ a partir de estos porcentajes. Basta la dirección. Al crecer K la reflexión sube y luego se aplana.

La columna de presión no sigue esa monotonía. En la ejecución con $q = 0.3$, $\sigma = 0$ deja una media de +0,3000, el mismo tamaño que la fuente. $\sigma = 0.25$ deja +0,0007 y $\sigma = 1$ deja +0,0006. Frente a 0,3000, 0,0007 es el 0,23%. Hasta ahí la relajación casi borra el desplazamiento. Con $\sigma = 4$ el desplazamiento es +0,0227, unas 32 veces 0,0007. Con $\sigma = 10$ es −0,1142. El signo se invierte. El valor absoluto es unas 160 veces 0,0007. Una frontera que refleja el 46,8% no sujeta la media: la sacude. El suelo de esa columna está en 0,25 y 1,0. El 0,25 que usaron Poinsot y Lele es el borde izquierdo de ese suelo. 0,58 corresponde al coeficiente de B1, así que no cae tal cual sobre esta columna.

Las ondas que inventó la malla viajan al revés#

En la frontera no solo se refleja la acústica física. Las componentes con longitud de onda menor que unos cuatro espaciados de malla no son solución de Navier–Stokes: son artefactos de la discretización. Poinsot y Lele las llamaron "ondas q" para separarlas de las "ondas p" físicas.

Su firma es la velocidad de grupo. Incluso en la ecuación de advección unidimensional a velocidad $V$, la velocidad de grupo $u_g$ de las longitudes de onda cortas tiene signo opuesto a $V$. El flujo avanza hacia la derecha y el error numérico repta aguas arriba, hacia la izquierda. Peor aún: $|u_g / V|$ crece con el orden del esquema espacial, así que los códigos de alto orden quedan más expuestos, no menos.

Por eso un tratamiento de frontera hay que juzgarlo con dos coeficientes de reflexión: $A_p / A_1$ para las ondas físicas y $A_q / A_1$ para las numéricas. Cualquier tratamiento utilizable necesita $A_q / A_1 \ll 1$ en toda circunstancia; si además presume de no reflectante, necesita $A_p / A_1 \ll 1$. Basta arrancar el cálculo desde un campo inicial con gradientes abruptos para generar ondas q, y en una DNS nada las elimina después.

Si el periodo es 2L/c, la sospecha cae primero sobre σ#

$\sigma$ no es una perilla de ajuste. Un extremo es una presión sin anclaje. El otro es una caja convertida en resonador. Cerca de 0,25 la frontera agarra las frecuencias bajas y deja pasar el resto.

Una oscilación sin causa visible empieza por un conteo: características entrantes frente a condiciones impuestas. Un periodo múltiplo de $2L/c$ es la frontera, no la física. Al bajar $\sigma$, si la oscilación se reduce, la causa era la frontera. Si la presión media se va, el cálculo cruzó al otro lado. Contar esas ondas antes de tocar la salida es el motivo de una pasada anterior que deja NSCBC como salida no reflectante.