Refinar la malla a la mitad adelantó el estallido al doble — autovalores complejos del modelo two-fluid de 6 ecuaciones

Refinar la malla a la mitad adelantó el estallido al doble#

Cuando un solver revienta, el primer movimiento habitual es refinar la malla. Si el error baja, el problema está en la discretización. Si nada cambia, el problema está en el modelo físico. Ese suele ser el orden del descarte.

Pero hay un caso en que la aguja se mueve al revés. Se reduce el tamaño de celda a la mitad y el estallido llega exactamente el doble de rápido. Se cuadruplica el número de celdas y llega cuatro veces antes. Reducir el paso temporal no altera la tasa de crecimiento.

Al subir sigma de 0 a 1.1 los puntos del plano complejo caen sobre el eje real. Al bajar slip u_r a 0 la perturbación deja de crecer.

El instante en que los puntos pasan a verde es el mismo en que el panel derecho deja de amplificar. El carril de malla fina llega antes al techo. Los números se cuadran más abajo.

Ese síntoma no es un error de discretización. Indica que las ecuaciones de gobierno son mal planteadas (ill-posed: el estado en que una perturbación crece a una tasa inversamente proporcional a su longitud de onda, sin cota). Cuando Pandare y Luo construyeron un solver two-fluid de volúmenes finitos basado en densidad en su artículo AIAA de 2018, este fue el primer punto que hubo que atender.

Los autovalores del modelo de seis ecuaciones con presión única se extrajeron directamente, para ver de dónde sale el par complejo y qué coeficiente necesita la presión interfacial antes de volver al eje real. La respuesta es exactamente 1.

Un par de ondas lentas abandona el eje real#

El modelo two-fluid trata ambas fases como continuos que se interpenetran. Masa, cantidad de movimiento y energía se resuelven por separado para cada fase. Al colapsar las dos presiones en una sola ($p_g = p_l \equiv p$) quedan seis EDP. Eso es el modelo de Wallis, también llamado modelo de seis ecuaciones con presión única.

Si se apaga por un momento la compresibilidad y se cuasi-linealiza en una dimensión con las variables primitivas $(\alpha_g, p, u_g, u_l)$, el par de ondas lentas tiene autovalor en forma cerrada.

Aquí $\alpha_k$ es la fracción volumétrica, $\rho_k$ la densidad de fase, $u_k$ la velocidad de fase y $\sigma$ el coeficiente del término de presión interfacial que se describe más abajo. El primer término es una velocidad media ponderada por densidad; el segundo es cuánto se separan las dos ondas.

Todo depende de lo que hay bajo el radical. Para $\sigma < 1$ se vuelve negativo y los dos autovalores forman un par complejo conjugado. Mientras el deslizamiento $|u_l - u_g|$ no sea nulo, esto ocurre sin falta. Es decir: en el instante en que las dos fases se mueven a velocidades distintas, el modelo queda mal planteado. Con deslizamiento nulo el radicando desaparece entero.

Una tabla — siete ecuaciones, seis, y seis con corrección#

Hay tres opciones que compiten por este lugar. Puestas en columnas se ve qué compra y qué vende cada una.

El modelo de siete ecuaciones da a la fracción volumétrica su propia ecuación de transporte. Eso compra hiperbolicidad, a costa de términos de relajación de presión que traen rigidez consigo. Esa estructura quedó ordenada en una entrada anterior sobre flux splitting para Baer–Nunziato. El inconveniente es que el modelo se justifica físicamente sobre todo para partículas densamente empaquetadas y suspensiones. Encaja mucho peor en una tubería con agua y aire estratificados.

Los autovalores 4×4, desde Python#

Antes de confiar en la forma cerrada, conviene resolver el sistema original tal cual. Se arma $A W_t + B W_x = 0$ con la compresibilidad retenida y se toman los autovalores de $A^{-1}B$. Aire y agua, $\alpha_g = 0.5$, con el gas adelantado a 10 m/s.

import numpy as np
 
def interfacial_dp(a, rg, rl, ur, sigma):
    """Correccion de Stuhmiller: p_int = p - dp"""
    return sigma * a * (1 - a) * rg * rl * ur**2 / (a * rl + (1 - a) * rg)
 
def two_fluid_matrices(a, rg, rl, cg, cl, ug, ul, sigma):
    """A W_t + B W_x = 0,  W = (alpha_g, p, u_g, u_l)"""
    dp = interfacial_dp(a, rg, rl, ul - ug, sigma)
    kg, kl = a / (rg * cg**2), (1 - a) / (rl * cl**2)
    A = np.array([[ 1.0, kg,  0.0,    0.0],
                  [-1.0, kl,  0.0,    0.0],
                  [ 0.0, 0.0, a * rg, 0.0],
                  [ 0.0, 0.0, 0.0,    (1 - a) * rl]])
    B = np.array([[ ug,  ug * kg, a,          0.0],
                  [-ul,  ul * kl, 0.0,        1 - a],
                  [ dp,  a,       a * rg * ug, 0.0],
                  [-dp,  1 - a,   0.0,        (1 - a) * rl * ul]])
    return A, B
 
def char_speeds(sigma, a=0.5, rg=1.2, rl=1000.0, cg=340.0, cl=1500.0, ug=10.0, ul=0.0):
    A, B = two_fluid_matrices(a, rg, rl, cg, cl, ug, ul, sigma)
    return np.linalg.eigvals(np.linalg.solve(A, B))
 
print("air/water, alpha_g=0.5, u_g=10, u_l=0 m/s")
print("sigma   max|Im lambda|   slow pair Re")
for s in [0.0, 0.5, 0.9, 1.0, 1.1, 1.5]:
    lam = char_speeds(s)
    slow = np.sort(lam.real)[1:3]
    print("%5.2f   %12.5f   %8.4f %8.4f" % (s, np.abs(lam.imag).max(), slow[0], slow[1]))

air/water, alpha_g=0.5, u_g=10, u_l=0 m/s
sigma   max|Im lambda|   slow pair Re
 0.00        0.34614     0.0120   0.0120
 0.50        0.24481     0.0120   0.0120
 0.90        0.10967     0.0120   0.0120
 1.00        0.00718     0.0120   0.0120
 1.10        0.00000    -0.0972   0.1212
 1.50        0.00000    -0.2326   0.2566

En $\sigma = 0$ la parte imaginaria vale 0.346 m/s. Las partes reales de las dos ondas lentas quedan juntas en 0.0120. El gas corre a 10 m/s y sin embargo la onda avanza a 0.012 m/s, por el peso de la densidad: el agua es 830 veces más pesada que el aire, así que la media se arrastra hacia el líquido.

Al subir $\sigma$ la parte imaginaria se encoge. En 1.1 llega a cero y el par se separa en $-0.097$ y $0.121$. El residuo de 0.00718 en $\sigma = 1.0$ proviene de la compresibilidad. La forma cerrada se dedujo en el límite incompresible, de modo que una velocidad del sonido finita empuja el umbral un pelo por encima de 1.

Al reescribir la media ponderada por densidad, el 0.0120 deja de ser opaco. Con $\alpha_g = 0.5$, $\rho_g = 1.2$, $\rho_l = 1000$ y $u_g = 10$, el denominador es $0.5\times 1.2 + 0.5\times 1000 = 500.6$. El numerador es $0.5\times 1.2\times 10 = 6$, así que la velocidad media queda en $6/500.6 = 0.0120$. El agua pesa $1000/1.2 \approx 833$ veces más que el aire, y esa frase es toda la división. El gas corre a 10 m/s. La parte real del par lento sigue pegada a 0.012 m/s.

Las partes imaginarias siguen la forma cerrada con la misma fidelidad. El cociente entre 0.34614 en $\sigma = 0$ y 0.24481 en $\sigma = 0.5$ es $\sqrt{2} \approx 1.414$. En $\sigma = 0.9$ el factor es $\sqrt{0.1}$ y la parte imaginaria vale 0.10967. $|\mathrm{Im}\,\lambda|$ escala como $\sqrt{1-\sigma}$. En $\sigma = 1.1$ las partes reales se abren a $-0.0972$ y $0.1212$, y el punto medio sigue en 0.012. La semianchura 0.109 coincide con la parte imaginaria de $\sigma = 0.9$. Solo cambió el signo dentro de $(\sigma-1)$.

La forma cerrada sitúa el umbral exactamente en 1#

Ahora toca medir los mismos números con la forma cerrada y bisecar el $\sigma$ crítico.

from math import sqrt, pi
 
def material_pair(a, sigma, rg=1.2, rl=1000.0, ug=10.0, ul=0.0):
    """Forma cerrada en el limite incompresible para las dos ondas lentas"""
    al = 1.0 - a
    den = al * rg + a * rl
    mean = (al * rg * ug + a * rl * ul) / den
    disc = (sigma - 1.0) * a * al * rg * rl * (ul - ug) ** 2 / den**2
    if disc >= 0.0:
        return (mean - sqrt(disc), mean + sqrt(disc)), 0.0
    return (mean, mean), sqrt(-disc)
 
print("closed form vs the 4x4 eigenvalues above")
for s in [0.0, 0.5, 0.9, 1.1, 1.5]:
    (r1, r2), im = material_pair(0.5, s)
    print("sigma=%4.2f  Re = %8.4f %8.4f   |Im| = %8.5f" % (s, r1, r2, im))
 
print()
print("growth rate of the shortest resolved mode, L = 1 m, sigma = 0")
_, im0 = material_pair(0.5, 0.0)
for n in [50, 100, 200, 400, 800]:
    k = pi * n          # k = pi / dx, dx = 1/n
    print("N=%4d  dx=%7.5f  k=%8.1f 1/m  growth=%8.2f 1/s" % (n, 1.0 / n, k, k * im0))
 
print()
print("critical sigma (incompressible limit) for a few states")
for a in [0.1, 0.5, 0.9]:
    for ur in [1.0, 30.0]:
        lo, hi = 0.0, 5.0
        for _ in range(60):
            mid = 0.5 * (lo + hi)
            _, im = material_pair(a, mid, ug=ur)
            if im > 0.0: lo = mid
            else: hi = mid
        print("alpha_g=%.1f  u_r=%4.1f  ->  sigma_c = %.6f" % (a, ur, hi))

closed form vs the 4x4 eigenvalues above
sigma=0.00  Re =   0.0120   0.0120   |Im| =  0.34599
sigma=0.50  Re =   0.0120   0.0120   |Im| =  0.24466
sigma=0.90  Re =   0.0120   0.0120   |Im| =  0.10941
sigma=1.10  Re =  -0.0974   0.1214   |Im| =  0.00000
sigma=1.50  Re =  -0.2327   0.2566   |Im| =  0.00000
 
growth rate of the shortest resolved mode, L = 1 m, sigma = 0
N=  50  dx=0.02000  k=   157.1 1/m  growth=   54.35 1/s
N= 100  dx=0.01000  k=   314.2 1/m  growth=  108.70 1/s
N= 200  dx=0.00500  k=   628.3 1/m  growth=  217.40 1/s
N= 400  dx=0.00250  k=  1256.6 1/m  growth=  434.79 1/s
N= 800  dx=0.00125  k=  2513.3 1/m  growth=  869.58 1/s
 
critical sigma (incompressible limit) for a few states
alpha_g=0.1  u_r= 1.0  ->  sigma_c = 1.000000
alpha_g=0.1  u_r=30.0  ->  sigma_c = 1.000000
alpha_g=0.5  u_r= 1.0  ->  sigma_c = 1.000000
alpha_g=0.5  u_r=30.0  ->  sigma_c = 1.000000
alpha_g=0.9  u_r= 1.0  ->  sigma_c = 1.000000
alpha_g=0.9  u_r=30.0  ->  sigma_c = 1.000000

Con sigma en 0 el carril de malla fina llega primero a la línea de blow-up. Al pasar de 1 los cuatro carriles se tumban juntos. Los 54.35 1/s de N = 50 y los 869.58 1/s de N = 800 se separan desde la misma salida.

La forma cerrada coincide con los autovalores 4×4 hasta la tercera cifra decimal. Barriendo la fracción volumétrica de 0.1 a 0.9 y el deslizamiento de 1 a 30 m/s, el umbral se queda en 1.000000. Por eso el coeficiente de la corrección de Stuhmiller,

no es una perilla de ajuste arbitraria. $\sigma = 1$ es el coeficiente más pequeño que lleva el radicando exactamente a cero. En producción se deja margen y se usa un valor algo mayor que 1.

Inestabilidad y mal planteamiento son cosas distintas#

Un esquema numéricamente inestable mejora al reducir el paso temporal. Un problema mal planteado no, porque la tasa de crecimiento escala con el número de onda.

Al reducir $\Delta x$ a la mitad, la longitud de onda más corta representable también se reduce a la mitad y la tasa de crecimiento se duplica. En la salida anterior, 54.35 1/s en $N = 50$ pasa a 869.58 1/s en $N = 800$: exactamente un factor 16. Cuanto más fina la malla, más rápido muere la respuesta. La reparación pertenece al cierre, no al ancho de malla.

La tabla multiplica del mismo modo. En $N = 50$, $k = 157.1$ 1/m y $|\mathrm{Im}\,\lambda| = 0.34599$, de modo que $157.1 \times 0.34599 \approx 54.35$ 1/s. Al duplicar $N$ también se duplica $k$: 54.35, 108.70, 217.40, 434.79, 869.58, cada uno el doble exacto del anterior. De $N = 50$ a $N = 800$ la malla es 16 veces más fina y el crecimiento es 16 veces mayor. Una perturbación que se amplifica 54 veces por segundo con $\Delta x = 0.02$ m queda cerca de 870 veces por segundo con $\Delta x = 0.00125$ m. Reducir el paso temporal deja intacto el número de onda 157.1. Lo que se encoge es $\Delta t$, no $k$.

En códigos reales el síntoma suele quedar oculto. La difusión numérica del esquema upwind de primer orden aporta un amortiguamiento de orden $O(k^2 \Delta x)$ que cancela el crecimiento lo bastante como para que la corrida siga adelante. Por eso un código que se porta bien a bajo orden revienta en cuanto se sube el orden. La difusión numérica estaba pagando en silencio una deuda del modelo.

El resto de la tabla — cómo sobrevive un solver basado en densidad a bajo Mach#

Recuperar la hiperbolicidad no cierra el asunto. Casi todas las aplicaciones reales de flujo multifásico viven a número de Mach muy bajo, donde un solver basado en densidad queda encadenado al CFL acústico y el paso temporal se derrumba.

Ese terreno ha pertenecido tradicionalmente a los métodos basados en presión. Suponer un campo de velocidad solenoidal borra la velocidad del sonido de las ecuaciones, así que el CFL depende solo de la velocidad del flujo. El precio es que la compresibilidad nunca se trata con rigor. Basta con introducir un fenómeno de alta temperatura como la ebullición para que los errores crezcan.

Pandare y Luo eligen otra ruta: mantener la formulación basada en densidad pero transformar a las variables primitivas $[p, v, T]$ y resolver de forma totalmente implícita. Poner la presión como incógnita acondiciona mucho mejor el sistema a bajo Mach. Los términos de fuerza interfacial como arrastre y masa virtual también se tratan implícitamente, relajando aún más el límite de paso temporal.

El lado del flujo numérico carga con un compromiso del mismo tipo. Cuando una onda de choque fuerte se encuentra con una interfaz material, AUSM$^+$-up produce presiones negativas. El remedio establecido era llamar a un solver de Riemann exacto solo en esas caras, pero las iteraciones de Newton salen caras. El artículo añade en su lugar un término de acoplamiento de fracción volumétrica al flujo de masa, que compra la misma robustez: en esencia, disipación tipo Lax–Friedrichs proporcional al salto de fracción volumétrica. La exigencia de no perturbar una interfaz en reposo es la misma restricción que devolver en cada paso una interfaz emborronada mediante antidifusión.

No subir el orden mientras σ esté por debajo de 1#

Si un solo 4×4 con deslizamiento residual ya muestra parte imaginaria, reducir la malla a la mitad es un diagnóstico, no una reparación. Una corrida cuyo crecimiento pasa de 54.35 1/s en $N = 50$ a 869.58 1/s en $N = 800$ está enferma en el modelo, no en el esquema. Lo primero es leer $\sigma$.

Por debajo de 1, la corrida vive de la difusión numérica $O(k^2 \Delta x)$ del upwind de primer orden. Al subir el orden ese colchón desaparece. Solo después de dejar $\sigma$ cerca de 1.1 el refinamiento de malla empieza a recortar el error en lugar de alimentarlo. Cuando un solver basado en densidad ya aguantó el régimen de bajo Mach, la decisión siguiente es si montar la presión de forma explícita, sin una Poisson.