[Revisión de artículo] Una interfaz agua-aire en flujo uniforme generó un salto de presión de 8×10⁸ Pa — la condición de no perturbación en un modelo BN basado en presión

Una interfaz en un flujo uniforme y la presión saltó 8×10⁸ Pa#

Considérese un tubo 1D con velocidad uniforme de 1 m/s y presión uniforme de 1 bar. Un bloque de agua viaja dentro, rodeado de aire. Físicamente no debería ocurrir nada. Sin embargo, al subir solo el flujo de masa a segundo orden y dejar la fracción volumétrica en primer orden, la presión en las celdas de interfaz se desvió $8\times10^{8}$ Pa tras 120 pasos. Este artículo lee el trabajo de Re y Abgrall de 2022 (IJNMF, A pressure-based method for weakly compressible two-phase flows under a Baer–Nunziato type model) en cinco preguntas y respuestas. Al final queda claro cómo discretizar los términos no conservativos de un código multifásico para que ese número vuelva a cero.

P1. ¿Qué problema aborda el artículo?#

El objetivo es el flujo bifásico débilmente compresible (weakly compressible: número de Mach bajo y cambios de densidad pequeños). Ejemplos: el chapoteo en un tanque o bolsas de aire atrapadas junto a una pared. El líquido es casi incompresible, pero una mezcla con burbujas es muy compresible. Por eso hay que conservar la compresibilidad de ambas fases.

El modelo es el sistema tipo Baer–Nunziato de Saurel–Abgrall (BN: siete ecuaciones donde cada fase tiene su propia presión, velocidad y temperatura). El artículo mantiene finitos los parámetros de relajación de presión y velocidad (términos fuente que empujan ambas fases hacia el equilibrio). La ecuación de estado (EOS) puede ser general: stiffened gas, Peng–Robinson y otras.

Los solvers basados en densidad necesitan precondicionamiento a bajo Mach. Pero una mezcla bifásica no tiene una única velocidad del sonido bien definida. Por eso los autores eligieron una formulación basada en presión. El escalado de Bijl–Wesseling divide la presión en dos variables y los términos acústicos se integran de forma implícita. Así, los cálculos fueron estables con CFL acústico mayor que 10. Una malla escalonada (staggered grid) evita las oscilaciones de tablero de ajedrez. La motivación es la misma que en la división de Toro–Vázquez que traslada la rigidez a un subsistema de presión.

P2. ¿Qué es la condición de no perturbación?#

La condición de no perturbación (non-disturbance condition) la formuló Abgrall en 1996. Un estado con presión y velocidad uniformes debe seguir uniforme cuando una interfaz lo atraviesa. En las ecuaciones continuas es trivial. En las discretas no lo es.

La ecuación de cantidad de movimiento de la fase $\sigma$ en el modelo BN contiene el término no conservativo $P_I\,\partial\alpha_\sigma/\partial x$.

$\alpha_\sigma$ es la fracción volumétrica y $P_I$ la presión de interfaz. Con $P$ uniforme, el término de presión de la izquierda, $P\,\partial\alpha/\partial x$, cancela exactamente el lado derecho. Para que la cancelación sobreviva a la discretización, ambos términos deben usar la misma diferencia. La ecuación (49) del artículo da la respuesta.

$k$ es la celda escalonada entre las celdas $i$ e $i+1$. El gradiente de presión se discretiza como $-[(\alpha P)_{i+1}-(\alpha P)_i]$. Es la misma diferencia sobre las mismas dos celdas, así que la suma es cero con $P$ uniforme.

En la simulación de abajo se puede cambiar directamente la discretización de $H_P$.

Con averaged (gradiente central en el centro de celda promediado a la cara) aparecen fuerzas rojas en las caras junto a la interfaz y la velocidad del agua se aleja de 1. Ensanchar la interfaz de 1 a 6 celdas reduce la fuerza, pero no la elimina. La ecuación (49) mantiene las barras verdes —cero— para cualquier ancho y presión.

P3. ¿Por qué la ecuación de fracción volumétrica debe seguir al flujo de masa?#

La misma lógica vale para la ecuación de fracción volumétrica $\partial_t\alpha_\sigma + u_I\,\partial_x\alpha_\sigma = \dots$. La densidad de fase se recupera como $\rho_\sigma = (\alpha\rho)_\sigma/\alpha_\sigma$. Si numerador y denominador reciben difusión distinta, una $\rho_\sigma$ uniforme cambia en las celdas de interfaz.

Con una EOS linealizada alrededor de la velocidad del sonido $c_\sigma$, el error de presión de mezcla cabe en una línea.

$\rho_{0,\sigma}$ es la densidad de referencia. El corchete es la diferencia entre "lo que transportó la ecuación de masa" y "lo que transportó la ecuación de fracción volumétrica". El agua tiene $c^2 = 2.25\times10^6$ m²/s², así que un desajuste de solo 1 kg/m³ se convierte en 2.25 MPa.

Por eso la ecuación (45) escribe $u_I\,\partial_x\alpha$ con la misma forma que el flujo de Rusanov de la ecuación de masa (un flujo de primer orden cuya difusión la fija la velocidad de onda máxima local).

Como el método es basado en presión, la velocidad de onda $S$ no incluye la del sonido. El flujo de Rusanov de la ecuación de masa también usa $S=\max|u|$. Ambas ecuaciones comparten el mismo $S$, de modo que una $\rho_\sigma$ uniforme se conserva. Es decir, $\hat F$ no es un flujo de cara verdadero. Es un operador que imita al flujo de masa.

P4. ¿Qué muestra el código?#

Una interfaz agua (1000 kg/m³, 1500 m/s) – aire (1.2 kg/m³, 340 m/s) se desplaza 0.3 con $u=1$. La fracción volumétrica y las masas de ambas fases se resuelven en cuatro combinaciones. Después, con $P$ uniforme, se mide el residuo de cara $R_k$ de la cantidad de movimiento para dos elecciones de $H_P$.

import numpy as np
 
# interfaz agua-aire 1D en u, p uniformes (debe seguir uniforme)
N, L, u, CFL = 200, 1.0, 1.0, 0.5
dx = L / N
dt = CFL * dx / u
x = (np.arange(N) + 0.5) * dx
rho0 = np.array([1000.0, 1.2])     # densidades de referencia agua, aire [kg/m^3]
c0 = np.array([1500.0, 340.0])     # velocidad del sonido por fase [m/s]
p0, eps = 1.0e5, 1e-6
 
def minmod_slope(q):
    a = q - np.roll(q, 1)
    b = np.roll(q, -1) - q
    return np.where(a * b > 0, np.sign(a) * np.minimum(abs(a), abs(b)), 0.0)
 
def face_flux(q, order):
    # flujo de Rusanov con S=|u| (basado en presión: sin sonido), u>0
    qL = q + 0.5 * minmod_slope(q) if order == 2 else q
    return u * qL
 
def central_flux(q):
    return 0.5 * u * (q + np.roll(q, -1))
 
def mixture_dp(a1, m1, m2):
    a = np.array([a1, 1 - a1])
    rho = np.array([m1, m2]) / a
    p = p0 + (c0**2)[:, None] * (rho - rho0[:, None])
    return np.max(abs((a * p).sum(0) - p0))
 
def march(alpha_rule, mass_order, t_end=0.3):
    a1 = np.where(abs(x - 0.3) < 0.15, 1 - eps, eps)
    m1, m2 = rho0[0] * a1, rho0[1] * (1 - a1)
    n = int(round(t_end / dt))
    for _ in range(n):
        Fa = central_flux(a1) if alpha_rule == "central" else face_flux(a1, alpha_rule)
        F1, F2 = face_flux(m1, mass_order), face_flux(m2, mass_order)
        a1 = np.clip(a1 - dt / dx * (Fa - np.roll(Fa, 1)), eps, 1 - eps)
        m1 = m1 - dt / dx * (F1 - np.roll(F1, 1))
        m2 = m2 - dt / dx * (F2 - np.roll(F2, 1))
    return n, mixture_dp(a1, m1, m2)
 
cases = [("A  alpha 1st / mass 1st", 1, 1),
         ("B  alpha 1st / mass 2nd", 1, 2),
         ("C  alpha 2nd / mass 2nd", 2, 2),
         ("D  alpha central / mass 1st", "central", 1)]
print("uniform u=1, p=1e5 Pa, N=200, CFL=0.5")
for name, ar, mo in cases:
    n, dp = march(ar, mo)
    print(f"{name:30s} steps={n:4d}  max|p_mix - p0| = {dp:10.3e} Pa")
 
def hp_residual(a, P, rule):
    # fuente de cantidad de movimiento en la cara escalonada k (entre i e i+1)
    grad = -(np.roll(a, -1) * P - a * P)          # -[(aP)_{i+1} - (aP)_i]
    if rule == "jump":
        HP = P * (np.roll(a, -1) - a)             # ecuación (49) del artículo
    else:
        g = 0.5 * (np.roll(a, -1) - np.roll(a, 1))   # gradiente central en el centro de celda
        HP = P * 0.5 * (g + np.roll(g, -1))          # promediado a la cara
    return grad + HP
 
print()
print("momentum source on faces, uniform P=1e5 Pa, water-rich faces (alpha>0.5)")
for w in (1, 3, 6):
    aw = 0.5 * (1 + np.tanh((0.15 - abs(x - 0.3)) / (w * dx)))
    aw = np.clip(aw, eps, 1 - eps)
    af = 0.5 * (aw + np.roll(aw, -1))
    for rule in ("jump", "averaged"):
        R = hp_residual(aw, p0, rule)
        du = dt * R / (dx * af * rho0[0])
        wet = af > 0.5
        print(f"width={w}  H_P={rule:8s}  max|R| = {np.max(abs(R)):9.3e} Pa"
              f"  du_water/step = {np.max(abs(du[wet])):9.3e} m/s")

uniform u=1, p=1e5 Pa, N=200, CFL=0.5
A  alpha 1st / mass 1st        steps= 120  max|p_mix - p0| =  5.116e-07 Pa
B  alpha 1st / mass 2nd        steps= 120  max|p_mix - p0| =  8.024e+08 Pa
C  alpha 2nd / mass 2nd        steps= 120  max|p_mix - p0| =  1.769e-05 Pa
D  alpha central / mass 1st    steps= 120  max|p_mix - p0| =  1.665e+09 Pa
 
momentum source on faces, uniform P=1e5 Pa, water-rich faces (alpha>0.5)
width=1  H_P=jump      max|R| = 1.319e-11 Pa  du_water/step = 6.777e-15 m/s
width=1  H_P=averaged  max|R| = 1.203e+04 Pa  du_water/step = 1.874e+00 m/s
width=3  H_P=jump      max|R| = 8.185e-12 Pa  du_water/step = 4.670e-15 m/s
width=3  H_P=averaged  max|R| = 8.326e+02 Pa  du_water/step = 4.142e-01 m/s
width=6  H_P=jump      max|R| = 1.283e-11 Pa  du_water/step = 6.439e-15 m/s
width=6  H_P=averaged  max|R| = 1.126e+02 Pa  du_water/step = 8.670e-02 m/s

El código anterior puede ejecutarse directamente en el navegador.

Con los valores por defecto, tras 120 pasos los cuatro números coinciden con el bloque de salida. Al bajar c_water de 1500 a 300, los errores de B y D caen unas 25 veces. Es justo el escalado con $c^2$ de la P3. A y C se quedan en nivel de redondeo con cualquier ajuste.

A y C quedan en $5\times10^{-7}$ y $2\times10^{-5}$ Pa. Lo que separó los resultados no fue el orden, sino si las dos ecuaciones usan la misma diferencia. B subió solo la masa a MUSCL. En las celdas de interfaz, $(\alpha\rho)_1$ se difumina menos y $\alpha_1$ más. En el lado del aire queda menos masa de agua que volumen de agua, y $\rho_1$ se aleja mucho de su referencia. La diferencia central de D hace oscilar $\alpha$, y esa oscilación se convierte en presión.

En la cantidad de movimiento, la ecuación (49) da $10^{-11}$ Pa: precisión de máquina. El gradiente promediado produce un residuo de cara de $1.2\times10^4$ Pa con una interfaz de una celda. La velocidad del agua cambia 1.87 m/s en un solo paso. El flujo mismo avanza a 1 m/s, así que el primer paso lo invierte.

P5. ¿Entonces hay que renunciar al segundo orden?#

No. Como muestra C, basta con usar la misma reconstrucción en ambas ecuaciones. Minmod es homogénea de grado uno. Así, si $(\alpha\rho)_\sigma = \rho_0\alpha_\sigma$, la pendiente también se escala por $\rho_0$ y el cociente se mantiene. Un limitador no lineal aplicado por separado a cada masa de fase puede romper esta propiedad. En ese caso conviene reconstruir variables primitivas como $\alpha$ y $\rho_\sigma$ y formar el flujo de masa con su producto.

Hay dos salvedades. Primera: en este experimento $u$ es constante, por lo que $u_I\,\partial_x\alpha$ coincide con $\partial_x(\alpha u)$. En un código real $u_I$ varía de celda a celda. Por eso el artículo usa el valor de celda $(u_I)_i$ y no uno de cara: así se garantiza $H_u=0$ con $\alpha$ uniforme. Segunda: el propio artículo califica esta discretización de "robusta, pero dependiente del esquema". Si se cambia el flujo de masa, hay que volver a derivar $H_u$ y $H_P$. Los intentos de cambiar solo la reconstrucción de la fracción volumétrica, como un limitador compresivo, son arriesgados justo aquí.

Un resumen breve de los demás resultados del artículo. Los tubos de choque agua–aluminio y agua–aire y una mezcla vapor–líquido de CO₂ coincidieron con las soluciones de referencia dentro de unos pocos por ciento. Las diferencias cerca de los choques son las esperables en un modelo no conservativo. Con parámetros de relajación finitos, el acercamiento de ambas fases al equilibrio de presión y velocidad aparece directamente en la solución.

Una prueba de interfaz uniforme antes de cambiar el próximo limitador#

Se coloca una interfaz nítida en $u$ y $p$ uniformes y se avanzan 100 pasos. Se imprime la desviación máxima de la presión de mezcla. Cualquier valor por encima de la precisión de máquina indica que la ecuación de $\alpha$ y la de masa están discretizadas de forma inconsistente. Para la cantidad de movimiento, basta revisar el residuo de cara $R_k$ tras el primer paso. La prueba debe repetirse cada vez que cambie un flujo o un limitador. Son pocas líneas, pero atrapan un error de 8×10⁸ Pa el primer día.