Lo que se sacrifica al renunciar a la solución exacta de Riemann — el solver de Roe y el arreglo de entropía#

El problema de Riemann de Euler 1-D tiene solución exacta. Una ecuación no lineal para la presión de la región estrella, un lazo de Newton, y se cierra. Los códigos compresibles de producción no usan esa solución en cada cara. La respuesta está en la mano. Se descarta. Costo. Fuera del gas ideal, la solución exacta puede no existir.

Al mover los controles de densidad y velocidad a izquierda y derecha, los pesos sqrt(rho) y las tres velocidades $\hat{u}-\hat{c}$, $\hat{u}$, $\hat{u}+\hat{c}$ se recalculan en la cara.

Cuando $\hat{u}-\hat{c}<0<\hat{u}+\hat{c}$, el abanico queda a horcajadas sobre la cara. Subsónico. Si ambas velocidades se empujan con fuerza hacia el mismo lado, las tres ondas comparten signo y el estado pasa a supersónico.

Por qué la solución exacta desapareció de la práctica#

Un solver exacto de Riemann (Godunov) encuentra la raíz de una ecuación no lineal en cada cara. Con millones de celdas, ese bucle de búsqueda de raíces corre millones de veces por paso. Peor aún, la solución exacta queda atada a la ecuación de estado del gas ideal. Al pasar a un gas real o a una mezcla bifásica, la solución exacta deja de existir.

Los solvers aproximados de Riemann esquivan este problema. Linealizan localmente el problema de Riemann no lineal. Un solo conjunto de fórmulas algebraicas entrega el flujo sin ninguna iteración. El precio es una pequeña pérdida de exactitud.

Condensar el jacobiano en un único promedio#

La idea de Roe es simple. Se aproxima la diferencia de flujo a través de una cara mediante una única matriz constante $\hat{A}$ que actúa sobre los dos estados $q_L$, $q_R$.

Aquí $q$ es el vector de variables conservadas y $f$ es el flujo físico. Esta condición es decisiva. Si los dos estados están conectados por una sola onda (un choque o una discontinuidad de contacto), entonces $\hat{A}$ propaga esa onda de forma exacta, incluso con amplitud finita. Así, el solver es aproximado, pero exacto frente a un choque o un contacto puros.

El detalle está en cómo elegir $\hat{A}$. Un promedio descuidado rompe la condición anterior.

De dónde viene el promedio con sqrt(rho)#

La respuesta de Roe es un promedio ponderado por la raíz cuadrada de la densidad.

Aquí $u$ es la velocidad, $H$ es la entalpía específica total, y los circunflejos marcan los valores promediados en la cara. La velocidad del sonido de Roe se obtiene como $\hat{c}=\sqrt{(\gamma-1)(\hat{H}-\hat{u}^2/2)}$.

¿Por qué la raíz cuadrada de la densidad? Si se escriben el estado $q$ y el flujo $f$ en las componentes del vector de parámetros $W=\sqrt{\rho}\,(H,u,1)^T$, ambos resultan perfectamente cuadráticos. Al ser cuadráticos, las derivadas de $q$ y $f$ respecto a $W$ son lineales, y la integral de trayectoria que construye $\hat{A}$ se evalúa de forma exacta. Lo que emerge es precisamente este promedio ponderado por sqrt(rho).

El denominador del peso es $\sqrt{\rho_L}+\sqrt{\rho_R}$. Si una densidad es cuatro veces la otra, su raíz es el doble, y esa velocidad y esa entalpía tiran del promedio con el doble de fuerza. Un promedio aritmético no cumple $\hat{A}(q_R-q_L)=f(q_R)-f(q_L)$ en un salto finito. Sin sqrt(rho) las velocidades de onda ya fallan delante de un choque puro. La velocidad del sonido es $\hat{c}=\sqrt{(\gamma-1)(\hat{H}-\hat{u}^2/2)}$. Con $\gamma=1.4$, $\gamma-1=0.4$. Lo que queda al restar la energía cinética $\hat{u}^2/2$ de la entalpía total es el cuadrado de la velocidad del sonido. El lugar donde cambian de signo las tres velocidades es el mismo corte en tres ondas que la estructura propia de Euler en una normal arbitraria.

Descomponer en ondas y volver a sumarlas#

Ensamblar el flujo lleva tres pasos. Primero se descompone el salto $\Delta q = q_R - q_L$ en una suma de los tres autovectores $\hat{K}_k$; el tamaño de cada componente es la intensidad de onda $\alpha_k$. Luego se transporta cada onda contra la corriente según el signo de su autovalor $\hat{\lambda}_k$.

El primer término es un promedio centrado; el segundo es la disipación upwind ponderada por las magnitudes de los autovalores. Esa estructura se traduce directamente a código. La prueba es el problema de Shu–Osher: un choque a Mach 3 se lanza contra un campo de densidad sinusoidal y deja tras de sí estructura de alta frecuencia. Es justo donde la baja disipación numérica de Roe rinde sus frutos.

import numpy as np
 
gamma = 1.4
 
def phys_flux(U):                       # conservadas U=(rho, rho*u, E) -> flujo físico
    rho = U[0]; u = U[1] / rho; E = U[2]
    p = (gamma - 1) * (E - 0.5 * rho * u * u)
    return np.array([rho * u, rho * u * u + p, u * (E + p)])
 
def roe_flux(UL, UR, delta):            # flujo aproximado de Roe (+ arreglo de entropía de Harten)
    rhoL, rhoR = UL[0], UR[0]
    uL, uR = UL[1] / rhoL, UR[1] / rhoR
    pL = (gamma - 1) * (UL[2] - 0.5 * rhoL * uL * uL)
    pR = (gamma - 1) * (UR[2] - 0.5 * rhoR * uR * uR)
    HL = (UL[2] + pL) / rhoL; HR = (UR[2] + pR) / rhoR
    sL, sR = np.sqrt(rhoL), np.sqrt(rhoR)      # pesos sqrt(rho)
    u = (sL * uL + sR * uR) / (sL + sR)        # velocidad promediada de Roe
    H = (sL * HL + sR * HR) / (sL + sR)        # entalpía promediada de Roe
    c = np.sqrt((gamma - 1) * (H - 0.5 * u * u))   # velocidad del sonido promediada de Roe
    rr = sL * sR
    drho, dp, du = rhoR - rhoL, pR - pL, uR - uL
    alpha = np.array([(dp - rr * c * du) / (2 * c * c),   # intensidades de onda alpha_k
                      drho - dp / (c * c),
                      (dp + rr * c * du) / (2 * c * c)])
    lam = np.array([u - c, u, u + c])          # autovalores (velocidades de onda)
    K = np.array([[1, u - c, H - u * c],       # autovectores derechos
                  [1, u,     0.5 * u * u],
                  [1, u + c, H + u * c]])
    al = np.abs(lam)
    small = al < delta                         # arreglo de entropía: piso sobre |lambda|
    al[small] = (lam[small] ** 2 + delta ** 2) / (2 * delta)
    diss = (al * alpha) @ K
    return 0.5 * (phys_flux(UL) + phys_flux(UR)) - 0.5 * diss
 
def run_shu_osher(N=400, tmax=1.8, cfl=0.4, delta=0.1):
    x = np.linspace(0, 10, N); dx = x[1] - x[0]
    rho = np.where(x < 1, 3.857143, 1 + 0.2 * np.sin(5 * x))   # choque + densidad sinusoidal
    u   = np.where(x < 1, 2.629369, 0.0)
    p   = np.where(x < 1, 10.33333, 1.0)
    U = np.array([rho, rho * u, p / (gamma - 1) + 0.5 * rho * u * u])
    t = 0.0
    while t < tmax:
        r = U[0]; v = U[1] / r; pp = (gamma - 1) * (U[2] - 0.5 * r * v * v)
        dt = cfl * dx / np.max(np.abs(v) + np.sqrt(gamma * pp / r))
        dt = min(dt, tmax - t)
        F = np.zeros((3, N + 1))
        for i in range(1, N):
            F[:, i] = roe_flux(U[:, i - 1], U[:, i], delta)
        F[:, 0] = phys_flux(U[:, 0]); F[:, N] = phys_flux(U[:, N - 1])
        U[:, 1:N - 1] -= dt / dx * (F[:, 2:N] - F[:, 1:N - 1])
        t += dt
    return x, U[0]
 
x, rho = run_shu_osher()
print(f"t=1.8  min rho={rho.min():.3f}  max rho={rho.max():.3f}")   # -> min~0.81  max~4.08

Unas cuarenta líneas bastan para un solver compresible completo. Sin iteración de Newton, sin solución exacta de Riemann. delta es el parámetro del arreglo de entropía que se trata a continuación.

El campo inicial de Shu–Osher es $\rho=3.857143$, $u=2.629369$, $p=10.33333$ para $x<1$, y $\rho=1+0.2\sin(5x)$, $u=0$, $p=1$ a la derecha. $N=400$ en $[0,10]$, CFL 0.4, $t=1.8$, $\delta=0.1$. La salida deja la densidad mínima cerca de 0.81 y la máxima cerca de 4.08. La densidad izquierda 3.857 sube un poco hacia 4.08, y un valle de la senoide baja hacia 0.8. Ese intervalo, sin Newton, es lo que sostiene esta corrida.

Roe rompe la condición de entropía#

El solver de Roe trata cada onda como un salto. Incluso una rarefacción se aproxima como un tren de pequeños choques. Por lo general eso es inofensivo. Pero cuando una rarefacción contiene un punto sónico, la situación cambia. En ese punto un autovalor $\hat{\lambda}_k$ cambia de signo y pasa por $0$.

Cuando $|\hat{\lambda}_k|\to 0$, la disipación upwind de esa onda desaparece. El esquema fija un choque de expansión estacionario en lugar de abrir un abanico suave. Satisface la condición de Rankine–Hugoniot pero viola la condición de entropía: una solución no física.

El remedio es el arreglo de entropía de Harten. Coloca un piso bajo la magnitud del autovalor cerca de $0$.

Con $|\hat{\lambda}_k|=0$ la fórmula del arreglo vale $\delta/2$, porque $(0+\delta^2)/(2\delta)=\delta/2$. Si $\delta$ es 0, el piso también es 0 y la disipación upwind desaparece. El delta=0.1 del código pone 0.05 sobre ese piso. Un 5–10 % de $(|\hat{u}|+\hat{c})$ es el ancho antes de que un contacto se difumine de más.

Al subir el control de $\delta$ desde 0 en la expansión de Burgers $u_t+(u^2/2)_x=0$ con $u_L=-0.5$ y $u_R=1.0$, se ve si el quiebre azul en $x=0$ se mueve. El punto sónico $u=0$ queda en el centro del intervalo.

Con $\delta=0$ aparece un quiebre azul —el choque de expansión estacionario— fijo en $x=0$. Al subir $\delta$, la curva numérica se desliza hasta el abanico de rarefacción exacto en ámbar. Un 5–10 % del $(|\hat{u}|+\hat{c})$ local suele ser seguro. Si el valor es demasiado grande, los contactos se difuminan.

Más barato: HLL y HLLC#

Si Roe resulta pesado, conviene conservar menos ondas. HLL (Harten–Lax–van Leer) mantiene solo las dos ondas acústicas. Descarta la onda de contacto central.

Aquí $S_L$, $S_R$ son estimaciones de las velocidades de onda más externas, izquierda y derecha. HLL es robusto y preserva bien la densidad positiva. Pero no logra mantener nítidos los contactos: no hay onda central.

La C de HLLC corresponde a la onda de Contacto central. Restaura la onda que HLL descartó: tres ondas, cuatro regiones de estado constante. Los tres solvers se sitúan sobre un espectro según el número de ondas.

El valor por defecto en producción suele ser HLLC. Mantiene vivos los contactos y a la vez facilita imponer la positividad de densidad y presión. Roe tiene una resolución excelente, pero el arreglo de entropía es obligatorio y hay que vigilar el fenómeno del carbúnculo en choques fuertes alineados con la malla.

Con δ en cero, el choque de expansión se queda en la cara#

La ponderación sqrt(rho) sale de la parametrización $W=\sqrt{\rho}\,(H,u,1)^T$, la que vuelve cuadráticos a $q$ y $f$. Sin esa integral de trayectoria, la condición de salto de $\hat{A}$ se rompe con amplitud finita. Si $|\hat{\lambda}_k|\to 0$ se deja pasar en el punto sónico, Rankine–Hugoniot sigue en pie y la entropía no. Roe sin el piso de Harten fija el choque de expansión como solución estacionaria.

Cuántas ondas conservar es la decisión siguiente. HLL deja dos ondas acústicas y difumina el contacto. HLLC devuelve la onda de contacto y cuatro regiones de estado. Roe llega a cinco ondas en 3-D, pero el arreglo de entropía es obligatorio, y un choque fuerte alineado con la malla sigue exigiendo vigilar el carbúnculo. A bajo Mach las escalas de la difusión artificial se separan por su cuenta. El registro de Hope–Collins que parte el Roe de bajo Mach en tres escalados es esa separación.