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

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 de cinco ecuaciones atado al equilibrio de presiones. En 2009 Saurel y sus colegas dieron un paso atrás. Las dos presiones se resuelven aparte y solo se pegan al final del paso.

Al mover los deslizadores de $\rho_k$ y $c_k$ cambia la separación entre la curva de Wood y la congelada.

Cerca del agua a 1500 m/s y del aire a 340 m/s, solo la rama de Wood cae al suelo. La rama congelada se queda entre las dos velocidades puras.

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)

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.

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.

Al meter la celda del código en la fórmula de Wood reaparece el 24 m/s del texto. Con $\alpha_1 = 0.7$, $\rho_1 = 1000$, $c_1 = 1500$ sale $\rho_1 c_1^2 = 2.25\times 10^9$. En el aire, $\alpha_2 = 0.3$, $\rho_2 = 1$, $c_2 = 340$, así que $\rho_2 c_2^2 = 1.156\times 10^5$. El término del aire a la derecha ronda $2.6\times 10^{-6}$; el del agua, $3\times 10^{-10}$, casi no se ve. La densidad de mezcla es $0.7\times 1000 + 0.3 = 700.3$. Cuando $1/(\rho c_{eq}^2)$ lo come el término del aire, $c_{eq}$ queda cerca de 24 m/s. La fórmula congelada en la misma celda, $c_f^2 = Y_1 c_1^2 + Y_2 c_2^2$, va al revés. Casi toda la masa es agua, así que $c_f$ se pega a 1500 m/s. En una celda las dos velocidades se separan unas sesenta veces. Ese hundimiento es lo que sigue una velocidad del sonido con agua y aire a medias.

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.

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

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

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

la integral cierra en forma analítica.

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

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$). La fórmula es $\hat{p}_I = (Z_2 p_1 + Z_1 p_2)/(Z_1 + Z_2)$. El peso grande cae sobre la presión pequeña. Con $p_1 = 5\times 10^5$ y $p_2 = 1\times 10^5$, en el numerador $Z_1 p_2$ vale $1.5\times 10^{11}$ y $Z_2 p_1$ se queda en $1.7\times 10^8$. Así $\hat{p}_I$ no se sienta junto a la presión del agua. Queda al lado de la presión del aire $1\times 10^5$, cerca de $1.001\times 10^5$ Pa. La fase rígida no se queda con la interfaz. La presión común tras la relajación la busca después la bisección de la saturación.

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$.

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 y la reconstrucción MUSCL-THINC-BVD corren sobre este esqueleto de seis ecuaciones.

Dejar el modelo de cinco ecuaciones en esa celda#

Entregar al modelo de cinco ecuaciones una celda cuya velocidad de Wood cayó a 24 m/s retrasa las llegadas acústicas y mueve la posición del choque. Esa celda retrocede a seis ecuaciones, lleva dos presiones y las pega con la media de impedancia solo al final del paso.

Si además discrepan las velocidades, el compromiso se acaba. Medios porosos y explosivos granulares vuelven a siete ecuaciones. Cerrar la interfaz con aún menos ecuaciones es otra pregunta, la de un equilibrio de cuatro ecuaciones.