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

En 1986 Baer y Nunziato querían el instante en que una llama dentro de un propulsante sólido pasa a detonación. La pólvora era un montón de granos finos. El gas caliente se filtraba por los huecos y los quemaba. En un mismo punto, sólido y gas estaban a presiones distintas y a velocidades distintas. Un modelo de mezcla con una presión y una velocidad no cabía en ese instante. Así que cada fase recibió sus propias leyes de conservación.

Poner las tasas de relajación μ y λ en cero impide que las dos curvas de presión se encuentren, y subirlas las arrastra hacia un solo punto.

Conviene mirar si la presión y la velocidad se juntan por separado. Al subir un solo deslizador, la otra curva tiene que quedarse en su sitio.

Siete ecuaciones nacidas de granos de pólvora#

Los modelos de interfaz difusa, que dan a la frontera entre dos fases un espesor finito, tienen ramas. El modelo más simple, de cuatro ecuaciones, supone presión, velocidad y temperatura iguales. Baer–Nunziato está en el otro extremo. No se supone equilibrio. Cada fase guarda su densidad, su velocidad y su presión. Una ecuación más transporta la fracción de volumen $\alpha$, la parte del volumen que ocupa una fase.

Esa generalidad protege la hiperbolicidad. Permitir el no equilibrio vuelve el modelo incondicionalmente hiperbólico. Un modelo de seis ecuaciones que fuerza una sola presión puede perder la hiperbolicidad y dar un cuadrado de la velocidad del sonido negativo. BN evita esa trampa. El precio son siete ecuaciones, más términos que no entran en la forma de conservación estándar.

Al desplegar las siete ecuaciones#

Las fases se escriben $k=1,2$. Primero va la fracción de volumen.

El lado izquierdo lleva la fracción de volumen a la velocidad de interfaz $u_I$. El lado derecho es la relajación de presión. Un $\mu$ mayor significa un equilibrio más rápido. Las otras seis ecuaciones son masa, momento y energía de cada fase.

$\rho_k$ es la densidad de fase, $u_k$ la velocidad de fase, $p_k$ la presión de fase. $E_k=\rho_k e_k+\tfrac12\rho_k u_k^2$ es la densidad de energía total de la fase. $k^*$ es la otra fase. $\lambda$ es la tasa de relajación de velocidad. Cada fase se cierra con su ecuación de estado. A menudo es un gas rígido.

$\gamma_k$ es el cociente de calores específicos. $p_{\infty,k}$ es una constante de presión para la repulsión entre moléculas. Un líquido que apenas se comprime, como el agua, tiene $p_\infty$ de varios GPa.

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

Las ecuaciones de momento y energía llevan $\partial\alpha_k/\partial x$ a la derecha. Ese término no se pliega en una divergencia. No se puede escribir como $\partial_t U+\partial_x F(U)=0$. El volumen finito de tipo Godunov tropieza delante de un término no conservativo así. Cuando $\alpha$ salta en una discontinuidad, no está claro qué valor multiplicar. Una discretización mala lanza oscilaciones falsas a través de la interfaz aunque la presión y la velocidad empiecen uniformes.

La prueba es el criterio de Abgrall. Si la presión y la velocidad empiezan uniformes, tienen que seguir uniformes sea cual sea el reparto de la fracción de volumen. Dos fases que corren a una presión y una velocidad solo deben transportar la interfaz, no fabricar una onda nueva. El registro de un salto de presión del 21 por ciento en un flujo que debía seguir uniforme es lo que ocurre cuando el término no conservativo se multiplica mal.

Los valores de interfaz $u_I$ y $p_I$ también son elecciones del modelo. La elección habitual de Saurel–Abgrall es una velocidad ponderada por masa y una presión ponderada por volumen.

Esa elección vuelve simétrico al modelo. Cada celda puede usar las mismas ecuaciones y la misma numérica.

Relajación: el camino al equilibrio#

Las fuentes $\mu(p_1-p_2)$ y $\lambda(u_{k^*}-u_k)$ son los operadores de relajación. Un $\mu$ y un $\lambda$ grandes arrastran las fases a una presión y una velocidad casi de inmediato. Mandar las tasas a infinito es el modelo de equilibrio más simple.

En la figura de arriba, una tasa de cero significa que las dos curvas no se encuentran. Ese es el no equilibrio que permiten las siete ecuaciones. Un $\mu$ mayor junta las curvas de presión en un punto con más pendiente. La presión y la velocidad se relajan por separado.

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

Suponer la relajación instantánea, forzar el equilibrio, y una ecuación desaparece dentro de una restricción algebraica. Atar las velocidades deja seis ecuaciones. Atar también las presiones deja el modelo de cinco ecuaciones de Kapila. Atar las temperaturas deja cuatro ecuaciones. BN se sienta en la cima de esa jerarquía.

El modelo de cinco ecuaciones es el más usado para capturar interfaces. Su velocidad del sonido de mezcla no es monótona en la fracción de volumen, la curva de Wood, así que la fuente de fracción de volumen es rígida. BN despega esa rigidez y la deja en una fuente de relajación. Numéricamente es el más flexible de los dos.

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

Se deja la relajación en el hueso. Las masas parciales y las energías internas de ambas fases se quedan fijas. Solo se mueve la fracción de volumen $\alpha$, por $d\alpha/dt=\mu(p_1-p_2)$. Se integra cómo una presión de gas rígido llega al equilibrio al cambiar $\alpha$.

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 corrida 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

Al pulsar un modelo cambia qué equilibrio quedó forzado y dónde la curva de velocidad del sonido deja de ser monótona. La fuente rígida de fracción de volumen aparece en el momento en que la jerarquía baja a cinco ecuaciones.

Los 324 bar impresos son dos términos enormes que se cancelan. Con $\alpha=0.5$ la densidad del líquido es 960. La energía interna específica es $3.93\times10^{8}/480=8.1875\times10^{5}$. Entonces $(\gamma-1)\rho e$ vale $2.6724\times10^{9}$ Pa y $\gamma p_\infty$ vale $2.64\times10^{9}$ Pa. La resta deja $3.24\times10^{7}$ Pa, 324 bar, el 1.2 por ciento restante.

En $\alpha=0.5$ la presión del líquido cambia cerca de $-5.34\times10^{9}$ Pa por unidad de $\alpha$. Un desplazamiento de 0.001 son 53 bar. El movimiento impreso va de 0.500 a 0.506. La presión del gas $4.4\times10^{5}/(1-\alpha)$ vale 8.8 bar en 0.5 y cerca de 8.9 bar junto a 0.506. El equilibrio es el líquido bajando hasta el gas.

El primer paso usa $\mu=4\times10^{-8}$, paso de tiempo $10^{-3}$ y brecha $3.152\times10^{7}$ Pa. El producto ronda $1.26\times10^{-3}$, un quinto del 0.006 impreso. Los 75 pasos son la cola que se encoge. La tolerancia de 1 Pa importa al final.

La fracción de volumen se movió solo un 0.6 por ciento, de 0.500 a 0.506, mientras la presión del líquido cayó de 324 bar a 8.9 bar. El agua modelada como gas rígido apenas se comprime, así que un cambio mínimo de volumen fabrica un cambio enorme de presión. El equilibrio se sienta cerca de la presión del gas por la misma razón. El líquido rígido puede igualar la presión sin expandirse mucho.

El primer paso de relajación que ya mueve 0.001 en fracción de volumen#

El primer paso ya mueve la fracción de volumen en $1.26\times10^{-3}$. La presión del líquido baja decenas de bar solo con ese movimiento. La presión interfacial del término no conservativo multiplica ese colapso por el gradiente de la interfaz. Poner la relajación y el flujo en la misma etapa fija el tamaño de una onda falsa en este primer paso.

¿A dónde va ese producto si la advección y la presión se separan del todo dentro del flujo? La separación TV que aparta la advección de la presión en el flujo de Baer–Nunziato deja abierto si esa separación cubre el término no conservativo. Cuando 324 bar se derrumban en un paso, ¿el flujo de presión ya separado deja quieta la interfaz?