[Reseña de artículo] La presión pasó a implícita y Δt se redujo 100 veces — el segundo 1/M² oculto en la velocidad relativa

La presión pasó a implícita y el paso temporal solo se hizo más pequeño#

Hay un fracaso por el que pasa cualquiera que escriba un resolvedor bifásico de bajo Mach. Se detecta que los términos de presión son rígidos (stiff) y se tratan de forma implícita. Creyendo que ya no hay dependencia de la velocidad del sonido, se sube $\Delta t$ a la escala de la velocidad material. Y estalla. En $M=0.01$ el paso temporal utilizable es cien veces menor que el que funcionaba en $M=1$.

El artículo de 2022 de Lukáčová-Medviďová y colaboradores señala con precisión dónde vive ese fracaso. El modelo bifásico isentrópico lleva $1/M^2$ en dos lugares. Uno es la presión, ya conocida. El otro va montado sobre la diferencia de entalpía en la ecuación de velocidad relativa. Este texto desentierra el segundo y mide, sobre una malla de 64 celdas, si tratarlo de forma implícita libera de verdad la cota del paso temporal del número de Mach.

Al empujar log10 M a la izquierda, los abanicos rápidos se abren juntos. La barra de $\Delta t$ admisible suelta el número de Mach solo con f^c implicit marcado.

El abanico azul es la presión; el rosa, la diferencia de entalpía. Con solo f^a implicit, la barra inferior sigue encogiéndose con $M$.

Dos copias de 1/M² en un mismo modelo#

El vector de estado del artículo es $W = (\rho,\ \alpha\rho,\ \alpha_1\rho_1,\ \rho u,\ u_1-u_2)^T$: densidad de mezcla, fracción volumétrica (la parte del volumen de celda que ocupa la fase 1), densidad parcial, cantidad de movimiento de mezcla y la diferencia de las dos velocidades de fase. Esa última componente distingue al modelo. En lugar de llevar una velocidad por fase, lleva una velocidad de mezcla y una diferencia de velocidad.

Al adimensionalizar, $M_k$ entra en el flujo. En la ecuación de cantidad de movimiento la presión aparece como

donde $\hat{p}_k$ es la presión truncada alrededor de un estado de referencia, $M_k$ el número de Mach de la fase $k$ y $\varrho_k$ una razón de densidades de referencia. Este término diverge cuando $M\to 0$: es el primer término rígido y se trata de forma implícita igual que en los esquemas IMEX monofásicos.

El problema está en la quinta ecuación. La diferencia de velocidad la impulsa una diferencia de entalpía:

Aquí $\hat{h}_k$ es la entalpía truncada de la fase $k$. En el denominador está el mismo $M_k^2$. Pasar la presión a implícita deja este término intacto.

Tres subsistemas y sus características en un solo gráfico#

El splitting del artículo corta el flujo en $f^a$ (acústica rápida), $f^b$ (transporte material) y $f^c$ (términos rígidos de mezcla), y añade al final el término fuente de relajación $r(W)$. Cada pieza tiene sus propios autovalores. El dial de arriba abrió juntos el abanico azul y el rosa por este splitting. Llevar también $f^c$ a implícito es todo lo que hace la sección 3.2.3.

Escribir los autovalores del subsistema de velocidad relativa#

En $f^c$ quedan tres ecuaciones: densidad parcial, cantidad de movimiento y diferencia de velocidad. Los autovalores de este sistema $3\times 3$ son

con el discriminante

donde $\chi$ es la fracción másica y el superíndice $RS$ marca valores de la solución de referencia. El segundo término de $\mathcal{A}$ conserva intacto $1/M_k^2$, de modo que $\lambda^{\pm} \sim \mathcal{O}(1/M)$.

Como $\mathcal{A} > 0$, los autovalores son reales: el subsistema es hiperbólico y por eso se puede resolver por separado. Pero su velocidad característica va como $1/M$. Tratado de forma explícita, la condición CFL queda $\Delta t \le \Delta x \cdot M / \text{const}$. Todo lo ganado al hacer implícita la presión se devuelve aquí.

La ecuación de cantidad de movimiento corresponde a $\lambda^0=0$, así que no restringe el paso temporal. Por eso el artículo acopla solo dos ecuaciones — densidad parcial y diferencia de velocidad — en la resolución implícita. En un IMEX que separa acústica y transporte bastaba con pasar a implícita la acústica rápida. Con dos fases hay dos términos rígidos. El $1/M^2$ de la diferencia de entalpía es el segundo.

La restricción verdadera al elegir el splitting#

El artículo enuncia dos condiciones sobre el splitting: cada subsistema debe ser hiperbólico por sí mismo y debe ser fácil de resolver. La primera es la exigente. Si los términos se agrupan a la ligera, los autovalores se vuelven complejos y el paso del subsistema se desestabiliza por su cuenta.

Cómo se resuelve el paso implícito también forma parte de la condición. Al sustituir la actualización de la diferencia de velocidad en la de la densidad parcial aparece un problema elíptico no lineal para $(\alpha\rho_1)^{\star\star\star}$, que el artículo resuelve con una iteración linealizada. La matriz de coeficientes es simétrica definida positiva porque las dos cantidades

son ambas positivas. Para $s$ esto se deduce de $\varrho_1(\alpha\rho_1)/\rho = \chi \in (0,1)$. Ese único signo sostiene la frase "se resuelve de forma eficiente con resolvedores lineales estándar". Un paso temporal que sobrevive cuando la velocidad del sonido sale de la ecuación no lo hereda sola la discretización. El AP-IMEX para gas no ideal anota esa grieta desde el lado de la velocidad del sonido.

Explícito e implícito sobre la misma malla, en Python#

Al linealizar en torno a $u_1-u_2 = 0$, el núcleo de $f^c$ se reduce a un sistema de ondas $2\times 2$,

con $m = \alpha\rho_1$ y $w = u_1-u_2$. La velocidad de onda es $\sqrt{s\tilde{s}}$ y, como $\tilde{s} \propto 1/M^2$, escala como $1/M$. Se corrieron ambos tratamientos sobre la misma malla periódica de 64 celdas. El caso explícito usa un flujo de Rusanov cuya viscosidad numérica emplea el $\lambda_c$ correcto; el implícito elimina $w$ y resuelve el problema elíptico resultante con un resolvedor tridiagonal periódico. $\Delta t$ se mantiene fijo en un CFL material de 0,5.

import math
 
N, L = 64, 1.0
dx = L / N
RHO, ALPHA, CHI, U_MAT = 1.0, 0.5, 0.5, 0.2
M0 = CHI * RHO                 # alpha*rho_1, estado de fondo
S = M0 * (1.0 - M0 / RHO)      # s(W) de la Ec. (3.25)
 
def stiff_coeff(mach):
    """s-tilde: coeficiente de entalpia, O(1/M^2)"""
    return 1.0 / (mach**2 * ALPHA) + 1.0 / (mach**2 * (1.0 - ALPHA))
 
def relative_wave_speed(mach):
    return math.sqrt(S * stiff_coeff(mach))
 
def initial_state():
    m = [M0 + 1e-3 * math.sin(2*math.pi*(i+0.5)*dx) for i in range(N)]
    return m, [0.0] * N
 
def explicit_enthalpy(mach, dt, nstep):
    """Rusanov en (3.19a)+(3.19c), termino de entalpia explicito"""
    k, lam = stiff_coeff(mach), relative_wave_speed(mach)
    m, w = initial_state()
    for _ in range(nstep):
        fm = [0.0]*N; fw = [0.0]*N
        for i in range(N):
            j = (i+1) % N
            fm[i] = 0.5*S*(w[i]+w[j]) - 0.5*lam*(m[j]-m[i])
            fw[i] = 0.5*k*(m[i]+m[j]) - 0.5*lam*(w[j]-w[i])
        m = [m[i] - dt/dx*(fm[i]-fm[i-1]) for i in range(N)]
        w = [w[i] - dt/dx*(fw[i]-fw[i-1]) for i in range(N)]
        if not all(abs(v) < 1e12 for v in m):
            return float('inf')
    return max(abs(v-M0) for v in m)
 
def periodic_tridiag(a, b, c, d):
    """Resuelve el sistema tridiagonal periodico por Sherman-Morrison"""
    n = len(d)
    gam = -b[0]
    bb = b[:]; bb[0] -= gam; bb[n-1] -= c[n-1]*a[0]/gam
    def thomas(rhs):
        cp = [0.0]*n; dp = [0.0]*n
        cp[0] = c[0]/bb[0]; dp[0] = rhs[0]/bb[0]
        for i in range(1, n):
            den = bb[i] - a[i]*cp[i-1]
            cp[i] = c[i]/den
            dp[i] = (rhs[i] - a[i]*dp[i-1])/den
        x = [0.0]*n; x[n-1] = dp[n-1]
        for i in range(n-2, -1, -1):
            x[i] = dp[i] - cp[i]*x[i+1]
        return x
    u = [0.0]*n; u[0] = gam; u[n-1] = c[n-1]
    y = thomas(d); z = thomas(u)
    fac = (y[0] + a[0]*y[n-1]/gam) / (1.0 + z[0] + a[0]*z[n-1]/gam)
    return [y[i] - fac*z[i] for i in range(n)]
 
def implicit_enthalpy(mach, dt, nstep):
    """Ec. (3.23): elimina w y resuelve el problema eliptico para m"""
    k = stiff_coeff(mach)
    beta = dt*dt*S*k/(dx*dx)
    a = [-beta]*N; b = [1.0 + 2*beta]*N; c = [-beta]*N
    m, w = initial_state()
    for _ in range(nstep):
        rhs = [m[i] - dt/(2*dx)*S*(w[(i+1) % N]-w[i-1]) for i in range(N)]
        m = periodic_tridiag(a, b, c, rhs)
        w = [w[i] - dt/(2*dx)*k*(m[(i+1) % N]-m[i-1]) for i in range(N)]
    return max(abs(v-M0) for v in m)
 
def max_stable_dt(mach, nstep=200):
    lo, hi = 1e-9, 1.0
    for _ in range(40):
        mid = math.sqrt(lo*hi)
        if explicit_enthalpy(mach, mid, nstep) == float('inf'):
            hi = mid
        else:
            lo = mid
    return lo
 
dt_mat = 0.5*dx/U_MAT
print(f"dx={dx:.5f}  material CFL dt={dt_mat:.5f}  (u={U_MAT})")
print(f"{'M':>8} {'lambda_c':>10} {'dt_exp_max':>11} {'ratio':>7} {'|dm|_exp':>10} {'|dm|_imp':>10}")
for mach in (1.0, 0.3, 0.1, 0.03, 0.01):
    lam = relative_wave_speed(mach)
    dte = max_stable_dt(mach)
    e = explicit_enthalpy(mach, dt_mat, 200)
    i = implicit_enthalpy(mach, dt_mat, 200)
    es = "blow-up" if e == float('inf') else f"{e:.3e}"
    print(f"{mach:8.3g} {lam:10.2f} {dte:11.2e} {dte/dt_mat:7.3f} {es:>10} {i:10.3e}")

dx=0.01562  material CFL dt=0.03906  (u=0.2)
       M   lambda_c  dt_exp_max   ratio   |dm|_exp   |dm|_imp
       1       1.00    1.87e-02   0.479    blow-up  1.718e-06
     0.3       3.33    5.61e-03   0.144    blow-up  1.771e-14
     0.1      10.00    1.87e-03   0.048    blow-up  4.731e-13
    0.03      33.33    5.60e-04   0.014    blow-up  6.318e-12
    0.01     100.00    1.87e-04   0.005    blow-up  5.195e-11

Al desmarcar enthalpy implicit y empujar log10 M por debajo de $-1$, con el $\Delta t$ material fijo la amplitud logarítmica toca el techo.

Al marcarlo de nuevo, la amplitud se apaga con ese mismo $\Delta t$. Al bajar material CFL hasta que el caso explícito revive, ese valor es la columna ratio de arriba.

Por qué la cota sigue a M exactamente#

En la columna dt_exp_max: cuando $M$ baja cien veces, de 1 a 0,01, la cota baja de $1.87\times 10^{-2}$ a $1.87\times 10^{-4}$, exactamente el mismo factor. No es casualidad. En $\lambda_c = \sqrt{s\tilde{s}}$, $s$ no depende de $M$ mientras $\tilde{s} \propto M^{-2}$, de modo que $\lambda_c \propto M^{-1}$. Al sustituirlo en la condición de estabilidad de Rusanov, $\lambda_c \Delta t/\Delta x \le 1$, resulta $\Delta t_{\max} \propto M$.

La columna ratio es la que duele en la práctica. En $M=0.01$ solo puede usarse el 0,5 % del $\Delta t$ que pedía el CFL material, lo que significa doscientas veces más pasos para el mismo tiempo físico. La columna implícita, en cambio, completó los 200 pasos con el $\Delta t$ fijo sin divergir; en $M=0.01$ eso corresponde a un CFL acústico de 256.

Conviene leer $\lambda_c$ al revés. Vale $1.00$ en $M=1$, $10.00$ en $M=0.1$ y $100.00$ en $M=0.01$. El $s=0.25$ y el $\tilde{s}=4/M^2$ del código dan $\sqrt{s\tilde{s}}=1/M$, los mismos casilleros impresos. Al meter el paso material $\Delta t=0.03906$ en $M=1$ sale $\lambda_c\Delta t/\Delta x=0.03906/0.01562=2.50$. Eso pasa el límite de Rusanov, 1, y por eso ratio imprime $0.479$. La cuenta $1.87\times10^{-2}/0.03906=0.479$ cierra. Una entalpía explícita no aguanta el CFL material de 0,5 ni en $M=1$. En $M=0.01$ la razón es $0.005$. Entonces $0.479/0.005\approx 96$, pegado al factor 100 del arranque.

|dm|_imp es otra columna. Frente a una amplitud inicial de $10^{-3}$, el residuo en $M=1$ es $1.718\times10^{-6}$ y en $M=0.01$ es $5.195\times10^{-11}$. Las cinco filas explícitas dicen blow-up. Solo la columna implícita termina 200 pasos con el $\Delta t$ fijo. El CFL acústico impreso, 256, muestra adónde se fue ese límite explícito. Si el $1/M^2$ de $\tilde{s}$ no pasa a implícito, ese número vuelve como CFL explícito.

Incondicionalmente estable no es incondicionalmente preciso#

Conviene leer con honestidad la columna |dm|_imp. La amplitud inicial era $10^{-3}$ y el caso implícito dejó $10^{-11}$ o menos. El Euler retrógrado prácticamente aniquiló el modo rápido. No es un error: es lo que un esquema AP (asymptotic preserving, que preserva el límite asintótico) está diseñado para hacer. Cuando $\Delta t$ supera con mucho el tiempo acústico, las ondas de esa escala no se resuelven, se filtran. Lo que interesa conservar es la dinámica en la escala de la velocidad material.

Aun así, ese amortiguamiento no debe confundirse con precisión. Para usar el mismo esquema en regímenes compresibles ($M \sim 1$) hay que recuperar los términos truncados $\bar{p}_k$ y $\bar{h}_k$ y sumarlos de forma explícita, que es la tarea de la sección 3.2.4. Esos términos de orden superior se anulan en el límite de bajo Mach con datos iniciales bien preparados, así que añadirlos explícitamente no reintroduce dependencia de $M$ en la condición CFL. La divergencia y el amortiguamiento ya se separan en la curva logarítmica justo después de la tabla impresa.

Dos términos de relajación, dos tratos distintos#

Terminado el sistema homogéneo, queda el término fuente de relajación: relajación de presión sobre la fracción volumétrica y fricción sobre la diferencia de velocidad.

Aquí $\tau$ es el tiempo de relajación de presión y $\zeta$ el coeficiente de fricción. Ambos pueden ser rígidos y ambos se integran de forma implícita, pero cuestan cosas muy distintas.

El término de fricción es lineal en $u_1-u_2$, así que el Euler retrógrado se cierra con una sola división:

La relajación de presión es otra cosa: las presiones de fase son no lineales en $\alpha$. El artículo aplica el método de Newton a $g(\alpha)=0$ y el jacobiano arrastra las velocidades del sonido:

Esta derivada diverge cuando $\alpha \to 0$ o $1$, justo donde Newton empieza a sufrir cerca de una fase pura. La costumbre de recortar $\alpha$ lejos de sus extremos queda explicada por esta sola expresión. Cabe notar que $\tau \to 0$ da relajación instantánea y $\tau \to \infty$ recupera el modelo homogéneo.

Dejar explícita la diferencia de entalpía vuelve a atar Δt a M#

Si solo el $1/M^2$ de la presión pasa a implícito, la ecuación de velocidad relativa conserva explícito $\hat{h}_1/M_1^2-\hat{h}_2/M_2^2$. Su velocidad característica va como $1/M$, así que en $M=0.01$ solo es usable el $0.5\%$ del CFL material. Ahí están el ratio impreso $0.005$ y el CFL acústico 256.

La relajación de presión es otro paso implícito. $\partial g/\partial\alpha$ diverge cuando $\alpha\to 0$ o $1$, y Newton se atasca cerca de una fase pura. Resolver con dos presiones y luego fundirlas en una ya tiene libro propio en la relajación de presión de Saurel. Sacar la diferencia de entalpía de esa lista hace que la acústica implícita no devuelva el paso temporal.