[Reseña] Segundo orden en una sola etapa — un solver GRP para el modelo multifásico de Kapila

Cómo un solver GRP logra precisión temporal de segundo orden y doma el rígido término fuente de fracción de volumen a la vez

Nos enseñaron que la precisión de segundo orden exige recorrer varias etapas de Runge-Kutta. Este artículo llega allí en una sola etapa. El truco es un solver GRP (Generalized Riemann Problem) que en la interfaz devuelve no solo el valor sino también su derivada temporal. Aquí confirmamos sobre un problema escalar por qué GRP es de segundo orden de un tiro, y luego reproducimos en código cómo esa derivada temporal doma el rígido término fuente del modelo multifásico de Kapila.

Artículo: Tuowei Chen, Zhifang Du, Generalized Riemann Problem Method for the Kapila Model of Compressible Multiphase Flows, arXiv:2310.08241 (2023).

Segundo orden en una etapa — el problema que resuelve GRP#

Un paso de volúmenes finitos lo decide el flujo de interfaz $\mathbf{F}_{i+1/2}$ . El método de Godunov resuelve un problema de Riemann en la interfaz y obtiene un valor. Ese valor es constante en el tiempo. Por eso, para alcanzar el segundo orden temporal hay que repetir el mismo operador espacial varias veces, al estilo Runge-Kutta.

GRP cambia la premisa. Evalúa el flujo en el nivel temporal intermedio $t^{n+1/2}$ .

\mathbf{U}(\mathbf{x}_{i+1/2}, t^{n+1/2}) = \mathbf{U}^{n,*}_{i+1/2} + \frac{\Delta t}{2}\left(\frac{\partial \mathbf{U}}{\partial t}\right)^{n,*}_{i+1/2}

Aquí $\mathbf{U}^{n,*}$ es la solución de Riemann (valor de interfaz) y $(\partial \mathbf{U}/\partial t)^{n,*}$ es la derivada temporal instantánea en esa interfaz. Combinar ambos mediante una expansión de Taylor empuja el valor de interfaz $\Delta t/2$ hacia el futuro. Una sola evaluación, segundo orden en el tiempo.

De constante a tramos a lineal a tramos — Riemann frente a GRP#

Un problema de Riemann tiene datos constantes a tramos a cada lado de la interfaz. Un problema GRP tiene datos lineales a tramos (suaves a tramos). Las ondas que emanan de la discontinuidad inicial ya no viajan en línea recta, y el solver devuelve la derivada temporal junto con la solución de Riemann.

Construyamos intuición con la advección escalar $u_t + a u_x = 0$ . Con pendiente de celda $\sigma_i$ , la solución de Riemann en la interfaz para $a>0$ es el valor extrapolado desde la celda izquierda.

U^{*}_{i+1/2} = u_i + \sigma_i \frac{\Delta x}{2}

Aquí $u_i$ es el promedio de celda, $\sigma_i$ la pendiente interna y $\Delta x$ el ancho de celda. Conviene manipularlo en la simulación de abajo.

advection speed a: 1.00left slope uₓ: 0.50CFL: 0.80GRP (off = Godunov)

Al desmarcar la casilla GRP se obtiene Godunov: el valor de interfaz (punto cian) queda congelado mientras pasa el tiempo. Al volver a marcarla, el valor fluye a lo largo de la característica (punto ámbar); el valor que alimenta al flujo es el de $\Delta t/2$ (anillo ámbar). Con la pendiente $u_x$ en cero, GRP se reduce exactamente a Godunov.

El procedimiento de Lax-Wendroff — convertir una derivada temporal en espacial#

¿Cómo obtenemos esa derivada temporal? La respuesta vive dentro de la propia ecuación de gobierno. Ese es el corazón del procedimiento de Lax-Wendroff.

\frac{\partial u}{\partial t} = -a \frac{\partial u}{\partial x} = -a\,\sigma_i

La derivada temporal de la izquierda se convierte en la derivada espacial de la derecha. La pendiente espacial $\sigma_i$ ya la conocemos de la reconstrucción. Así, el valor de interfaz a medio paso es

U^{n+1/2}_{i+1/2} = U^{*}_{i+1/2} - a\,\sigma_i \frac{\Delta t}{2}

y el flujo es $F = a\,U^{n+1/2}$ . Para un sistema (Euler, Kapila) esto se generaliza a $\partial u/\partial t = -A\,\partial u/\partial x$ con el jacobiano $A$ . El artículo usa este procedimiento para linealizar las ondas acústicas y resolver el GRP de forma aproximada.

Python — midiendo el orden de convergencia de un solver GRP escalar#

Advectamos una condición inicial suave $u_0=\sin(2\pi x)$ durante un período y comparamos el orden de convergencia del error $L_1$ de Godunov y GRP. El problema de juguete es advección lineal en un dominio periódico.

import numpy as np
 
# Adveccion lineal 1D  u_t + a u_x = 0,  dominio periodico [0,1]
# Godunov de primer orden vs GRP de segundo orden en una etapa (Lax-Wendroff)
 
def cell_gradient(u, dx):
    # pendiente central (sin limitar) para un orden limpio en solucion suave
    return (np.roll(u, -1) - np.roll(u, 1)) / (2.0 * dx)
 
def godunov_step(u, a, dx, dt):
    # constante a tramos: valor de interfaz aguas arriba (a > 0)
    flux = a * u                       # F_{i+1/2} = a u_i
    return u - dt / dx * (flux - np.roll(flux, 1))
 
def grp_step(u, a, dx, dt):
    g = cell_gradient(u, dx)           # pendiente espacial sigma_i
    u_star = u + g * dx / 2.0          # solucion de Riemann en interfaz U*
    u_half = u_star - a * g * dt / 2.0 # medio paso de Lax-Wendroff U^{n+1/2}
    flux = a * u_half                  # F_{i+1/2}^{n+1/2}
    return u - dt / dx * (flux - np.roll(flux, 1))
 
def l1_error(num, exact, dx):
    return np.sum(np.abs(num - exact)) * dx
 
def convergence_study(a=1.0, cfl=0.4, t_end=1.0):
    print(f"{'N':>5} {'Godunov':>16} {'GRP':>16}")
    prev = {}
    for N in [40, 80, 160, 320]:
        dx = 1.0 / N
        x = (np.arange(N) + 0.5) * dx
        u0 = np.sin(2 * np.pi * x)     # condicion inicial suave
        dt = cfl * dx / abs(a)
        nsteps = int(round(t_end / dt))
        dt = t_end / nsteps            # aterrizar exactamente en t_end
        ug = u0.copy(); ur = u0.copy()
        for _ in range(nsteps):
            ug = godunov_step(ug, a, dx, dt)
            ur = grp_step(ur, a, dx, dt)
        exact = np.sin(2 * np.pi * (x - a * t_end))
        eg = l1_error(ug, exact, dx)
        er = l1_error(ur, exact, dx)
        og = f"{np.log2(prev['g']/eg):.2f}" if 'g' in prev else "  - "
        orr = f"{np.log2(prev['r']/er):.2f}" if 'r' in prev else "  - "
        print(f"{N:>5} {eg:.2e}({og}) {er:.2e}({orr})")
        prev = {'g': eg, 'r': er}
 
convergence_study()

La salida es:

    N          Godunov              GRP
   40 1.63e-01(  - ) 1.31e-03(  - )
   80 8.76e-02(0.90) 2.69e-04(2.28)
  160 4.54e-02(0.95) 6.32e-05(2.09)
  320 2.31e-02(0.97) 1.55e-05(2.03)

Godunov ronda el primer orden (0.9–0.97) mientras GRP converge cerca de 2.0. No añadimos ninguna etapa extra: solo la pendiente y la derivada temporal. En la misma malla el error es cien veces menor.

El rígido término fuente de fracción de volumen del modelo de Kapila#

Ahora el plato fuerte. La ecuación de fracción de volumen del modelo de cinco ecuaciones de Kapila (una reducción del modelo de siete ecuaciones BN) carga con un término fuente rígido.

\frac{\partial \alpha_1}{\partial t} + \mathbf{u}\cdot\nabla\alpha_1 = \alpha_1 K\,\nabla\cdot\mathbf{u}

Aquí $\alpha_1$ es la fracción de volumen de la fase 1 y $K$ es el coeficiente de compactación/expansión que sale de la velocidad del sonido de Wood. El problema es que esta fuente puede empujar $\alpha_1$ fuera de $[0,1]$ . Una integración explícita rebasa la cota y diverge en cuanto crece la rigidez.

Aquí brilla la derivada temporal de GRP. Como conocemos el valor de interfaz en el nuevo nivel temporal $t^{n+1}$ , la fórmula de Gauss-Green permite estimar el promedio de celda de $\nabla\cdot\mathbf{u}$ en ese nuevo tiempo. Entonces podemos integrar la fuente con Crank-Nicolson (semi-implícito) e incluso fijar una restricción de paso temporal que preserva la cota de la fracción de volumen. Comparemos la relajación rígida con un modelo simple.

stiffness λΔt: 1.60α₀ (initial): 0.85α_eq (target): 0.10

Sube $\lambda\Delta t$ (rigidez × paso temporal) por encima de 1.6 y el Euler explícito (rojo) rompe la banda $[0,1]$ y oscila. Bajo la misma condición, Crank-Nicolson (cian) se mantiene dentro de la banda. Lleva la rigidez hasta 3 y el cian sigue convergiendo monótonamente al equilibrio. Eso es lo que el artículo llama "preservación de la cota de fracción de volumen para cualquier $\Delta t$ ".

Lo que resultó dudoso al reproducirlo#

El GRP escalar salió limpiamente de segundo orden, pero el sistema real de Kapila no es tan dócil. Primero, una pendiente central sin limitar oscila en los choques. El segundo orden apareció porque la prueba era suave; en discontinuidades es obligatorio un limitador tipo minmod, y entonces el orden cae localmente a primero.

Segundo, la velocidad del sonido de Wood es no monótona en la fracción de volumen. Esa es la fuente de las oscilaciones del número de Mach en la zona de difusión numérica de una interfaz. La derivada temporal de GRP captura bien la compactación y la expansión, pero no elimina esa no monotonicidad en sí. El artículo también la esquiva con robustez en vez de afrontarla de frente.

Tercero, los códigos de producción (OpenFOAM y similares) se construyen en su mayoría sobre división de Strang más Runge-Kutta, así que un GRP de una sola etapa es difícil de insertar. Para disfrutar la ventaja de GRP hay que reescribir la propia función de flujo en forma Lax-Wendroff. El costo de portarlo no es pequeño.

Qué leer a continuación#

La casa del GRP es el GRP de dinámica de gases de Ben-Artzi y Falcovitz. Para ver de dónde viene la aproximación acústica linealizada de este artículo, ese es el lugar al que remontarse. En el lado multifásico, el método de relajación-proyección de Saurel (la entrada de este blog del 2026-05-09) resuelve la misma rigidez de otra manera. Pon ambos lado a lado y la bifurcación de diseño se afila: ¿separar la fuente, o disolverla dentro del flujo?

El punto es uno solo. Si calculas la derivada temporal junto al valor en la interfaz, tanto la precisión temporal de segundo orden como una fuente rígida pueden manejarse en una sola evaluación.