Implementación directa del modelo de 5 ecuaciones#

El agua detrás del choque está a presión $10^{9}$. El agua por delante está a $10^{5}$. En medio hay una burbuja de aire, y la fracción de volumen $\alpha_1$ salta entre casi 0 y casi 1. El modelo de interfaz difusa de 5 ecuaciones en 1D transporta ese salto con cinco ecuaciones. Los dos fluidos comparten una velocidad y una presión. Solo la fracción de volumen viaja aparte.

Al bajar a 0 el control de difusión artificial, el frente de fracción de volumen se corta en una sola celda.

Al subir el control, el frente se reparte en varias celdas. El $\alpha$ de esas celdas entra directo en el denominador de la presión de mezcla.

Resumen de las ecuaciones de gobierno#

En 1D las variables conservadas y el flujo son estos.

La quinta componente del flujo es 0. $\alpha_1$ es no conservativa, así que se actualiza aparte.

Cinco ecuaciones es lo que queda al reducir ambas fases a una velocidad y una presión. Devolver un sistema a cada fase es el modelo de siete ecuaciones de Baer–Nunziato.

Ecuación de estado: Stiffened Gas EOS#

Los dos fluidos usan un gas endurecido (stiffened gas).

La presión de mezcla sale de esta regla.

$p_k$ es la presión del fluido $k$, $e_k$ su energía interna, y $\gamma_k$, $p_{\infty,k}$ las constantes de la EOS. El denominador es $\sum \alpha_k/(\gamma_k-1)$.

Implementación en Python#

Configuración inicial y EOS#

import numpy as np
import matplotlib.pyplot as plt
 
# Dominio
N = 1000
x = np.linspace(0, 1, N)
dx = x[1] - x[0]
 
# Parámetros de la EOS
# Fluido 1: Aire (gas ideal)
gamma1, pinf1 = 1.4, 0.0
# Fluido 2: Agua (gas endurecido)
gamma2, pinf2 = 4.4, 6e8
 
def pressure(alpha1, arho1, arho2, E, u):
    """Calcula la presión de la mezcla a partir de las variables conservadas."""
    alpha2 = 1.0 - alpha1
    rho = arho1 + arho2
    ke = 0.5 * rho * u**2
    rhoe = E - ke
 
    num = rhoe - alpha1*gamma1*pinf1/(gamma1-1) - alpha2*gamma2*pinf2/(gamma2-1)
    den = alpha1/(gamma1-1) + alpha2/(gamma2-1)
    return num / den
 
def sound_speed(alpha1, arho1, arho2, p):
    """Calcula la velocidad del sonido de la mezcla (fórmula de Wood)."""
    alpha2 = 1.0 - alpha1
    rho = arho1 + arho2
    rho1 = arho1 / np.maximum(alpha1, 1e-12)
    rho2 = arho2 / np.maximum(alpha2, 1e-12)
 
    c1sq = gamma1 * (p + pinf1) / rho1
    c2sq = gamma2 * (p + pinf2) / rho2
 
    # Fórmula de Wood
    inv_rhoc2 = alpha1/(rho1*c1sq) + alpha2/(rho2*c2sq)
    return np.sqrt(1.0 / (rho * inv_rhoc2))

pressure es la regla de mezcla reescrita sobre variables conservadas. En sound_speed, $c_k^2=\gamma_k(p+p_{\infty,k})/\rho_k$ es la velocidad del sonido del fluido puro.

Cálculo del flujo HLLC#

HLLC cubre las cuatro leyes de conservación. La velocidad de contacto $S^\star$ vuelve junto al flujo y luego transporta $\alpha$.

def hllc_flux(UL, UR, alpha1L, alpha1R):
    """Flujo numérico HLLC para las 4 ecuaciones de conservación."""
    arho1L, arho2L, momL, EL = UL
    arho1R, arho2R, momR, ER = UR
 
    rhoL = arho1L + arho2L
    rhoR = arho1R + arho2R
    uL = momL / rhoL
    uR = momR / rhoR
    pL = pressure(alpha1L, arho1L, arho2L, EL, uL)
    pR = pressure(alpha1R, arho1R, arho2R, ER, uR)
    cL = sound_speed(alpha1L, arho1L, arho2L, pL)
    cR = sound_speed(alpha1R, arho1R, arho2R, pR)
 
    # Estimaciones de la velocidad de onda
    SL = min(uL - cL, uR - cR)
    SR = max(uL + cL, uR + cR)
    Sstar = (pR - pL + rhoL*uL*(SL-uL) - rhoR*uR*(SR-uR)) \
          / (rhoL*(SL-uL) - rhoR*(SR-uR))
 
    FL = np.array([arho1L*uL, arho2L*uL, momL*uL + pL, (EL+pL)*uL])
    FR = np.array([arho1R*uR, arho2R*uR, momR*uR + pR, (ER+pR)*uR])
 
    if SL >= 0:
        return FL, Sstar
    elif Sstar >= 0:
        coeff = rhoL * (SL - uL) / (SL - Sstar)
        UstarL = coeff * np.array([
            arho1L/rhoL,
            arho2L/rhoL,
            Sstar,
            EL/rhoL + (Sstar - uL)*(Sstar + pL/(rhoL*(SL-uL)))
        ])
        return FL + SL*(UstarL - UL), Sstar
    elif SR > 0:
        coeff = rhoR * (SR - uR) / (SR - Sstar)
        UstarR = coeff * np.array([
            arho1R/rhoR,
            arho2R/rhoR,
            Sstar,
            ER/rhoR + (Sstar - uR)*(Sstar + pR/(rhoR*(SR-uR)))
        ])
        return FR + SR*(UstarR - UR), Sstar
    else:
        return FR, Sstar

Integración temporal (1er orden)#

El avance temporal es Euler hacia adelante. El paso es $\mathrm{CFL}\,\Delta x/\max(|u|+c)$, con CFL 0.5.

def solve(U, alpha1, dx, t_end, cfl=0.5):
    """Integración temporal de Euler hacia adelante."""
    t = 0.0
    while t < t_end:
        rho = U[0] + U[1]
        u = U[2] / rho
        p = pressure(alpha1, U[0], U[1], U[3], u)
        c = sound_speed(alpha1, U[0], U[1], p)
 
        dt = cfl * dx / np.max(np.abs(u) + c)
        if t + dt > t_end:
            dt = t_end - t
 
        # Calcular flujos en las interfaces de las celdas
        flux = np.zeros((4, N+1))
        uface = np.zeros(N+1)
 
        for i in range(N+1):
            iL = max(i-1, 0)
            iR = min(i, N-1)
            UL_local = U[:4, iL]
            UR_local = U[:4, iR]
            flux[:, i], uface[i] = hllc_flux(
                UL_local, UR_local, alpha1[iL], alpha1[iR]
            )
 
        # Actualizar variables conservadas
        for k in range(4):
            U[k, 1:-1] -= dt/dx * (flux[k, 2:-1] - flux[k, 1:-2])
 
        # Actualizar alpha1 (no conservativo, upwind)
        for i in range(1, N-1):
            if uface[i] >= 0:
                alpha1[i] -= dt/dx * uface[i] * (alpha1[i] - alpha1[i-1])
            else:
                alpha1[i] -= dt/dx * uface[i] * (alpha1[i+1] - alpha1[i])
 
        alpha1 = np.clip(alpha1, 1e-8, 1-1e-8)
        t += dt
 
    return U, alpha1

La actualización de $\alpha_1$ no es un flujo conservativo. El signo de la velocidad de cara uface elige el lado upwind.

Condiciones iniciales: Onda de choque - Burbuja de gas#

Antes de la posición 0.25 el agua ya cruzó el choque: $\rho=1100$, $u=50$, $p=10^{9}$. Entre 0.4 y 0.6 está la burbuja de aire, $\rho=1$, $u=0$, $p=10^{5}$. El resto es agua ambiente, $\rho=1000$, $u=0$, $p=10^{5}$.

# Condición inicial: choque golpeando una burbuja de aire en agua
U = np.zeros((4, N))
alpha1 = np.zeros(N)
 
for i in range(N):
    if x[i] < 0.25:
        # Agua post-choque (alta presión)
        rho2 = 1100.0; u_val = 50.0; p_val = 1e9
        alpha1[i] = 1e-8
        arho1 = alpha1[i] * 1.0  # aire insignificante
        arho2 = (1 - alpha1[i]) * rho2
    elif 0.4 < x[i] < 0.6:
        # Burbuja de aire
        rho1 = 1.0; u_val = 0.0; p_val = 1e5
        alpha1[i] = 1 - 1e-8
        arho1 = alpha1[i] * rho1
        arho2 = (1 - alpha1[i]) * 1000.0
    else:
        # Agua ambiente
        rho2 = 1000.0; u_val = 0.0; p_val = 1e5
        alpha1[i] = 1e-8
        arho1 = alpha1[i] * 1.0
        arho2 = (1 - alpha1[i]) * rho2
 
    rho = arho1 + arho2
    E = compute_total_energy(alpha1[i], arho1, arho2, p_val, u_val)
 
    U[0, i] = arho1
    U[1, i] = arho2
    U[2, i] = rho * u_val
    U[3, i] = E

compute_total_energy no está en este fragmento. La EOS ya dice cómo volver de la presión a la energía: $\rho e=(p+\gamma p_\infty)/(\gamma-1)$.

Consideraciones importantes#

1. Limitación del rango de $\alpha$#

Si $\alpha_1$ vale exactamente 0 o 1, el producto conservado $\alpha_1\rho_1$ no se puede dividir entre $\alpha_1$. El código recorta $\alpha_1$ entre $10^{-8}$ y $1-10^{-8}$.

2. Velocidad del sonido de la mezcla: Fórmula de Wood#

Una celda mezclada puede ser más lenta que cualquiera de los dos fluidos puros. La fórmula lo dice, y el CFL lee esa misma velocidad.

Los estados puros de la condición inicial ya fijan el número. El aire de la burbuja tiene $\rho=1$, $p=10^{5}$, $\gamma=1.4$, $p_\infty=0$. Entonces $c=\sqrt{1.4\times 10^{5}}\approx 374\,\mathrm{m/s}$. El agua ambiente tiene $\rho=1000$, $p=10^{5}$, $\gamma=4.4$, $p_\infty=6\times 10^{8}$, así que $c\approx 1625\,\mathrm{m/s}$. El agua detrás del choque tiene $\rho=1100$, $p=10^{9}$, así que $c\approx 2530\,\mathrm{m/s}$. El grupo $p+p_\infty$ es cerca de 2.67 veces el del agua ambiente, mientras la densidad solo es $1100/1000=1.10$ veces mayor. La energía cinética a $u=50$ es $0.5\times 1100\times 50^{2}=1.375\times 10^{6}$. Eso es cerca del 0.13% de la densidad de energía interna $\rho e\approx 1.07\times 10^{9}$ en la misma celda. La energía de este choque está en $p_\infty$, no en los 50 m/s.

Al mezclar aire de la burbuja y agua ambiente con $\alpha_1=0.5$ y una presión común $10^{5}$, en la fórmula de Wood $1/(\rho c^2)$ es casi por completo el término del aire. $\rho_1 c_1^2=1.4\times 10^{5}$, mientras $\rho_2 c_2^2\approx 2.64\times 10^{9}$. La densidad de mezcla es $0.5\times 1+0.5\times 1000=500.5$, y $c_\mathrm{mix}\approx 24\,\mathrm{m/s}$. Ni 374 ni 1625.

El denominador de la presión es sensible en el mismo sitio. El agua pura tiene $\sum\alpha/(\gamma-1)=1/3.4\approx 0.294$. Si se filtra $\alpha_1=0.01$ en esa agua, el denominador pasa a $0.01/0.4+0.99/3.4\approx 0.316$, cerca de un 7.5% más. La presión se divide por ese denominador, así que un error de $\alpha$ en una interfaz emborronada aparece como error de presión.

El paso CFL es $\Delta t=0.5\,\Delta x/\max(|u|+c)$, y $\Delta x$ es el espaciado de linspace(0, 1, 1000), o sea $1/999$. No lo fija la celda mezclada. Detrás del choque, $|u|+c\approx 50+2530=2580$, de modo que $\Delta t\approx 1.9\times 10^{-7}$. Una celda a 24 m/s permitiría un paso cerca de 110 veces mayor. Una velocidad de mezcla baja no aprieta el paso de tiempo. Lo que hace es dejar la acústica de la interfaz casi quieta dentro de ese mismo $\Delta t$.

3. Extensión a alta precisión#

El integrador de arriba es de primer orden, y la actualización espacial usa el valor de la celda tal cual. El segundo orden pide MUSCL más un slope limiter. Por encima, WENO. La precisión temporal sube con Runge–Kutta. Nada de eso quita la actualización no conservativa de $\alpha$ ni la velocidad de Wood.

Interpretación de resultados#

Cuando el choque entra en la burbuja, un choque transmitido sale por el lado de aguas abajo y una rarefacción reflejada vuelve a la izquierda. La burbuja se aplasta y sube su presión interna. Un chorro en la pared de aguas abajo pide 2D o 3D. Las explosiones submarinas y la litotricia por ondas de choque se apoyan en esta figura.

Un α solo upwind deja que la presión de contacto se mueva primero#

El quinto flujo es 0, y $\alpha$ se pone upwind según el signo de la velocidad de cara. El salto del 7.5% en el denominador de la presión cuando $\alpha_1=0.01$ se filtra en el agua significa que una celda mala de interfaz puede mover la presión de equilibrio. Esta actualización de primer orden no impone ese equilibrio por su cuenta. La discretización del término no conservativo es lo que decide si la presión se queda quieta.

La comprobación de ese equilibrio es un contacto donde no debería pasar nada más. Dónde se escapa la presión en ese contacto es la pregunta del criterio de Abgrall y el término no conservativo.