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. Hoy toca ver qué puede calcular realmente una frontera, qué se ve obligada a inventar, y qué cuesta el coeficiente que gobierna esa invención.

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#

Sitúese la frontera 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$.

Todo el juego está en 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.

Conviene mover el número de Mach en el diagrama siguiente. La dirección con que cada una de las cinco características cruza la frontera cambia en tiempo real.

Al subir $M$ de $+0.3$ a $+1.4$, la $\mathcal{L}_1$ que entraba cambia de sentido y el número de condiciones necesarias cae de uno a cero. Con el signo invertido, esa misma cara pasa a ser entrada y la cuenta salta a cuatro. Cuando un código que corría sin problemas con salida subsónica diverge al volverse supersónica, casi siempre está imponiendo una condición más de las que este registro permite.

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$.

Vale la pena probarlo en la simulación de abajo. El conducto está cerrado a la izquierda y abierto a la derecha. Con fire pulse se lanza una onda de presión; después basta mover el deslizador sigma.

Hay dos cosas que observar. Primero: cuando el pulso llega a la salida, ¿se levanta la curva roja, la amplitud entrante $A^-$? Ese bulto es la reflexión. Segundo: al subir source q, conviene mirar el historial de presión inferior. Con $\sigma = 0$, la curva blanca se estaciona por encima de la línea $p_\infty$ y ahí se queda. La reflexión desapareció, pero la presión se perdió.

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

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.

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.

Una línea para elegir la condición de salida#

$\sigma$ no es una perilla de ajuste, sino una coordenada entre dos fracasos: una presión sin anclaje en un extremo, una caja de cálculo convertida en resonador en el otro. Si 0.25 se recomienda una y otra vez es porque en ese punto agarra las frecuencias bajas y deja pasar el resto.

Cuando una salida empieza a oscilar sin motivo aparente, el orden es este. Primero, contar las características entrantes en esa cara y comprobar que el número de condiciones impuestas coincide. Luego, mirar si el periodo de la oscilación es múltiplo de $2L/c$: si lo es, se trata de reflexión en la frontera y no de física. Por último, bajar $\sigma$. Si la oscilación se reduce, la culpa era de la frontera; si la presión media empieza a irse, se ha cruzado al fracaso del otro lado.