[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 termina poniendo la zancadilla. Hay que diagonalizar el jacobiano del flujo no lineal $f(u)$ en cada paso, y para una ecuación de estado complicada esa estructura de autovalores no tiene forma cerrada. En 1995, Shi Jin y Zhouping Xin dieron vuelta la idea por completo. En lugar de resolver la ecuación no lineal directamente, ascendieron el propio flujo a incógnita nueva, construyeron un sistema lineal alrededor de él y dejaron que ese sistema se relajara hacia la ecuación original. Esta entrada sigue la idea de relajación desde cero y luego observa cómo Thomann et al. (2019) la estiraron hasta un solucionador de Euler de todas las velocidades (all-speed). Al final quedará claro por qué desviarse por un sistema lineal resulta más rápido y robusto, y cuál es el precio.

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

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. Un esquema tipo Godunov resuelve un problema de Riemann en cada cara, y para ello necesita 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. La clave es que 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.

Conviene manipular los parámetros en la simulación de abajo. Resuelve la ecuación de Burgers ($f(u)=u^2/2$) mediante el sistema de relajación anterior.

Con la condición inicial shock, al bajar la velocidad de relajación a por debajo de 1.0 el perfil oscila y estalla. Al subir a lo suficiente, captura un choque limpio. La siguiente sección explica por qué.

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.

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. Nótese que nunca se diagonalizó un jacobiano no lineal: solo advección lineal repetida a las velocidades constantes $\pm a$.

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.

Lo que cambió este enfoque#

• Invirtió la dirección de la linealización. En lugar de aproximar una ecuación no lineal para volverla lineal, se construye un sistema lineal exacto y se lo relaja hacia la original. La diagonalización desaparece.
• La estabilidad se compra con una sola desigualdad. $a \ge |f'(u)|$: la condición subcaracterística es justamente lo que impide que la viscosidad artificial se vuelva negativa.
• La estructura linealmente degenerada abre la práctica. La relajación tipo Suliciu se extiende a bajo Mach, gas real y flujos multifásicos, con el esquema IMEX de todas las velocidades de Thomann et al. como ejemplo insignia.