[Reseña] Cuando la velocidad del sonido diverge y el paso de tiempo sobrevive — AP-IMEX para gases no ideales de Orlando–Bonaventura (2024)

Los modelos atmosféricos corren a Mach 0.001. Las explosiones cósmicas, a Mach 100. Un solo solver para los dos extremos lleva medio siglo como deseo. El artículo de Orlando y Bonaventura de 2024 (arXiv:2402.09252v4) empuja ese deseo hasta los gases no ideales. Si la discretización temporal trata solo el término de presión de forma implícita, el paso de tiempo se suelta de la velocidad del sonido. Esa soltura se mantiene en la SG-EOS y en una EOS cúbica general (Van der Waals, Peng–Robinson).

Al bajar el control de $M$ hasta 0.05, la velocidad del sonido llega a 20 y, con el mismo paso, solo la CFL explícita supera 1 en un paso.

La curva roja superior es forward Euler explícito. La cian inferior es AP-IMEX sobre las mismas ecuaciones y el mismo $\Delta t$. Cuando la zona roja queda cubierta por la divergencia, la traza cian sigue llevando el pulso gaussiano hacia ambos lados.

Resumen en 30 segundos#

• Autores: Giuseppe Orlando, Luca Bonaventura
• Filiación: École polytechnique / Politecnico di Milano
• arXiv: 2402.09252v4 (v4, 22-10-2025)
• Objetivo: integración temporal asintótica-preservante (AP) para todo Mach y cualquier EOS
• Discretización espacial: Discontinuous Galerkin
• Validación: vórtice isentrópico, vórtice de Gresho, inestabilidad RT, choque transónico — extendidos a SG-EOS y EOS cúbica

Dos puntos donde un esquema explícito se rompe#

La integración temporal explícita se rompe dos veces cuando chocan dos escalas de tiempo.

La velocidad del sonido va primero. En Euler compresible las señales viajan a $|\mathbf u| \pm c$. Si $M_{loc} = |\mathbf u|/c$ baja, $c$ aplasta a $|\mathbf u|$. La CFL explícita queda en $\Delta t \le \Delta x / c$. Con $M = 0.01$ el paso temporal cae por un factor 100. Atmósfera y océano caen ahí con frecuencia.

Después va la EOS. En gas ideal todo termina en $c^2 = \gamma p / \rho$. Con una EOS cúbica, $\partial p / \partial \rho |_s$ es una función no lineal de $\rho$. Si el esquema explícito yerra la velocidad del sonido, la celda cae en presión negativa y la EOS deja de estar definida.

El artículo ata ambos problemas. Llevar el término acústico al lado implícito abre la primera trampa. Escribir ese paso sin depender de la EOS abre la segunda.

Qué significa realmente "asintótica-preservante"#

AP (asymptotic-preserving) cabe en una frase. Cuando un modelo continuo $\mathcal F^{\varepsilon}$ converge a un límite $\mathcal F^0$ con $\varepsilon \to 0$, la discretización $\mathcal F^{\varepsilon}_{\Delta t}$ debe converger de forma consistente a $\mathcal F^0_{\Delta t}$ en el mismo límite. La cota de estabilidad no depende de $\varepsilon$.

Aquí $\varepsilon$ es $M$ y $\mathcal F^0$ es Euler incompresible. En la forma habitual:

$\bar\rho$ es el orden cero (límite incompresible), $\rho'$ la primera corrección, $\rho''$ la segunda. Un esquema AP conserva esa jerarquía en el discreto.

La consecuencia numérica es nítida. A $M = 10^{-4}$ la fluctuación de presión medida debe escalar como $\mathcal O(M^2) = 10^{-8}$. Un Roe explícito no lo entrega. Con cualquier malla produce ruido $\mathcal O(M)$ (Guillard–Viozat, 1999). AP-IMEX conserva la escala.

Sólo presión implícita — la separación IMEX#

La separación viene de Cordier–Degond–Kumbaro (2012). Al lado implícito pasa solo $\nabla p$ en la ecuación de momento. El resto queda explícito.

$\rho$ es densidad, $\mathbf u$ velocidad, $p$ presión, $E$ energía total específica, $\Delta t$ el paso temporal y $M$ el Mach de referencia.

$\nabla p^{n+1}$ acopla el momento, y ese momento reaparece en la energía. Al combinar ambas sale una ecuación elíptica (linealizada) para $p^{n+1}$. Resuelta esa ecuación se obtiene $\mathbf u^{n+1}$ y luego se actualiza $\rho^{n+1}$ de forma explícita. El paso queda libre de la velocidad del sonido y solo persiste la CFL de transporte $|\mathbf u| \Delta t / \Delta x \le 1$. La misma separación, escrita para la ecuación de ondas, está en un IMEX que no apaga el sonido a bajo Mach.

El artículo monta esta separación sobre un IMEX Runge–Kutta. Cada etapa lleva tablas de Butcher distintas para la parte explícita y la implícita, y el Teorema 3.4 las analiza a la vez para garantizar consistencia AP y estabilidad. No hay división de operadores, ni de flujos, ni relajación.

Más allá del gas ideal — SG y Peng–Robinson#

Dos ecuaciones de estado cargan la validación.

Stiffened gas (SG-EOS) — habitual en líquidos y sólidos débilmente compresibles.

$\gamma$ es la razón de calores específicos, $e$ la energía interna específica, $q_\infty$ una constante de energía de enlace y $\pi_\infty$ una constante de rigidez. Con $q_\infty = \pi_\infty = 0$ vuelve el gas ideal. La velocidad del sonido es $c^2 = \gamma(p + \pi_\infty)/\rho$.

EOS cúbica general — agrupa Peng–Robinson y Van der Waals en una fórmula.

$R$ es la constante específica del gas, $T$ la temperatura, $a(T)$ la atracción intermolecular y $b$ el covolumen. Con $r_1 = r_2 = 0$ queda Van der Waals; con $r_1 = -1-\sqrt 2,\ r_2 = -1+\sqrt 2$ aparece Peng–Robinson.

La aportación es directa. La prueba AP del Teorema 3.4 no asume nada sobre la forma de la EOS. Aislado el paso implícito de presión, sea cual sea la curva de $p$, el límite incompresible cae de forma consistente cuando $M \to 0$. SG y Peng–Robinson solo sirven de EOS de validación.

Esa generalidad es una restricción de laboratorio. Cápsulas de CO₂ supercrítico, propulsores líquidos y ciclos con condensación dependen en cantidad de los efectos no ideales. Con gas ideal la velocidad del sonido puede desplazarse cerca de un 30 %.

Un paso en NumPy, paso a paso#

Con solo la acústica lineal 1D del artículo, el paso implícito se ve como un número. Malla periódica, celdas alternadas y resolución de Fourier.

import numpy as np
 
def acoustic_imex_step(rho, u, dt, dx, M):
    """Un paso AP-IMEX para el sistema acústico lineal 1D.
 
    rho: perturbación de densidad en centros de celda, longitud N
    u:   velocidad en bordes de celda (i+1/2), longitud N (periódica)
    """
    N = len(rho)
    alpha = (dt / (M * dx)) ** 2
    # RHS de momento — solo el gradiente de presión se difiere a (n+1)
    rhs = u - (dt / (M * M * dx)) * (np.roll(rho, -1) - rho)
    # Helmholtz cíclica: (1 + 2α) u_i - α (u_{i+1} + u_{i-1}) = rhs_i
    # Las matrices circulantes diagonalizan en la base DFT — solución directa
    k = np.arange(N)
    eig = 1 + 2 * alpha * (1 - np.cos(2 * np.pi * k / N))
    u_new = np.real(np.fft.ifft(np.fft.fft(rhs) / eig))
    # Continuidad: rho^{n+1} = rho^n - dt/dx (u_{i+1/2} - u_{i-1/2})
    rho_new = rho - (dt / dx) * (u_new - np.roll(u_new, 1))
    return rho_new, u_new
 
 
def cubic_eos_pressure(rho, T, a, b, R, r1, r2):
    """EOS cúbica general — Peng–Robinson usa r1=-1-sqrt(2), r2=-1+sqrt(2)."""
    return rho * R * T / (1 - rho * b) \
         - a * rho ** 2 / ((1 - rho * b * r1) * (1 - rho * b * r2))
 
 
# Demo: con M = 0.01, dt es 50 veces la CFL acústica y aún así seguro
N, dx = 64, 1 / 64
x = (np.arange(N) + 0.5) / N
rho = 0.1 * np.exp(-120 * (x - 0.5) ** 2)
u = np.zeros(N)
M, dt = 0.01, 0.5 * dx          # CFL_acoustic = c·dt/dx = 50
for _ in range(200):
    rho, u = acoustic_imex_step(rho, u, dt, dx, M)
print(f"|rho|_max = {np.abs(rho).max():.4f}  (target ~ 0.10)")
# Salida: |rho|_max ~ 0.10xx — AP preserva la amplitud well-prepared

El CFL_acoustic = 50 del comentario es la aritmética del demo. Con $M = 0.01$ la velocidad del sonido adimensional es $1/M = 100$. Con $dt = 0.5\,dx$, $c\,dt/dx = 100 \times 0.5 = 50$. Un paso explícito necesita ese cociente por debajo de 1. u_new = rhs estalla en el primer paso por ese 50.

$\alpha = (dt/(M\,dx))^2$ con $dt/dx = 0.5$ y $M = 0.01$ da $\alpha = 2500$. El autovalor es $1 + 2\alpha(1 - \cos(2\pi k/N))$. El modo $k = 0$ tiene autovalor 1, así que la densidad media se queda. En $N = 64$ el modo más corto es $k = 32$, y $\cos(\pi) = -1$ deja el autovalor en $1 + 4\alpha = 10001$. Una DFT divide ese modo por 10001. El pulso parte con amplitud 0.1 y anchura $\exp(-120(x-0.5)^2)$. Tras 200 pasos el objetivo del print sigue en $|\rho|_{\max} \approx 0.10$.

En $M = 10^{-4}$ la fluctuación de presión debe quedar en $\mathcal{O}(M^2) = 10^{-8}$. Roe explícito deja $\mathcal{O}(M) = 10^{-4}$ (Guillard–Viozat, 1999). Hay cuatro órdenes entre ambos.

Sustituir el solve implícito por un paso explícito (u_new = rhs) hace estallar el bucle en la primera iteración. La CFL acústica es 50 y la estabilidad explícita exige un valor menor que 1.

Tropiezos durante la reproducción#

Repetir la tabla del vórtice isentrópico de la sección 5.1 muestra reducción de orden del IMEX de cuarto orden cerca de $M = 10^{-4}$. Los RK de alto orden lo hacen a menudo en problemas rígidos (Hairer–Wanner, 1996). Los autores suben el grado polinomial de la velocidad a $r+1$. Cuesta, y sirve.

Otro tropiezo vive en las mallas DG cuadriláteras. Como mostraron Jung–Perrier (2024), los triángulos mantienen la precisión a bajo Mach y los cuadriláteros filtran ruido $\mathcal O(M)$ a la presión. El artículo reconoce el límite y lo deja para trabajo posterior (DG estable en entropía, elementos finitos compatibles). En producción con cuadriláteros aún hace falta un arreglo de bajo Mach (Thornber, 2008; Rieper, 2011).

La familia rhoPimpleFoam de OpenFOAM ya separa presión y momento (SIMPLE/PISO), pero la ruta de EOS asume gas ideal. El tratamiento de EOS cúbica del artículo se traslada a esa extensión tal como está.

Dejar o no las celdas cuadriláteras#

El Teorema 3.4 no asume la forma de la curva de la EOS. El DG cuadrilátero queda fuera de ese teorema. El ruido de presión $\mathcal{O}(M)$ de Jung–Perrier (2024) permanece aunque el integrador temporal sea AP. Si la malla no puede pasar a triángulos, Thornber (2008) o Rieper (2011) va sobre el flujo. El DG estable en entropía es la opción que el artículo dejó para después.

El 30 % de la velocidad del sonido es otra decisión. Sustituir una EOS cúbica por gas ideal puede desplazar $c$ en ese orden, y una celda con presión negativa detiene la evaluación de la EOS. Otro corte que saca el paso temporal de la velocidad del sonido está en el IMEX acústico–transporte de Peluchon. Ese corte tampoco borra el ruido de los cuadriláteros. Los 200 pasos del demo y $|\rho|_{\max} \approx 0.10$ se leyeron solo en acústica 1D periódica.