[Reseña de paper] Resolver con dos presiones, fundir en una sola — Saurel (2009) sobre relajación multifásica

Cómo el modelo de seis ecuaciones con relajación rígida de presión esquiva las trampas del modelo difuso de cinco ecuaciones.

En el laboratorio nacional de Idaho una simulación de seguridad nuclear se detuvo de golpe. En la celda donde una pequeña burbuja de gas quedaba atrapada entre sodio líquido, la velocidad del sonido de la mezcla colapsó hasta unos 25 m/s. Ése era el límite del modelo multifásico de cinco ecuaciones, atado a la hipótesis de equilibrio de presiones. En 2009 Richard Saurel y dos colegas propusieron una respuesta inesperada — dejar que las dos fases lleven presiones distintas y forzarlas a coincidir al final de cada paso temporal. Este artículo recorre por qué funciona ese algoritmo de seis ecuaciones de no equilibrio más relajación, de dónde sale la presión de interfaz $\hat{p}_I$ como media ponderada por impedancia, y qué muestran de forma distinta las dos curvas de velocidad del sonido (Wood y frozen).

La idea en una línea#

El modelo de cinco ecuaciones se rompe en la velocidad del sonido. Hay que retroceder al modelo de seis ecuaciones, resolverlo y luego fundir las dos presiones con una relajación rígida. La media por impedancia $\hat{p}_I = (Z_2 p_1 + Z_1 p_2)/(Z_1+Z_2)$ es la línea que sostiene todo el armazón.

La velocidad de Wood es un saco de arena#

El modelo de cinco ecuaciones de Kapila (2001) supone equilibrio de presiones: ambas fases comparten una sola presión dentro de cada celda. Parece limpio, pero la velocidad del sonido de mezcla se desordena. La velocidad de equilibrio $c_{eq}$ obedece la fórmula de Wood (1930)

\frac{1}{\rho c_{eq}^2} = \frac{\alpha_1}{\rho_1 c_1^2} + \frac{\alpha_2}{\rho_2 c_2^2}

con $\alpha_k$ la fracción volumétrica de la fase $k$ , $\rho_k$ la densidad y $c_k$ la velocidad del sonido pura. El problema: la curva no es monótona. Una mezcla de agua (1500 m/s) y aire (340 m/s) en una sola celda puede caer hasta unos 24 m/s — muy por debajo de cualquiera de las dos fases puras.

Conviene experimentar con la simulación de abajo. Mover los deslizadores de $\rho_k$ y $c_k$ revela cómo divergen las dos curvas.

ρ₁ = 1000 kg/m³c₁ = 1500 m/sρ₂ = 10.0 kg/m³c₂ = 340 m/scursor α₁ = 0.500log scale on y-axis

Wood (equilibrium) sound speed dips far below either pure-fluid value near intermediate α; frozen speed stays inside the [c₂, c₁] band.

En la zona de difusión numérica aparece una mezcla artificial entre los dos estados puros y las ondas acústicas se frenan bruscamente al atravesarla. Los tiempos de llegada se desfasan, el seguimiento del choque se rompe. Además, mantener $\alpha$ estrictamente positivo cerca de 0 o 1 también resulta delicado.

El desvío — dejar que las presiones discrepen#

Saurel da un paso atrás. Renuncia a la hipótesis de equilibrio y permite que cada fase lleve su propia presión $p_1$ , $p_2$ . Tomando solo el límite de relajación de velocidad sobre el modelo de siete ecuaciones de Baer–Nunziato, se obtiene el modelo de seis ecuaciones.

\partial_t \alpha_1 + u\,\partial_x \alpha_1 = \mu(p_1 - p_2)

\partial_t (\alpha_k \rho_k) + \partial_x(\alpha_k \rho_k u) = 0

\partial_t (\rho u) + \partial_x(\rho u^2 + \alpha_1 p_1 + \alpha_2 p_2) = 0

\partial_t (\alpha_k \rho_k e_k) + \partial_x(\alpha_k \rho_k e_k u) + \alpha_k p_k\,\partial_x u = \pm p_I \mu (p_1 - p_2)

Aquí $\mu$ es la tasa de relajación de presión (rígida), $p_I$ la presión de interfaz y $e_k$ la energía interna por fase. La velocidad del sonido del modelo se vuelve

c_f^2 = Y_1 c_1^2 + Y_2 c_2^2

con $Y_k = \alpha_k \rho_k / \rho$ la fracción másica. La curva queda monótona: el saco de arena de Wood desaparece.

La presión de interfaz $\hat{p}_I$ — promedio ponderado por impedancia#

El $p_I$ que aparece a la derecha de las ecuaciones de energía interna proviene del término asimétrico del modelo de siete ecuaciones tras tomar el límite de equilibrio de velocidades, y cae en la forma compacta

p_I = \frac{Z_2 p_1 + Z_1 p_2}{Z_1 + Z_2}, \quad Z_k = \rho_k c_k

donde $Z_k$ es la impedancia acústica. La fase más rígida pesa más. Intuitivamente, cuando dos fases intercambian trabajo a través de una interfaz, el lado más duro se mueve menos y el más blando se mueve más en exactamente esa proporción. Esa única línea garantiza simultáneamente la conservación de energía y la desigualdad de entropía.

El paso de relajación — una ecuación de una incógnita#

Numéricamente no se integra el límite $\mu \to \infty$ directamente. Se aplica división de operadores: primero se resuelve la parte hiperbólica y después solo la fuente de relajación. Durante la relajación, $u$ , $\rho E$ y $\alpha_k \rho_k$ permanecen congelados; solo evolucionan los volúmenes específicos $v_k = 1/\rho_k$ . Con la EOS de gas rigidizado

p_k = (\gamma_k - 1)\rho_k e_k - \gamma_k p_{\infty,k}

la integral cierra en forma analítica.

v_k(p) = v_k^0\,\frac{p^0 + \gamma_k p_{\infty,k} + (\gamma_k - 1)\hat{p}_I}{p + \gamma_k p_{\infty,k} + (\gamma_k - 1)\hat{p}_I}

La condición de saturación $\sum_k (\alpha\rho)_k v_k(p) = 1$ fija una sola presión relajada $p$ . Una bisección basta. Ése es el corazón del algoritmo.

Reinicialización — reanclar con la energía conservada#

La $p$ y $\alpha_k$ que entrega la relajación son consistentes con cada EOS de fase, pero pueden separarse de la energía de mezcla conservada porque las ecuaciones no conservativas de energía interna acumulan pequeños errores cerca de los choques. Saurel añade un retoque más. A partir de la energía total de mezcla conservada, se recalcula la presión con la EOS de mezcla

p(\rho, e, \alpha_1, \alpha_2) = \frac{\rho e - \big(\frac{\alpha_1 \gamma_1 p_{\infty,1}}{\gamma_1 - 1} + \frac{\alpha_2 \gamma_2 p_{\infty,2}}{\gamma_2 - 1}\big)}{\frac{\alpha_1}{\gamma_1 - 1} + \frac{\alpha_2}{\gamma_2 - 1}}

y luego se recalcula $e_k$ a esa presión. El sistema arranca como no conservativo, pero en los choques se comporta como uno conservativo. Coste: una evaluación adicional de EOS por celda.

Código — el momento en que dos presiones se encuentran en una celda#

Solo el paso de relajación, en unas treinta líneas de NumPy. Gas rigidizado más bisección sobre la condición de saturación.

import numpy as np
 
def stiffened_gas_volume(p, p0, v0, gamma, p_inf, p_hat_I):
    num = p0 + gamma * p_inf + (gamma - 1.0) * p_hat_I
    den = p   + gamma * p_inf + (gamma - 1.0) * p_hat_I
    return v0 * num / den
 
def saturation_residual(p, alpha_rho, p0, v0, gamma, p_inf, p_hat_I):
    total = 0.0
    for k in range(2):
        v_k = stiffened_gas_volume(p, p0[k], v0[k], gamma[k], p_inf[k], p_hat_I)
        total += alpha_rho[k] * v_k
    return total - 1.0
 
def stiff_pressure_relax(p0, v0, alpha_rho, gamma, p_inf, Z, tol=1e-9):
    p_hat_I = (Z[1] * p0[0] + Z[0] * p0[1]) / (Z[0] + Z[1])
    lo, hi = 1.0, 1.0e10
    for _ in range(80):
        mid = 0.5 * (lo + hi)
        r = saturation_residual(mid, alpha_rho, p0, v0, gamma, p_inf, p_hat_I)
        if r > 0:
            lo = mid       # presion demasiado baja -> suma de volumenes > 1
        else:
            hi = mid
        if hi - lo < tol * mid:
            break
    return 0.5 * (lo + hi), p_hat_I
 
# Agua (alpha=0.7) y aire (alpha=0.3), dos presiones distintas
gamma  = np.array([4.4, 1.4])
p_inf  = np.array([6.0e8, 0.0])
rho0   = np.array([1000.0, 1.0])
alpha0 = np.array([0.7, 0.3])
p0     = np.array([5.0e5, 1.0e5])
Z      = np.array([1500.0 * 1000.0, 340.0 * 1.0])
 
alpha_rho = alpha0 * rho0
v0        = 1.0 / rho0
 
p_eq, p_hat_I = stiff_pressure_relax(p0, v0, alpha_rho, gamma, p_inf, Z)
print(f"media de impedancia ^pI = {p_hat_I:.3e} Pa")
print(f"presion relajada     p  = {p_eq:.3e} Pa")

La impedancia del agua ( $Z_1 = 1.5\times 10^6$ ) es unas 4400 veces la del aire ( $Z_2 = 340$ ); en la práctica $\hat{p}_I$ casi coincide con la presión del agua. La presión relajada final surge de una compresión y expansión mutuas mínimas entre las dos fases.

Observación — el momento en que dos presiones se encuentran#

Una sola celda no muestra la dinámica. En una rejilla 1D con un modelo simplificado de ODE, $p_1$ y $p_2$ convergen exponencialmente hacia la media por impedancia $\hat{p}_I$ — más rápido cuanto mayor sea $\mu \cdot \Delta t$ .

Z₁ = 1.50 MPa·s/mZ₂ = 3.40 kPa·s/mμ·Δt = 0.050p₁(left) = 5.00 barp₂(left) = 1.00 barp₁(right) = 1.00 barp₂(right) = 5.00 bar

Stiff relaxation drives both p₁ and p₂ toward the impedance-weighted average p̂ᴵ = (Z₂p₁ + Z₁p₂)/(Z₁+Z₂). Increase μ·Δt for faster pull-in.

Al aumentar $Z_1$ , $\hat{p}_I$ se desplaza más hacia la presión de la fase 1; el cociente de impedancias actúa como ponderación. Con $\mu \cdot \Delta t$ cerca de 1, el equilibrio se alcanza prácticamente en un solo paso — eso significa relajación rígida. Saurel aplica este procedimiento sin condiciones al final de cada paso temporal y recupera el límite rígido del modelo de cinco ecuaciones.

Límites — la sombra de las siete ecuaciones#

Seis ecuaciones no son la respuesta universal. El equilibrio de velocidades sigue siendo una hipótesis. Medios porosos o explosivos granulares, con velocidades muy distintas entre fases, exigen volver al modelo completo de siete ecuaciones de Baer–Nunziato. Bajo choques muy intensos los errores no conservativos de energía interna no son despreciables, y entra la corrección de intercambio térmico artificial de Petitpas (sección 5 del paper). Y $\mu \to \infty$ supone una respuesta de interfaz infinitamente rápida; cuando entran tensión superficial o transferencia de masa, hace falta un modelo aparte.

Aun así, seis ecuaciones más relajación sigue siendo la base más usada hoy en simulación multifásica compresible. La integración temporal AP-IMEX (artículo del 2026-05-03) y la reconstrucción MUSCL-THINC-BVD (artículo del 2026-05-02) se apoyan en este esqueleto de seis ecuaciones.

Tres hechos para llevarse#

La trampa del modelo de cinco ecuaciones está en la velocidad del sonido. La no monotonía de la curva de Wood desvía las ondas acústicas en las interfaces.
Seis ecuaciones más relajación recupera el límite de cinco. No hay que temer al sistema no conservativo: resolver una vez, fundir una vez.
La media por impedancia $\hat{p}_I = (Z_2 p_1 + Z_1 p_2)/(Z_1+Z_2)$ es la pieza clave. Conservación de energía, desigualdad de entropía, velocidad del sonido monótona — todo cabe en una sola línea.