[Reseña] Cómo hacer que un flujo no lineal se comporte como lineal — relajación de Jin–Xin y Suliciu

Al programar un solucionador de Riemann exacto (el dispositivo que resuelve la estructura de ondas en la cara de una celda) durante suficiente tiempo, la ecuación de estado detiene el trabajo. Hay que diagonalizar el jacobiano del flujo no lineal $f(u)$ en cada paso. Un gas real no entrega la estructura de autovalores en forma cerrada. En 1995 Shi Jin y Zhouping Xin invirtieron la dirección. La ecuación no lineal no se empuja tal cual. El flujo pasa a ser una incógnita nueva, se arma un sistema lineal y luego se relaja hacia la ecuación original. Eso es la relajación. Thomann et al. (2019) llevaron ese esqueleto hasta Euler de todas las velocidades (all-speed). El rodeo es más firme en un sitio y se quiebra en otro. Ese precio es toda la nota.

Artículo — Andrea Thomann, Markus Zenk, Gabriella Puppo, Christian Klingenberg, An all speed second order IMEX relaxation scheme for the Euler equations, arXiv:1907.08398 (2019).

Al bajar el control de velocidad de relajación a por debajo de 1, el choque de Burgers con estado izquierdo 1 y derecho $-0.4$ empieza a oscilar.

Cerca de 1.3 el frente se sostiene. Subir solo eps deja la velocidad quieta y aumenta únicamente el espesor.

Por qué un flujo no lineal complica el problema de Riemann#

Considérese la ley de conservación escalar $u_t + f(u)_x = 0$. Aquí $u$ es la cantidad conservada y $f(u)$ es el flujo. Resolver un tubo de choque en una sola celda en cada cara exige primero la velocidad de onda $f'(u)$.

El problema es el flujo real. En las ecuaciones de Euler compresibles el flujo lleva la presión, y la presión es una función no lineal de la densidad y la energía interna. Con una ecuación de estado de gas real (real gas), los autovalores y autovectores dejan de salir de forma analítica. Diagonalizar el jacobiano numéricamente en cada paso resulta caro. En el régimen de bajo Mach, las ondas acústicas escalan como $1/M$, de modo que el paso de tiempo explícito se ahoga.

De ahí la pregunta: ¿se puede posponer la no linealidad en vez de enfrentarla de frente en cada paso?

El ascenso de Jin–Xin — elevar el flujo a incógnita#

Su respuesta: en lugar de calcular el flujo $f(u)$ como un valor, se introduce una nueva variable $v$ que lo contenga, y se agrega una fuente que arrastra lentamente $v$ hacia $f(u)$.

Aquí $u$ es la cantidad conservada original, $v$ es la variable de relajación que reemplaza al flujo, $a$ es una velocidad de relajación constante (relaxation speed) y $\varepsilon$ es el tiempo de relajación. El lado izquierdo es completamente lineal. La matriz de coeficientes tiene autovalores fijos $\pm a$: la diagonalización no lineal desaparece.

Cuando $\varepsilon \to 0$, la fuente del lado derecho fuerza $v = f(u)$ y la primera ecuación vuelve a colapsar en la ley de conservación original. Es decir, este sistema lineal de $2\times 2$ es una aproximación viscosa de la ecuación no lineal. Como sus velocidades de onda son las constantes $\pm a$, el problema de Riemann se resuelve una sola vez, a mano. La figura de arriba es este sistema sobre Burgers, $f(u)=u^2/2$.

La condición subcaracterística — con a pequeña, el sistema explota#

La velocidad de relajación $a$ no se elige libremente. Para que la aproximación sea estable, debe cumplir la condición subcaracterística (condición de Whitham):

$f'(u)$ es la velocidad de onda verdadera de la ecuación original. La condición dice que las velocidades congeladas $\pm a$ que carga el sistema de relajación deben encerrar siempre esa velocidad verdadera. Para Burgers, $f'(u)=u$, así que basta con $a \ge \max|u|$.

En el gráfico de abajo conviene variar la velocidad de relajación $a$ y el rango de la solución $\max|u|$. El instante en que la curva naranja $f'(u)$ abandona la banda $\pm a$ es el umbral de inestabilidad.

Al bajar $a$ por debajo de la mayor velocidad de onda de la solución, la curva atraviesa la banda. Entonces la aproximación de relajación produce difusión negativa en vez de difusión física. La información fluye en la dirección equivocada y el sistema estalla. El motivo se ve con precisión en el desarrollo de la siguiente sección.

La difusión oculta, vía Chapman–Enskog#

¿Qué ocurre cuando $\varepsilon$ es pequeño pero no nulo? Al expandir $v = f(u) + \varepsilon v_1 + \cdots$ (una expansión de Chapman–Enskog) y sustituir en el sistema de relajación, surge una ecuación efectiva:

El lado derecho es la difusión que la relajación inyecta a escondidas. Su coeficiente es $a^2 - f'(u)^2$. Para que no sea negativo hace falta exactamente $a \ge |f'(u)|$. Así que la condición subcaracterística era, en realidad, la condición de que la viscosidad artificial de la relajación nunca se vuelva negativa.

El compromiso queda a la vista. Una $a$ grande da estabilidad holgada, pero un $a^2 - f'^2$ grande difumina la solución. Ceñir $a$ al límite subcaracterístico da nitidez pero es riesgoso. Y $\varepsilon$ fija la magnitud global de la difusión: por eso subir $\varepsilon$ en la simulación anterior engrosa el choque.

El margen es la viscosidad. Los datos de Riemann del código son $u=1$ a la izquierda y $u=-0.4$ a la derecha, así que $\max|u|$ parte de 1 y el piso también es 1. Con $a=1.3$, $a^2=1.69$. A la izquierda, $f'=1$ y el coeficiente de difusión vale $1.69-1=0.69$. A la derecha, $f'=-0.4$, de modo que $1.69-0.16=1.53$. El mismo $\varepsilon=10^{-4}$ emborrona más el estado derecho. Con $a=0.8$, $a^2=0.64$ y el coeficiente izquierdo es $0.64-1=-0.36$. El signo se invierte. La información vuelve a través del frente y la amplitud no se queda en el valor inicial 1.

Suliciu — relajar solo la presión, no todo el flujo#

Jin–Xin relaja todas las componentes del flujo. Elegante, pero la difusión es excesiva. Lo que se prefiere más en la práctica es la relajación tipo Suliciu. En las ecuaciones de Euler, la única no linealidad genuina es la presión, así que se relaja solo la presión a una nueva variable $\pi$.

Aquí $\rho$ es la densidad, $u$ la velocidad, $p$ la presión real y $\pi$ la presión relajada. La condición subcaracterística pasa a ser $a > \rho\sqrt{\partial_\rho p}$, encerrando la velocidad del sonido. Los autovalores del sistema relajado son $u,\ u\pm a/\rho$, todos linealmente degenerados (linearly degenerate): solo quedan ondas tipo contacto, fáciles de manejar.

Aquí empieza el aporte de Thomann et al. (2019). Apuntando al régimen de bajo Mach, dividen la presión en una parte lenta y una parte acústica rápida. La parte lenta se trata de forma explícita; la parte acústica rápida se resuelve de forma implícita (implicit) sobre el sistema de relajación. Agregan una nueva variable de velocidad $\hat u$ para asegurar una difusión independiente del número de Mach, y el resultado es un esquema IMEX de segundo orden, asymptotic-preserving, que converge a Euler incompresible en el límite de bajo Mach. Gracias a la estructura linealmente degenerada, todo esto marcha sin una sola diagonalización no lineal.

Python — resolver Burgers mediante el sistema de relajación#

Conviene reproducir el esqueleto de la idea de relajación en numpy: la ecuación de Burgers a través del sistema de Jin–Xin. Al separar en variables características $r,s$, el lado izquierdo se convierte en dos advecciones simples.

import numpy as np
 
def burgers_flux(u):
    return 0.5 * u * u                     # f(u) = u^2 / 2
 
def relaxed_step(u, v, a, dx, dt, eps):
    # separacion caracteristica: r se mueve a +a, s a -a
    r = 0.5 * (u + v / a)
    s = 0.5 * (u - v / a)
    nu = a * dt / dx                       # numero CFL a*dt/dx
    # upwind de primer orden (periodico, con np.roll)
    r_new = r - nu * (r - np.roll(r, 1))   # onda hacia la derecha
    s_new = s + nu * (np.roll(s, -1) - s)  # onda hacia la izquierda
    u_new = r_new + s_new
    v_new = a * (r_new - s_new)
    # fuente de relajacion: arrastrar v hacia f(u) (integrador exponencial)
    kappa = 1.0 - np.exp(-dt / eps)
    v_new += (burgers_flux(u_new) - v_new) * kappa
    return u_new, v_new
 
def run_relaxation(u0, a, eps, cfl=0.9, t_end=0.3):
    n = u0.size
    dx = 2.0 / n
    u = u0.copy()
    v = burgers_flux(u)                    # partir de la variedad de equilibrio
    dt = cfl * dx / a                      # la subcaracteristica fija el CFL
    t = 0.0
    while t < t_end:
        u, v = relaxed_step(u, v, a, dx, dt, eps)
        t += dt
    return u
 
# datos iniciales de Riemann: izquierda 1.0, derecha -0.4 (un choque)
n = 400
x = np.linspace(-1, 1, n, endpoint=False) + 1.0 / n
u0 = np.where(x < 0.0, 1.0, -0.4)
 
u_ok = run_relaxation(u0, a=1.3, eps=1e-4)   # a >= max|u|=1.0  -> estable
u_bad = run_relaxation(u0, a=0.8, eps=1e-4)  # a <  max|u|      -> violada
 
print("a=1.3  peak |u| =", round(float(np.max(np.abs(u_ok))), 3))   # ~1.0
print("a=0.8  peak |u| =", round(float(np.max(np.abs(u_bad))), 3))  # estalla

La salida da testimonio directo de la condición subcaracterística. Con $a=1.3$ la amplitud máxima se queda cerca del valor inicial y se forma un choque limpio. Con $a=0.8$ la amplitud salta varias veces. No hubo diagonalización de un jacobiano no lineal. Solo advección lineal repetida a las velocidades constantes $\pm a$.

El número de pasos también lo fija $a$. Con $n=400$ en un intervalo de longitud 2, $\Delta x=0.005$. A CFL 0.9, $\Delta t=0.9\times 0.005/a$. Para $a=1.3$, $\Delta t\approx 3.46\times 10^{-3}$ y la razón con $t=0.3$ ronda 87. Para $a=0.8$, $\Delta t\approx 5.63\times 10^{-3}$, una marcha más gruesa. Bajar $a$ alarga el paso de tiempo. Con $\varepsilon=10^{-4}$, $\Delta t/\varepsilon$ supera 30 y el exponencial de $\kappa=1-\exp(-\Delta t/\varepsilon)$ desaparece. $\kappa=1$. Cada paso pega $v$ a $f(u)$. El espesor que queda es solo $\varepsilon(a^2-f'^2)$.

Una mirada crítica#

La relajación no es gratis. La difusión artificial $\varepsilon(a^2-f'^2)$ siempre acompaña, y cuanto más segura (grande) se hace $a$, más se difumina la solución. Imponer la condición subcaracterística exige una cota global de $\max|f'(u)|$, y cerca de choques fuertes o del vacío una cota conservadora sobredimensiona la difusión. La lógica que era limpia para escalares exige, para el Euler de gas real, una estimación local de $\rho\sqrt{\partial_\rho p}$ más maquinaria adicional para mantener positivas la presión y la densidad. El acople IMEX de Thomann et al. debe resolver una ecuación elíptica para la parte acústica, así que existe un cruce donde la ganancia de bajo Mach se compensa con el costo del solucionador lineal. Al reproducirlo, el ajuste conjunto de $\varepsilon$ y $a$ resulta más sensible de lo esperado.

Los solucionadores basados en densidad de OpenFOAM y Fluent no traen este solucionador de Riemann por relajación de forma directa, pero las familias HLLC y AUSM de solucionadores aproximados comparten la misma filosofía: simplificar la estructura de ondas para evitar la diagonalización. El solucionador de Suliciu está implementado en códigos de código abierto como SU2, como una opción que captura las discontinuidades de contacto de forma exacta.

Dónde se rompe una a pegada al límite inferior#

El piso de estos datos es 1. El código elige 1.3. Ese margen de 0.3 es el coeficiente de difusión 0.69 en el estado izquierdo. Al recortar el margen hasta cero, con $a=\max|u|$, el coeficiente sobre el frente vale 0. La captura se ve nítida. Un empujón de redondeo por debajo del piso invierte el signo.

El fallo que aparece después es congelar la cota en el primer paso. En cuanto $|u|$ pasa de 1, la $a=1.3$ que parecía segura ya está fuera de la condición subcaracterística. Suliciu traslada la misma prueba a $a>\rho\sqrt{\partial_\rho p}$. Si la velocidad del sonido salta de celda en celda y la cota llega un paso tarde, la difusión de la relajación ya es negativa, haya o no acústica implícita. Tocar solo la presión ocupa el mismo sitio que separar la advección del flujo de presión.