Darle a cada fase sus propias ecuaciones — El modelo de siete ecuaciones de Baer–Nunziato

En 1986, Mark Baer y John Nunziato querían explicar cómo una llama dentro de un propelente sólido se convierte en detonación. El explosivo era un montón de granos diminutos. El gas caliente se filtra por los huecos y quema los granos. En un mismo punto, sólido y gas coexistían a presiones distintas y velocidades distintas. Los modelos de mezcla de la época, que suponen una sola presión y una sola velocidad, no podían capturar ese instante. Por eso ambos construyeron un modelo que da a cada fase su propio conjunto completo de leyes de conservación. Este artículo recorre la estructura de ese modelo de siete ecuaciones: por qué hacen falta siete ecuaciones, por qué los términos no conservativos son un dolor de cabeza y cómo la relajación (relaxation) engendra modelos más simples. Al final integramos directamente el operador de relajación de presión y observamos cómo las dos fases marchan hacia el equilibrio.

Siete ecuaciones nacidas de granos de pólvora#

Los métodos de interfaz difusa (diffuse interface, que difuminan la frontera entre dos fases sobre un espesor finito) tienen ramas. La forma más simple, de cuatro ecuaciones, supone que ambas fases comparten presión, velocidad y temperatura por igual. Baer–Nunziato (BN) se sitúa en el extremo opuesto. No supone equilibrio alguno. Cada fase lleva su propia densidad, velocidad y presión. A eso se suma una ecuación de transporte más para la fracción de volumen (volume fraction) $\alpha$, la fracción de volumen que ocupa una de las fases.

¿Por qué importa esta generalidad? Admitir el no equilibrio hace que el modelo sea incondicionalmente hiperbólico (hyperbolic). El modelo de seis ecuaciones, que impone una sola presión, pierde la hiperbolicidad bajo ciertas condiciones y produce soluciones no físicas, como una velocidad del sonido al cuadrado negativa. BN evita esa trampa. Pero la paga. El recuento sube a siete y se cuelan términos que no pueden escribirse en forma conservativa estándar.

Al desplegar las siete ecuaciones#

Etiquetemos las dos fases como $k = 1, 2$. Primero, la ecuación de transporte de la fracción de volumen.

El lado izquierdo transporta la fracción de volumen a la velocidad de interfaz $u_I$; el derecho es la relajación de presión, con $\mu$ la tasa de relajación (cuanto mayor, más rápido el equilibrio). Las seis restantes son los balances de masa, cantidad de movimiento y energía de cada fase.

Aquí $\rho_k$ es la densidad fásica, $u_k$ la velocidad fásica, $p_k$ la presión fásica, $E_k = \rho_k e_k + \tfrac12 \rho_k u_k^2$ la densidad de energía total de la fase y $k^*$ la fase compañera. $\lambda$ es la tasa de relajación de velocidad. Cada fase se cierra con su propia ecuación de estado. Sea para propelente o cavitación, una elección habitual es el gas rigidizado (stiffened gas).

$\gamma_k$ es la razón de calores específicos y $p_{\infty,k}$ una constante de presión que representa la repulsión intermolecular. Un líquido casi incompresible como el agua tiene un $p_\infty$ grande, de varios GPa.

Términos no conservativos y variables de interfaz — la parte verdaderamente difícil#

Los lados derechos de las ecuaciones de cantidad de movimiento y energía llevan $\partial \alpha_k / \partial x$. Este término no se pliega en una divergencia. Es decir, el sistema no puede escribirse como una ley de conservación $\partial_t U + \partial_x F(U) = 0$. Los esquemas de volúmenes finitos tipo Godunov tropiezan ante esos términos no conservativos (non-conservative), porque resulta ambiguo por qué valor de $\alpha$ multiplicar cuando salta a través de una discontinuidad. Discretízalo mal y brotan oscilaciones espurias a través de la interfaz incluso en un campo uniforme de presión y velocidad.

La piedra de toque que filtra esta trampa es el criterio de Abgrall. Si la presión y la velocidad arrancan uniformes, esa uniformidad debe preservarse sin importar cómo se distribuya la fracción de volumen. Cuando dos fases fluyen a la misma presión y velocidad, la interfaz solo debería ser transportada, no generar nuevas ondas. La discretización de los términos no conservativos se diseña para satisfacer exactamente esta condición.

Las variables de interfaz $u_I$ y $p_I$ son también producto del modelado. La elección más común, de Saurel–Abgrall, es una velocidad ponderada por masa y una presión ponderada por volumen.

Esta elección hace simétrico el modelo, de modo que las mismas ecuaciones y el mismo método numérico se aplican en cada celda.

Relajación: el camino al equilibrio#

Las dos fuentes del lado derecho, $\mu(p_1-p_2)$ y $\lambda(u_{k^*}-u_k)$, son los operadores de relajación. Sube $\mu$ y $\lambda$ y las dos fases son arrastradas a una presión y velocidad comunes casi al instante. Ese es el encanto de BN: el límite de tasa de relajación infinita es exactamente un modelo de equilibrio más simple.

Prueba la simulación de abajo. Las dos fases parten de presiones y velocidades distintas, y las tasas de relajación $\mu$ y $\lambda$ las arrastran hacia el equilibrio.

Baja las tasas de relajación a cero y las dos curvas nunca se encuentran. Ese es el estado de no equilibrio que el modelo de siete ecuaciones permite. Cuanto mayor sea $\mu$, más abrupto será el colapso de las curvas de presión hacia un punto. Fíjate también en que la presión y la velocidad se relajan de forma independiente entre sí.

La jerarquía: 7 → 6 → 5 → 4#

Supón que la relajación es instantánea —impón el equilibrio— y una ecuación colapsa en una restricción algebraica y desaparece. Fija la velocidad al equilibrio y obtienes seis ecuaciones; añade la presión y obtienes las cinco de Kapila; añade la temperatura y obtienes cuatro. BN se sitúa en lo alto de esta jerarquía.

Abajo, pulsa cada modelo para comparar qué se supone en equilibrio y qué se pierde.

El modelo de cinco ecuaciones de Kapila es el más usado para la captura de interfaz, pero su velocidad del sonido de mezcla es no monótona respecto de la fracción de volumen (la curva de Wood), lo que vuelve rígida (stiff) la fuente de la fracción de volumen. BN es numéricamente más flexible porque despega esa rigidez como una fuente de relajación y la trata por separado.

Python — Integrando el operador de relajación de presión#

Aislemos el corazón de la relajación. Las masas parciales y las energías internas de las dos fases se conservan, y solo se mueve la fracción de volumen $\alpha$, según $d\alpha/dt = \mu(p_1 - p_2)$. Integramos cómo varían las presiones del gas rigidizado con $\alpha$ y alcanzan el equilibrio.

import numpy as np
 
# Operador de "relajación de presión" para dos fases de gas rigidizado.
# La masa parcial m_k y la energía interna E_k de cada fase se conservan;
# solo la fraccion de volumen alpha evoluciona por  da/dt = mu (p1 - p2).
 
def phase_pressures(alpha, mass, energy, gamma, p_inf):
    frac = np.array([alpha, 1.0 - alpha])
    rho  = mass / frac          # densidad fasica = masa parcial / fraccion de volumen
    e    = energy / mass        # energia interna especifica
    return (gamma - 1.0) * rho * e - gamma * p_inf
 
def relax_pressure(alpha0, mass, energy, gamma, p_inf,
                   mu=4.0e-8, dt=1.0e-3, nmax=200000, tol=1.0):
    alpha, hist = alpha0, []
    for n in range(nmax):
        p = phase_pressures(alpha, mass, energy, gamma, p_inf)
        hist.append((n * dt, alpha, p[0], p[1]))
        gap = p[0] - p[1]
        if abs(gap) < tol:                       # termina cuando se anula la brecha
            break
        alpha += dt * mu * gap                    # da/dt = mu (p1 - p2)
        alpha = min(max(alpha, 1e-4), 1.0 - 1e-4)
    return alpha, np.array(hist)
 
# Fase tipo agua (1, gas rigidizado) + fase tipo aire (2)
gamma  = np.array([4.4, 1.4])
p_inf  = np.array([6.0e8, 0.0])      # Pa
mass   = np.array([480.0, 0.55])     # masa parcial  (kg/m^3)
energy = np.array([3.93e8, 1.10e6])  # energia interna parcial (J/m^3)
 
alpha_eq, hist = relax_pressure(0.50, mass, energy, gamma, p_inf)
p0  = phase_pressures(0.50, mass, energy, gamma, p_inf)
peq = phase_pressures(alpha_eq, mass, energy, gamma, p_inf)
print(f"inicial: alpha=0.500  p1={p0[0]/1e5:8.3f} bar  p2={p0[1]/1e5:7.3f} bar")
print(f"equil.:  alpha={alpha_eq:.3f}  p1={peq[0]/1e5:7.3f} bar  p2={peq[1]/1e5:7.3f} bar")
print(f"{len(hist)} pasos,  brecha final de presion {abs(peq[0]-peq[1]):.3e} Pa")

La ejecución imprime esto.

inicial: alpha=0.500  p1= 324.000 bar  p2=  8.800 bar
equil.:  alpha=0.506  p1=  8.906 bar  p2=  8.906 bar
75 pasos,  brecha final de presion 9.012e-01 Pa

Lo llamativo es que la fracción de volumen se movió solo un 0,6%, de 0,500 a 0,506, y sin embargo la presión del líquido se desplomó de 324 bar a 8,9 bar. El agua modelada como gas rigidizado es casi incompresible, así que un cambio mínimo de volumen impulsa un cambio enorme de presión. Por eso también el equilibrio relajado se asienta cerca de la presión del gas: el líquido rígido puede igualar presiones apenas expandiéndose.

Vale la pena recordar#

• Siete ecuaciones = no equilibrio. Dar a cada fase su propia presión, velocidad y energía compra hiperbolicidad incondicional.
• Los términos no conservativos son la dificultad central. $p_I \partial_x \alpha$ no se pliega en forma conservativa. El criterio de Abgrall (preservar presión y velocidad uniformes) es la piedra de toque de la discretización.
• La relajación construye la jerarquía. El límite $\mu, \lambda \to \infty$ produce los modelos de 6 → 5 → 4 ecuaciones. BN se sitúa en lo alto y da el panorama más general.