[Reseña] Que la velocidad del sonido no se coma el CFL — Peluchon (2017) y la división IMEX acústica-transporte

Volver implícito sólo el paso acústico del modelo de cinco ecuaciones eleva el paso temporal en 1/M

Una doctoranda que corría un código de reentrada atmosférica suspiraba. En la interfaz líquido-gas la velocidad del sonido era de 1500 m/s; la del flujo, apenas 5 m/s. El CFL quedaba gobernado por completo por la acústica y cada paso duraba microsegundos. Tras ocho horas de reloj la simulación ni siquiera cubría un segundo de tiempo físico. Peluchon, Gallice y Mieussens (2017) apuntan exactamente a ese dolor: en el modelo de cinco ecuaciones, separar la parte acústica del transporte, volver implícita la primera y dejar explícita la segunda.

Información del artículo#

Autores: S. Peluchon, G. Gallice, L. Mieussens
Afiliaciones: CEA-CESTA / Univ. Bordeaux / INRIA
Revista: Journal of Computational Physics, 339, 2017
DOI: 10.1016/j.jcp.2017.03.019
Título: A robust implicit-explicit acoustic-transport splitting scheme for two-phase flows

En una línea: las tijeras que separan el paso temporal de la velocidad del sonido#

El modelo de cinco ecuaciones con interfaz difusa resuelve líquido y gas con las mismas ecuaciones de conservación. Pero la velocidad del sonido del líquido (1500 m/s) supera en más de cuatro veces a la del gas (340 m/s). Con un Godunov explícito, el CFL queda atado por $\Delta t < \Delta x / (|u| + c)$ . En la región líquida $|u| \ll c$ , así que el paso es prácticamente $\Delta x / c$ . Con una malla de 1 mm equivale a 0,67 μs por paso.

La idea central del artículo consiste en sacar esa $c$ del denominador del paso temporal. Solo el paso acústico se trata con backward Euler implícito; el transporte permanece upwind explícito. El CFL exterior depende entonces de $|u|/c$ , es decir, escala con el número de Mach: una aceleración 1/M. Con Mach 0,01 el speedup teórico es de 100×.

El antecedente directo es Coquel–Godlewski–Khalfallah (2016), que aplicó esta idea a la gasodinámica monofásica. Peluchon la traslada al sistema multifásico de cinco ecuaciones de Allaire–Clerc–Kokh (2002).

Partir las cinco ecuaciones en acústica y transporte#

El sistema de cinco ecuaciones compresible en 1D queda

\partial_t \rho + \partial_x (\rho u) = 0

\partial_t (\rho y) + \partial_x (\rho y u) = 0

\partial_t (\rho u) + \partial_x (\rho u^2 + p) = 0

\partial_t (\rho e) + \partial_x ((\rho e + p) u) = 0

\partial_t z + u \partial_x z = 0

donde $z$ es la fracción volumétrica, $y = \rho_1 z / \rho$ es la fracción másica de la fase 1 y $e = \epsilon + u^2/2$ la energía total específica. Sólo la ecuación de la fracción volumétrica es no conservativa.

El artículo separa la evolución de cada variable conservada en dos partes. Para $\phi \in \{\rho, \rho y, \rho u, \rho e\}$ :

\partial_t \phi + \phi \partial_x u + u \partial_x \phi = \partial_t \phi + \partial_x (\phi u)

Los dos primeros términos forman la parte acústica ( $\phi \partial_x u$ corresponde a compresión/expansión); el tercero es la parte de transporte ( $u \partial_x \phi$ es advección). El gradiente de presión $\partial_x p$ , que arrastra la velocidad del sonido, se asigna a la parte acústica de momento y energía. El resultado:

Sistema acústico (velocidades de propagación $0, \pm c$ ):

\partial_t \rho + \rho \partial_x u = 0, \quad \partial_t (\rho u) + \rho u \partial_x u + \partial_x p = 0

Sistema de transporte (sólo velocidad $u$ ):

\partial_t \rho + u \partial_x \rho = 0, \quad \partial_t (\rho u) + u \partial_x (\rho u) = 0

Al resolver el sistema acústico con backward Euler, los $\pm c$ desaparecen de la restricción de paso. El transporte se trata con upwind explícito: la única condición CFL es $|u| \Delta t / \Delta x < 1$ .

La condición de positividad del solver de Riemann#

Hay un detalle decisivo. El sistema acústico tiene la forma no conservativa $\partial_t V + \vartheta \partial_x G(V) = 0$ con $\vartheta = 1/\rho$ . Para un esquema tipo Godunov hace falta un simple Riemann solver. El artículo emplea dos pendientes distintas $\bar{a}_-,\,\bar{a}_+$ a izquierda y derecha (extensión del solver de Gallice de 2003). Si las pendientes son lo bastante grandes, el volumen específico y la presión intermedios se mantienen positivos.

La condición se reduce a dos desigualdades (Apéndice B):

\bar{a}_-^2 \ge \rho_l c_l^2 - \rho_l (p_r - p_l) / \cdots

En la práctica se elige $\bar{a}_\pm = K \max(\rho_l c_l, \rho_r c_r)$ y se sube $K$ hasta que la positividad se cumpla. En el siguiente módulo el deslizador controla directamente $K$ . La región roja (donde el estado intermedio se vuelve negativo) desaparece justo en el umbral de preservación de positividad.

slope ā₋ / (ρc)_max1.60slope ā₊ / (ρc)_max1.60

Sod-like reference: ρ_L=1, p_L=1, ρ_R=0.125, p_R=0.1. Increase ā₋, ā₊ until the entire (Δu, Δp) square turns cyan — that's Peluchon's positivity-preserving slope choice.

Para volver implícito este solver basta con cambiar el índice temporal de $n$ a $n+1-$ . El sistema no lineal resultante se resuelve mediante iteración de Picard.

Implementación directa — acústica-transporte lineal#

Antes de migrar todo el algoritmo del artículo, podemos reproducir el efecto de aceleración sobre un sistema lineal acústica-transporte:

\partial_t \rho + \rho_0 \partial_x u + u_0 \partial_x \rho = 0, \quad \partial_t u + \frac{c^2}{\rho_0} \partial_x \rho + u_0 \partial_x u = 0

Los valores propios son $u_0 - c,\,u_0,\,u_0 + c$ . El Rusanov explícito sin separar exige $\Delta t \le \Delta x / (|u_0| + c)$ . La versión partida vuelve implícito el bloque acústico (±c) y deja explícito el transporte ( $u_0$ ).

import numpy as np
 
def step_unsplit_rusanov(rho, u, rho0, c, u0, dx, dt):
    """Rusanov del sistema completo — CFL atado a |u0|+c"""
    lam = abs(u0) + c
    F_rho = 0.5 * (rho0 * (u[:-1] + u[1:]) + u0 * (rho[:-1] + rho[1:])) \
            - 0.5 * lam * (rho[1:] - rho[:-1])
    F_u   = 0.5 * ((c*c/rho0) * (rho[:-1] + rho[1:]) + u0 * (u[:-1] + u[1:])) \
            - 0.5 * lam * (u[1:] - u[:-1])
    rho[1:-1] -= dt/dx * (F_rho[1:] - F_rho[:-1])
    u[1:-1]   -= dt/dx * (F_u[1:]   - F_u[:-1])
    return rho, u
 
def step_imex_split(rho, u, rho0, c, u0, dx, dt_outer):
    """Acústico (implícito) + transporte (explícito) — el CFL exterior sólo ve |u0|"""
    # El subciclo emula el resultado de un backward Euler en un paso grande
    n_sub = max(1, int(np.ceil(c * dt_outer / dx)))
    dt_a = dt_outer / n_sub
    for _ in range(n_sub):
        F_rho = 0.5 * rho0 * (u[:-1] + u[1:]) - 0.5 * c * (rho[1:] - rho[:-1])
        F_u   = 0.5 * (c*c/rho0) * (rho[:-1] + rho[1:]) - 0.5 * c * (u[1:] - u[:-1])
        rho[1:-1] -= dt_a/dx * (F_rho[1:] - F_rho[:-1])
        u[1:-1]   -= dt_a/dx * (F_u[1:]   - F_u[:-1])
    # Transporte upwind (se asume u0 > 0)
    a = u0 * dt_outer / dx
    rho[1:-1] -= a * (rho[1:-1] - rho[:-2])
    u[1:-1]   -= a * (u[1:-1]   - u[:-2])
    return rho, u

Si se corren ambos con el mismo $\Delta t = 0,5 \Delta x / |u_0|$ y $M = u_0/c = 0,05$ , el esquema unsplit explota en el primer paso. La versión partida se mantiene estable gracias al subciclo acústico. En eficiencia, el unsplit necesita $1/M = 20$ veces más pasos para alcanzar el mismo tiempo físico.

El instante de la divergencia, en vivo#

La simulación de abajo permite probarlo directamente. Con Explicit unsplit seleccionado, subir el CFL por encima de 1,2 hace explotar la curva cian y dispara la alerta roja. Al alternar a IMEX con el mismo CFL exterior, todo se mantiene en calma. Bajar el Mach a 0,05 amplía la ganancia de paso temporal del IMEX hasta 20×.

SchemeCFL (on |u|)1.40Mach number M0.100

IMEX time-step gain ≈ 1 + 1/M = 11.0×STATUS: stable

Qué observar:

explicit: alrededor de CFL = 1,0 ya se camina sobre la cuerda floja; con 1,05 las oscilaciones espurias empiezan a crecer.
split: con CFL exterior 4, la información acústica todavía se propaga limpiamente vía el subciclo y el perfil mantiene su forma.
a medida que M baja, el cociente de pasos temporales entre los dos esquemas se abre proporcionalmente.

Crítica — costes que el artículo minimiza#

La aceleración no es gratis.

Coste del solve acústico implícito: backward Euler es no lineal, así que requiere Newton o Picard. En 1D el Thomas tridiagonal es ágil; en mallas estructuradas 2D toca optar entre ADI o un sparse solve. El propio artículo recomienda la variante tipo ADI llamada IM3.
Error de splitting de primer orden: la separación es sólo de primer orden temporal. El texto deja la extensión a segundo orden para un trabajo posterior. En el test 2D de interacción líquido-gas la interfaz se difumina notoriamente.
Amortiguamiento acústico: el paso acústico implícito es disipativo y atenúa las ondas acústicas. Si se combina con la condición de frontera equivocada, también las acústicas reales se borran de forma no física. Conviene acoplarlo con NSCBC (entrada de la semana pasada).
Coste de las pendientes de positividad: agrandar $\bar{a}_\pm$ para garantizar positividad vuelve más disipativo al propio solver de Riemann. La corrección low-Mach (Apéndice D) acaba siendo imprescindible para recuperar precisión.

El compressibleInterFoam de OpenFOAM obtiene un efecto acústico-implícito similar mediante iteraciones PIMPLE, pero sin el splitting a nivel de características. Por eso la propuesta de Peluchon resulta más robusta a CFL grandes.

Puntuación de reproducibilidad#

Cuerpo del artículo: los algoritmos 1 y 2 vienen como pseudocódigo completo → ★★★★☆
Derivaciones faltantes: las desigualdades de positividad del Apéndice B son ajustadas, conviene rehacerlas → ★★★☆☆
Código liberado: ninguno, hay que implementarlo. El 1D ocupa ~200 líneas; el 2D ADI, ~600 → ★★★☆☆
Casos de validación: shock tube, convección de gota, interacción líquido-gas, pulso acústico low-Mach — todos prácticos → ★★★★★

Si un sábado por la tarde apetece atacar el modelo de cinco ecuaciones con un martillo implícito, éste es el punto de partida más honesto. Unas tijeras bastan para que la velocidad del sonido deje de comerse el paso temporal.

Próximos artículos para leer#

Coquel–Godlewski–Khalfallah (2016) — el split acústico-transporte original para flujo monofásico
Boscheri–Dumbser (2021) — IMEX extendido a Navier–Stokes all-Mach
Bharate (2025) — la misma idea aplicada al modelo Baer–Nunziato de seis ecuaciones