Dejar que el sonido y la materia fluyan por separado — División acústico-convectiva (Lagrange–Projection)#

El sonido sigue a 340 m/s junto a quien camina a 5 km/h. Un solver compresible mete ambas velocidades en un solo paso. La razón es el problema. En un flujo lento, las ondas acústicas adelantan a la materia por cientos. Un Godunov directo queda atado a las ondas rápidas. La convección lenta se difumina.

Al bajar el control de Mach cambian a la vez la velocidad a la que se parte el pulso de presión y el sitio donde se queda el grumo de entropía.

En 0.02 el pulso ámbar se abre en dos ondas que escapan a los lados. El grumo teal casi no se mueve. Los tres triángulos quedan en $u\!-\!c$, $u$ y $u\!+\!c$. Con menos Mach solo se detiene el del centro. La separación exterior sigue en $2c$.

ten Eikelder et al. (2017) separan ambas. Parte acústica, parte convectiva, cada una con su esquema. Estas notas ponen esa división Lagrange–Projection sobre las ecuaciones de Euler monofásicas (single-phase) y la contrastan con un pulso acústico en un flujo lento.

“

Artículo: M.F.P. ten Eikelder, F. Daude, B. Koren, A.S. Tijsseling, "An acoustic-convective splitting-based approach for the Kapila two-phase flow model", Journal of Computational Physics 331 (2017) 188–208. DOI: 10.1016/j.jcp.2016.11.031

El jacobiano se parte en dos#

Al escribir las ecuaciones de Euler 1D en variables primitivas $\mathbf{W}=(\rho,u,p)^T$, toman forma cuasilineal.

Aquí $\rho$ es la densidad, $u$ la velocidad y $p$ la presión. La matriz de coeficientes $\mathbf{B}$ se parte en dos piezas.

$\mathbf{A}$ contiene todos los términos de presión: la parte acústica (acoustic). $\mathbf{C}=u\mathbf{I}$ simplemente arrastra la materia: la parte convectiva (convective). En realidad, la división equivale a arrancar el término $u\partial_x$ de la derivada lagrangiana $\mathrm{D}/\mathrm{D}t = \partial_t + u\partial_x$.

La descomposición aditiva de los autovalores#

Por qué esta división es tan potente se ve al mirar los autovalores. Las velocidades de onda del sistema completo son $u-c,\ u,\ u+c$. Se suman exactamente a partir de una contribución acústica y una convectiva.

$c=\sqrt{\gamma p/\rho}$ es la velocidad del sonido. Las velocidades acústicas son $\pm c$, independientes de la velocidad del flujo. La velocidad convectiva es $u$ en todos los casos. En el límite de bajo Mach $M=u/c\to 0$, la banda acústica conserva su ancho $2c$ mientras la velocidad convectiva se encoge a cero. Esa brecha es la rigidez (stiffness). La separación de los triángulos de la figura de arriba es esta misma suma. Cuando $M\to 0$ solo desaparece $\lambda^c$. $\lambda^a=\pm c$ se queda.

Paso acústico — HLLC en coordenadas lagrangianas#

La parte acústica es la expansión y compresión que impulsa la presión. El artículo la lleva a coordenadas de masa (lagrangianas) y la resuelve con un solver de Riemann tipo HLLC. El estado con asterisco en la cara $j+1/2$ se reduce a dos fórmulas.

$a_{j+1/2}=\max(\rho_j c_j,\ \rho_{j+1}c_{j+1})$ es la impedancia acústica (densidad por velocidad del sonido). Solo con esas velocidades y presiones con asterisco se realiza la actualización de las variables eulerianas. La tasa de expansión de cada celda se comprime en un único coeficiente.

$R_j$ indica cuánto estiró o comprimió el paso acústico el volumen de la celda. La densidad sale de inmediato como $\rho^{n+1-}_j = \rho^n_j / R_j$.

El fondo del cálculo es $c_0=1$, $\rho_0=1$, $\gamma=1.4$, de modo que $p_0=\rho_0 c_0^2/\gamma=1/1.4$. El pulso de presión tiene amplitud $10^{-3}$. La anomalía de densidad tiene amplitud $5\times 10^{-2}$. Los centros están en $x=0.5$ y $x=0.25$, así que al inicio no se solapan. El máximo de densidad arranca en el grumo de entropía, cerca de $1.05$. Si $R_j$ supera 1, la celda se expande y la densidad baja. Si queda por debajo de 1, comprime. Como $\Delta t/\Delta x$ multiplica el salto de $u^*$ en las caras, un CFL de 0.8 mantiene $R_j$ cerca de 1. Con 400 celdas en un dominio de longitud 1, $\Delta x=1/400$.

Paso convectivo — proyección upwind#

Ahora la velocidad con asterisco $u^*$ del paso acústico arrastra la materia. Se actualizan las cantidades conservadas $\varphi\in\{\rho,\rho u,\rho E\}$ con un simple esquema upwind.

$\varphi^{\text{up}}_{j+1/2}$ es $\varphi_j$ cuando $u^*\!\ge\!0$, y $\varphi_{j+1}$ en caso contrario. El $R_j$ del paso anterior regresa aquí, y la masa, el momento y la energía se conservan de forma exacta. Al recorrer los dos pasos en orden, un paso de tiempo queda completo.

Con $u^*$ positivo, el valor de cara sale de la celda izquierda. Con $u^*$ negativo, de la derecha. $R_j\,\varphi^{n+1-}$ deshace primero el cambio de volumen. El flujo de cara mueve la materia después. Si el orden se invierte, las sumas de masa, $\rho u$ y $\rho E$ se separan. Las caras periódicas no dejan escapar nada. Con la realimentación de $R_j$ en su sitio, la suma de momento conserva solo el redondeo. Cuánto se esparce un Roe directo a bajo Mach lo mide aparte las tres escalas de difusión artificial de Hope–Collins.

Python — pulso acústico en un flujo lento#

Aquí está toda la tubería en numpy: fronteras periódicas, un gas ideal, tres cantidades conservadas por celda. La condición inicial coloca un pequeño pulso de presión y una anomalía de densidad sobre un flujo de fondo $u_0=Mc_0$.

import numpy as np
 
GAMMA = 1.4  # razon de calores especificos del gas ideal
 
def primitives(rho, mom, Ene):
    """Conservadas -> variables primitivas (velocidad, presion, velocidad del sonido)."""
    u = mom / rho
    e = Ene / rho - 0.5 * u * u            # energia interna especifica
    p = (GAMMA - 1.0) * rho * e
    c = np.sqrt(GAMMA * p / rho)
    return u, p, c
 
def acoustic_faces(rho, u, p, c):
    """Estado acustico HLLC u*, p* en la cara j+1/2 (ec. 34 del articulo)."""
    rp, up, pp, cp = (np.roll(a, -1) for a in (rho, u, p, c))
    a = np.maximum(rho * c, rp * cp)       # impedancia acustica a = max(rho*c) (ec. 32)
    ustar = 0.5 * (u + up) + (p - pp) / (2 * a)
    pstar = 0.5 * (p + pp) + 0.5 * a * (u - up)
    return ustar, pstar
 
def lagrange_projection_step(rho, mom, Ene, dx, dt):
    """Paso acustico -> paso convectivo (ec. 38, 40 del articulo)."""
    u, p, c = primitives(rho, mom, Ene)
    uf, pf = acoustic_faces(rho, u, p, c)          # j+1/2
    uf_m, pf_m = np.roll(uf, 1), np.roll(pf, 1)    # j-1/2
    lam = dt / dx
    # 1) paso acustico: expansion/compresion impulsada por la presion
    R = 1.0 + lam * (uf - uf_m)                     # ec. (39)
    rho1 = rho / R
    mom1 = (mom - lam * (pf - pf_m)) / R
    Ene1 = (Ene - lam * (pf * uf - pf_m * uf_m)) / R
    # 2) paso convectivo: arrastrar la materia a u* (proyeccion upwind)
    def project(phi1):
        phi_f = np.where(uf >= 0, phi1, np.roll(phi1, -1))
        phi_f_m = np.roll(phi_f, 1)
        return R * phi1 - lam * (uf * phi_f - uf_m * phi_f_m)
    return project(rho1), project(mom1), project(Ene1)
 
def run_acoustic_pulse(mach, n=400, cfl=0.8, tmax=0.25):
    x = (np.arange(n) + 0.5) / n
    dx = 1.0 / n
    c0, rho0 = 1.0, 1.0
    p0 = rho0 * c0**2 / GAMMA
    u0 = mach * c0
    dp = 1e-3 * np.exp(-((x - 0.5) / 0.03)**2)     # pulso de presion acustico
    ds = 5e-2 * np.exp(-((x - 0.25) / 0.03)**2)    # anomalia de entropia (densidad) -> onda de contacto
    rho = rho0 + dp / c0**2 + ds
    u = np.full(n, u0)
    p = p0 + dp
    mom = rho * u
    Ene = p / (GAMMA - 1) + 0.5 * rho * u * u
    t, m0 = 0.0, mom.sum()
    while t < tmax:
        _, _, c = primitives(rho, mom, Ene)
        dt = min(cfl * dx / np.max(np.abs(mom / rho) + c), tmax - t)
        rho, mom, Ene = lagrange_projection_step(rho, mom, Ene, dx, dt)
        t += dt
    return rho, mom, m0, mom.sum()
 
for M in (0.02, 0.2):
    rho, mom, m0, m1 = run_acoustic_pulse(M)
    print(f"M={M}: error de momento={abs(m1-m0)/abs(m0):.1e}, rho_max={rho.max():.4f}")
# M=0.02: error de momento=2.2e-16, rho_max=1.0497
# M=0.2:  error de momento=0.0e+00, rho_max=1.0452

El momento se conserva a precisión de máquina. La presión se mantiene acotada y la densidad permanece positiva. La división no rompe la conservación porque el paso convectivo reintroduce el $R_j$.

El $\rho_{\max}$ impreso es 1.0497 con $M=0.02$ y 1.0452 con $M=0.2$. La cresta inicial, cerca de 1.05, apenas se erosiona. El pulso de presión solo suma $dp/c_0^2=10^{-3}$ a la densidad, así que el máximo lo pone el contacto. Errores de momento de $2.2\times 10^{-16}$ y 0 indican que la división en sí no es lo que rompe la conservación. El paso sigue limitado por $|u|+c$, de modo que subir $M$ de 0.02 a 0.2 casi no cambia el número de pasos. Con $c_0=1$, $|u|+c$ vale 1.02 frente a 1.2. La convección puede acelerarse por diez y la acústica sigue ocupando la mayor parte del paso.

Lo que de verdad se separa a bajo Mach#

La división rinde a bajo Mach. El paso estable de un método directo sigue encadenado a $|u|+c$, aunque la física de interés viva sobre $u$.

Al mover el control de Mach, la banda ámbar queda fija en un ancho $2c$ y solo se encoge la velocidad teal $u$.

El ancho acústico no depende de $M$. Solo $u$ converge a cero. Con $c=1$ y $M=0.05$, $u=0.05$ y $u+c=1.05$, así que $(u+c)/u=21$. El esquema dividido puede avanzar la convección con un paso distinto del acústico. La física lenta deja de ir a rastras de las ondas rápidas. La misma división evita la difusión numérica de más que un Godunov directo acumula a bajo Mach.

Una mirada crítica#

Hay tres tropiezos. La división usada aquí es de primer orden en el tiempo. El segundo orden pide un envoltorio de Strang —acústico, convectivo, acústico— y otra pasada en el costo. El escenario real del artículo es el modelo bifásico de cinco ecuaciones de Kapila. El término no conservativo $K\partial_x u$ de la fracción volumétrica y la garantía de positividad son el obstáculo de verdad. Este recorte monofásico los salta. La impedancia $a=\max(\rho c)$ es robusta y también disipativa. En un choque fuerte esa disipación puede difuminar el contacto. Por eso el artículo compara aparte precisión y costo frente al enfoque directo.

En un código, la bifurcación no es exótica. El rhoPimpleFoam basado en presión de OpenFOAM y la familia all-Mach ya tratan acústica y convección de forma implícita y explícita desde la misma raíz. Lagrange–Projection es la versión escrita en el lenguaje de los solvers de Riemann.

Sin R_j la conservación muere primero#

Si el $R_j$ acústico no vuelve a la proyección convectiva, las actualizaciones de masa, momento y energía dejan de coincidir. La ejecución de arriba llega a $2.2\times 10^{-16}$ solo con esa realimentación puesta. Un envoltorio de Strang puede subir el orden. La secuencia crece a acústico–convectivo–acústico y aparece una pasada más.

Al pasar al modelo de cinco ecuaciones de Kapila, $K\partial_x u$ en la fracción volumétrica sigue ahí. El recorte monofásico no toca ese término no conservativo. Con $a=\max(\rho c)$ el choque aguanta y el contacto se engrosa. Otra forma de sacar la acústica del integrador temporal es el IMEX acústico-transporte de Peluchon. En $M=0.05$, $(u+c)/u=21$. Si ese factor 21 se deja en el solver acústico o se queda en el CFL convectivo es la bifurcación de la implementación.