Menos ecuaciones, y dejó de reventar — el modelo de cuatro ecuaciones y la condición de equilibrio de interfaz

La presión uniforme se mantiene solo si Γ misma se advecta; dividir un ρΓ conservado entre ρ es lo que la dobla.

Al borrar una ecuación, el código dejó de reventar. El flujo multifásico compresible suele enseñarse al revés: un modelo de siete ecuaciones, con presión, velocidad y temperatura propias de cada fase, es el seguro, y menos ecuaciones significa más riesgo. Collis (2025) invierte ese orden. El modelo más corto, cuatro ecuaciones, es el que aguanta. La bisagra es el equilibrio de interfaz.

Con el interruptor Naive encendido y γ del gas derecho en 1.66 o en 1.09, la línea roja de presión en la interfaz sube y la verde se queda plana.

right gas γ = 1.66CFL = 0.60

Uniform velocity, uniform pressure, only γ jumps across the shaded interface. The naive scheme grows a pressure spike at the interface; the consistent one keeps P flat to machine precision.

El pico rojo crece con el salto de γ. Que la línea verde siga en el piso en el paso 900 dice si la condición de equilibrio sigue viva.

La escalera de los equilibrios — 7, 6, 5, 4 ecuaciones#

Un método de interfaz difusa necesita una regla para las dos fases dentro de una celda emborronada. Esa regla es lo que cuenta como equilibrio.

Siete ecuaciones (Baer–Nunziato): cada fase tiene su velocidad, su presión y su temperatura. Lo más general y lo más caro.
Seis ecuaciones: una velocidad. Las presiones siguen separadas y se encuentran por relajación.
Cinco ecuaciones: velocidad y presión compartidas. La temperatura sigue aparte.
Cuatro ecuaciones: velocidad, presión y temperatura compartidas. El más corto.

Una interfaz física mide nanómetros. La numérica es esa lámina inflada a unas pocas celdas. Dentro de esa banda, el equilibrio térmico y mecánico es prácticamente inmediato. Cerca de la interfaz, cuatro ecuaciones son la física más cercana.

T^{(1)} = T^{(2)}, \quad P^{(1)} = P^{(2)}, \quad \mathbf{u}^{(1)} = \mathbf{u}^{(2)}

Los superíndices $(1),(2)$ son las dos fases. El modelo de cuatro ecuaciones impone las tres igualdades dentro de la celda.

Por qué menos ecuaciones son más robustas#

Cinco ecuaciones o más suelen añadir un transporte extra, de fracción de volumen o de algo como $1/(\gamma-1)$ . A nivel de EDP esas ecuaciones repiten información. Discretizadas, se separan. La misma cantidad, dos caminos, dos valores.

Cuatro ecuaciones quitan el extra. Solo se conservan masa, momento y energía. La presión cierra con una ecuación de estado de mezcla. Menos ecuaciones, menos sitios donde desviarse. También más barato.

Queda un problema. El margen que daban esas ecuaciones extra desaparece, y una discretización conservativa aparta la presión en la interfaz.

El momento en que se rompe la condición de equilibrio de interfaz#

La condición de equilibrio de interfaz (IEC) es simple. Si presión, temperatura y velocidad empiezan uniformes, el paso de la interfaz no debe romper esa uniformidad.

Un esquema que la rompe hace crecer una oscilación de presión que no estaba. La oscilación crece con el tiempo. Con un cociente de densidad grande aparecen densidad negativa o energía interna negativa, y la corrida se detiene.

Por qué lo hace un esquema conservativo. En una mezcla de gases ideales, con $\Gamma \equiv 1/(\gamma-1)$ , la energía total es

E = P\,\Gamma + \tfrac{1}{2}\,\rho\,u^2

$E$ es la energía total, $P$ la presión, $\rho$ la densidad y $u$ la velocidad. $\Gamma$ salta en la interfaz porque los fluidos difieren.

Se avanza una celda con upwind de primer orden bajo $u$ y $P$ uniformes. $\nu = u\,\Delta t/\Delta x$ es el número CFL. La actualización conservativa de la energía cancela el término $+P$ y deja

E_i^{n+1} = E_i - \nu\,(E_i - E_{i-1})

Al sustituir la definición de energía, el lado derecho se parte.

E_i^{n+1} = P_0\big[\,\Gamma_i - \nu(\Gamma_i - \Gamma_{i-1})\,\big] + \tfrac{1}{2}\,u_0^2\,\rho_i^{n+1}

Para que la presión nueva siga en $P_0$ , $\Gamma$ tiene que actualizarse justo como ese corchete. Eso es un upwind no conservativo. El camino conservativo advecta $\rho\Gamma$ entero y vuelve a dividir entre $\rho$ . Las dos respuestas se separan en la interfaz. Una presión uniforme que aun así se agita cuando el camino no conservativo no cuadra es el sitio que marca el criterio de Abgrall. Esa separación es la oscilación de presión.

La línea roja que creció al mover γ es esta separación. Al pasar a IEC-consistent, la verde sigue plana en el paso 900.

Python — la oscilación y su cura#

El mismo cálculo. $u$ y $P$ uniformes, solo salta $\gamma$ , upwind de primer orden. Solo cambia el esquema.

import numpy as np
 
def advect_interface(scheme, gamma_R=1.66, nu=0.6, steps=300, N=200):
    x = np.linspace(0.0, 1.0, N, endpoint=False)
    ramp = np.clip((x - 0.28) / 0.06, 0.0, 1.0)   # interfaz difuminada en unas celdas
    gamma = 1.4 + (gamma_R - 1.4) * ramp          # aire (1.4) -> gas elegido
    rho = 1.0 + (0.2 - 1.0) * ramp                # la densidad tambien salta
    u0, P0 = 1.0, 1.0
    Gam = 1.0 / (gamma - 1.0)                      # rigidez de la mezcla 1/(gamma-1)
    mom = rho * u0
    E = P0 * Gam + 0.5 * rho * u0**2               # E = P*Gamma + 0.5*rho*u^2
 
    for _ in range(steps):
        rhoG = rho * Gam                           # rho*Gamma antes de actualizar
        r_up, m_up, e_up = np.roll(rho, 1), np.roll(mom, 1), np.roll(E, 1)
        rho = rho - nu * (rho - r_up)              # masa: upwind conservativo
        mom = mom - nu * (mom - m_up)              # momento: conservativo
        E   = E   - nu * (E   - e_up)              # energia: conservativo
        if scheme == "naive":
            rhoG = rhoG - nu * (rhoG - np.roll(rhoG, 1))  # transporta rho*Gamma conservativo
            Gam = rhoG / rho                              # y vuelve a dividir por rho
        else:  # consistent
            Gam = Gam - nu * (Gam - np.roll(Gam, 1))      # transporta Gamma directo (no conservativo)
 
    P = (E - 0.5 * mom**2 / rho) / Gam             # recupera la presion
    return float(np.max(np.abs(P - P0)))
 
for s in ("naive", "consistent"):
    print(f"{s:11s}  max|P-P0| = {advect_interface(s):.3e}")

La salida se separa por trece cifras, no por dos.

naive        max|P-P0| = 3.83e-02
consistent   max|P-P0| = 4.44e-15

Naive aparta la presión cerca de un 4% en la interfaz. Consistent queda plano hasta el redondeo de punto flotante. Lo hizo una línea. Recuperar $\Gamma$ dividiendo entre rho, o advectar $\Gamma$ misma.

Puestos los números en $\Gamma = 1/(\gamma-1)$ , el aire con $\gamma=1.4$ tiene $\Gamma=2.5$ . El helio con $\gamma=1.66$ tiene $\Gamma=1/0.66 \approx 1.515$ . El SF6 con $\gamma=1.09$ tiene $\Gamma=1/0.09 \approx 11.11$ . De aire a helio el salto es cerca de 0.985. De aire a SF6, cerca de 8.61. El salto del SF6 es unas 8.74 veces el del helio. Ese factor es por qué la línea roja crece cuando γ va a 1.09.

Un solo paso del $\nu=0.6$ del listado, en una celda abrupta. La celda derecha es $\rho=0.2$ , $\Gamma\approx 1.515$ ; la vecina upwind es aire, $\rho=1$ , $\Gamma=2.5$ . Si se advecta $\rho\Gamma$ de forma conservativa y se divide entre $\rho$ , $\Gamma_{\mathrm{naive}}\approx 2.384$ . Si se advecta $\Gamma$ directo, $\Gamma_{\mathrm{cons}}\approx 2.106$ . Al invertir la relación de energía de arriba, $P/P_0 = \Gamma_{\mathrm{cons}}/\Gamma_{\mathrm{naive}} \approx 0.883$ , así que un paso da $|P-P_0|\approx 0.117$ . La misma celda abrupta con SF6 da $|P-P_0|\approx 0.692$ . El $3.83\times 10^{-2}$ impreso no es esa celda abrupta. El ancho de la rampa $0.06$ son 12 celdas con $\Delta x=1/200$ , y en 300 pasos la interfaz recorre $\nu \cdot 300 \cdot \Delta x = 0.9$ . Una interfaz emborronada, empujada tanto, deja el máximo por debajo del salto de un paso, en $0.0383$ . El $0.117$ del paso abrupto partido por ese $0.0383$ da cerca de $3.05$ . Repartido en 12 celdas, el error sigue en unos pocos por ciento. El $4.44\times 10^{-15}$ consistent es el redondeo que queda cuando los dos caminos coinciden.

Lo que el artículo realmente hizo#

El juguete de arriba solo empuja una $\gamma$ de gas ideal. Collis (2025) es más amplio.

Una ecuación de estado NASG (Noble–Abel stiffened gas, una EOS que llega a líquidos) cierra agua, aire, helio y SF6 juntos. La presión de mezcla es una cuadrática cerrada por la ley de Amagat.
Una reconstrucción ENO de alto orden (WENO/TENO) sustituye al upwind de primer orden. La reconstrucción va atada a la hipótesis de equilibrio de cuatro ecuaciones, así que la IEC se cumple a precisión de máquina.
Un limitador que preserva la positividad, extendido a cuatro ecuaciones, cubre un choque fuerte contra una interfaz de gran cociente de densidad. Densidad y cuadrado de la velocidad del sonido no pasan a negativo.
Nada de esto añade una ecuación de fracción de volumen o de calor específico, ni suelta la conservación ni un filtro escalar.

También hay límites. El líquido es un solo componente y el gas son componentes ideales. Más componentes exigen un solver iterativo para el equilibrio de presión y temperatura. El artículo sigue siendo un preprint.

Recuperar Γ dividiendo entre ρ dobla una presión plana en un paso#

En una celda abrupta, el helio ya tiene $|P-P_0|\approx 0.117$ . El SF6 tiene $0.692$ . La salida de 300 pasos, $0.0383$ , es ese salto repartido en una rampa de 12 celdas. Basta dejar la línea que advecta $\rho\Gamma$ y divide entre $\rho$ para que un $P_0$ uniforme se doble en la interfaz de inmediato. Cambiar las ecuaciones de más por esa línea no conservativa es toda la decisión de cuatro ecuaciones.

Cinco ecuaciones comparten velocidad y presión y aun así guardan una ecuación extra. Si esa pieza es la que derrumba la velocidad del sonido en la interfaz, queda abierto si borrarla para salvar la IEC borra también la enfermedad o solo la esconde en la EOS de mezcla. Si la velocidad del sonido de cinco ecuaciones que se derrumba al mezclar agua y aire a medias sigue dentro del cierre de cuatro ecuaciones es el corte siguiente.