[Reseña] Cómo no matar el sonido a bajo Mach — esquema IMEX-AP de ondas

Meter un flujo de bajo Mach (low Mach) en un código compresible hace pensar que una velocidad del sonido mayor deja la figura más nítida. Ocurre lo contrario. Al bajar el número de Mach, un esquema upwind estándar borra de la malla una onda que estaba sana. Arun, Das Gupta y Samantaray (2019) aíslan ese colapso con el sistema lineal de ondas y arman un esquema asintóticamente preservante (AP) con IMEX-RK, estable sin depender del Mach. La razón de integrar el operador acústico de forma implícita está en los números de abajo.

Al mover $\varepsilon$ hacia la izquierda, la curva roja upwind sube con pendiente $-1$ y la curva cian del AP permanece plana.

En $\varepsilon=10^{-3}$ el cociente entre las dos curvas llega a cientos. Lo que crece al bajar el Mach no es la onda. Es la viscosidad numérica upwind.

Cuando ε tiende a cero, la malla se come la onda#

Las ecuaciones de Euler isentrópicas escaladas llevan un parámetro pequeño, el número de Mach $\varepsilon$. El término de presión entra en la cantidad de movimiento como $\nabla p / \varepsilon^2$.

$\rho$ es la densidad, $\mathbf{u}$ es la velocidad y $p=\rho^\gamma$ es la presión. $\varepsilon$ es la velocidad de referencia dividida por la velocidad del sonido de referencia.

Cuando $\varepsilon \to 0$ la solución va al límite incompresible. Es una perturbación singular (singular perturbation). Los autores acortan el análisis con el sistema lineal de ondas.

$\varrho$ es la perturbación linealizada de densidad y $a$ es la velocidad del sonido de ese estado. La velocidad efectiva del sonido es $c=a/\varepsilon$. Al reducir el Mach, esa velocidad acústica estalla.

El espacio well-prepared — las soluciones que deben sobrevivir#

En el límite $\varepsilon\to 0$ la solución no se escapa a cualquier lado. Se junta en el espacio well-prepared: densidad constante, velocidad de divergencia nula.

Ese espacio es el núcleo del operador acústico $L(U)=(a\nabla\cdot\mathbf{u},\,a\nabla\varrho)$. La solución discreta tiene que caer ahí también. Dellacherie (2010) mostró que un esquema que rompe esa invariancia pierde precisión a bajo Mach. Conservar el núcleo es ser asintóticamente preciso (AA). Dejar el límite de estabilidad independiente de $\varepsilon$ es ser asintóticamente preservante (AP).

La viscosidad upwind crece en proporción a la velocidad del sonido#

El culpable es la viscosidad numérica upwind. Diagonalizar el sistema de ondas en variables características lo parte en dos advecciones escalares a velocidades $\pm c$. Un paso upwind de Rusanov en cada modo tiene esta amplificación.

$\nu$ es el CFL acústico, $\theta$ es la fase por celda y $k$ es el número de onda.

Un paso duele poco. Lo que duele es la pila. Llegar a un tiempo físico fijo $T$ pide $M = T/\Delta t = T c /(\nu\Delta x)$ pasos acústicos. Con $\theta$ pequeño el decaimiento por paso es $\tfrac{1}{2}\nu(1-\nu)\theta^2$, y $M$ de ellos dan este decaimiento logarítmico.

Bajar el Mach un orden de magnitud sube la viscosidad efectiva otro orden. La pendiente $-1$ de la curva roja de arriba es esa proporción.

División de operadores — convección explícita, acústica implícita#

Al separar las ecuaciones por escala de tiempo aparecen dos operadores.

$H$ es convección, escala $O(1)$. $L/\varepsilon$ es acústica, escala $O(\varepsilon)$. Toda la rigidez vive en $L/\varepsilon$. Solo esa parte se integra de forma implícita. La convección queda explícita. La mezcla es IMEX (Implicit-Explicit). Si esa misma separación aguanta el paso de tiempo hasta Mach 0 es la pregunta de una reproducción IMEX TVD.

Con la acústica implícita, las diferencias centradas se quedan estables. La viscosidad upwind sale. $\Delta t$ lo fija la convección lenta, no el sonido rápido. Dentro de este modelo, un $\varepsilon$ menor no es motivo para encoger el paso.

IMEX-RK absorbe la rigidez#

La integración temporal es IMEX Runge-Kutta (IMEX-RK). El tableau explícito $(\tilde A,\tilde b)$ lleva $H$. El tableau implícito $(A,b)$ lleva $L/\varepsilon$. En cada etapa la acústica se condensa en una ecuación elíptica para la densidad, parecida a una Poisson de presión.

$\beta$ es el coeficiente diagonal del tableau IMEX y $\nabla^2$ es el laplaciano. La invertibilidad de ese problema elíptico sale de la teoría de puntos de silla (saddle point). En la malla, las matrices circulantes dan la misma conclusión. Solución única, estabilidad uniforme independiente de $\varepsilon$, e invariancia del espacio well-prepared.

Python — reproducir el decaimiento de amplitud de un modo#

El código de abajo calcula la amplificación upwind de un modo de Fourier e imprime la amplitud que queda tras el tiempo físico $T$, un Mach cada vez. El esquema implícito conserva el modo, así que su amplitud se queda en 1.

import numpy as np
 
a  = 1.0        # velocidad del sonido linealizada
k  = 2 * np.pi  # único modo de Fourier en [0,1]
T  = 1.0        # tiempo físico final fijo
nu = 0.5        # número de CFL acústico
 
def upwind_amp_factor(eps, N):
    """|g| de un modo característico bajo Rusanov en el sistema de ondas."""
    c  = a / eps              # velocidad acústica c = a/ε
    dx = 1.0 / N
    th = k * dx               # fase por celda θ = kΔx
    re = 1.0 - nu * (1.0 - np.cos(th))
    im = nu * np.sin(th)
    return np.hypot(re, im)   # |g| <= 1
 
def retained_amplitude(eps, N):
    """Amplitud del modo tras tiempo físico T (upwind vs implícito)."""
    c  = a / eps
    dx = 1.0 / N
    dt = nu * dx / c          # paso de CFL acústico
    M  = T / dt               # número de pasos explícitos
    explicit = upwind_amp_factor(eps, N) ** M
    implicit = 1.0            # central-implícito preserva el modo
    return M, explicit, implicit
 
if __name__ == "__main__":
    N = 64
    print(f"{'eps':>8} {'steps M':>10} {'explicit':>12} {'AP':>6}")
    for eps in [0.4, 0.2, 0.1, 0.05, 0.02, 0.01]:
        M, ex, ap = retained_amplitude(eps, N)
        print(f"{eps:8.3f} {M:10.0f} {ex:12.3e} {ap:6.2f}")

La salida es esta.

     eps    steps M     explicit     AP
   0.400        320    6.804e-01   1.00
   0.200        640    4.629e-01   1.00
   0.100       1280    2.143e-01   1.00
   0.050       2560    4.591e-02   1.00
   0.020       6400    4.430e-04   1.00
   0.010      12800    1.980e-07   1.00

Al bajar el control Mach $\varepsilon$, la amplitud roja llega al suelo antes de que termine el tiempo físico $T=1$. Al empujar el CFL $\nu$ hacia 1, $\nu(1-\nu)$ se reduce y el decaimiento afloja un momento. Bajar $\varepsilon$ todavía más entierra esa ganancia bajo la acumulación $1/\varepsilon$. Subir $N$ ayuda solo en parte: $\theta^2\propto N^{-2}$ y $\Delta x\propto N^{-1}$ no se cancelan.

Meter los números del código en la fórmula del conteo de pasos la cierra. Con $a=1$, $T=1$, $\nu=0.5$ y $N=64$, $\Delta x=1/64$ y $c=1/\varepsilon$. Entonces $\Delta t=\nu\Delta x/c$, de modo que $M=T/\Delta t=128/\varepsilon$. Los pasos impresos 320, 640, 1280, 2560, 6400 y 12800 son esa fórmula en $\varepsilon=0.4,0.2,0.1,0.05,0.02,0.01$. Partir $\varepsilon$ por dos duplica $M$. Partir $\varepsilon$ por diez multiplica $M$ por diez.

El decaimiento logarítmico $\tfrac12 T c(1-\nu)k^2\Delta x$ con los mismos números es el exponente $0.154/\varepsilon$. El coeficiente sale de $k=2\pi$. En $\varepsilon=0.4$ el exponente es 0.386 y $\mathrm{e}^{-0.386}=0.680$, igual al $6.804\times 10^{-1}$ impreso hasta el tercer dígito. En $\varepsilon=0.1$, $\mathrm{e}^{-1.542}=0.214$ frente al $2.143\times 10^{-1}$ impreso. En $\varepsilon=0.01$ el exponente es 15.4. $\mathrm{e}^{-15.4}$ queda cerca de $2.0\times 10^{-7}$, y la tabla imprime $1.980\times 10^{-7}$. La fase $\theta=k\Delta x=2\pi/64\approx 0.098$ no es diminuta, así que el paso de ángulo pequeño se queda corto al final de la tabla. Y $\nu=0.5$ maximiza $\nu(1-\nu)$ en 0.25. Esta columna es el decaimiento en el CFL donde un paso duele más.

AP frente a AA#

No son la misma pregunta. AP pregunta si el límite $\varepsilon\to 0$ del esquema discreto es una discretización válida de la ecuación límite, con una restricción de estabilidad independiente de $\varepsilon$. AA pregunta algo más estrecho. En una malla finita, con $\varepsilon$ pequeño pero no nulo, ¿la solución se queda cerca del espacio well-prepared? El upwind falla las dos. El implícito IMEX pasa las dos. La bifurcación es la discretización temporal, no el orden espacial.

Incluso en el peor ν, la columna AP se queda en 1.00#

Cada entrada AP de la tabla es 1.00. El paso implícito quita la viscosidad upwind, así que el exponente $0.154/\varepsilon$ no multiplica la amplitud. $\nu=0.5$ es el CFL donde un paso decae más, y tras $M=128/\varepsilon$ pasos el valor AP sigue en 1.00. Para este sistema lineal de ondas, tratar solo $L/\varepsilon$ de forma implícita es la decisión correcta.

Esa decisión no viaja a todo sistema de bajo Mach. Hay ecuaciones en las que la acústica ya es implícita y el paso de tiempo vuelve a encogerse. El sitio donde un segundo $1/M^2$ se esconde en la velocidad relativa está en la rigidez que queda en la entalpía.