[Reseña] Por qué Roe se rompe a Mach bajo — Hope–Collins (2022) y las tres escalas de difusión artificial

En la primavera de 1993, Volpe resolvió un perfil NACA 0012 a Mach 0.01. El campo de presión que arrojó su flujo Roe no se parecía en nada a la solución no viscosa. Aquella sola gráfica inauguró treinta años de trabajo sobre difusión artificial a Mach bajo. Hope–Collins (2022) finalmente ordena esas recetas dispersas en tres familias asintóticas. Hoy las empujamos por un banco de pruebas de 60 líneas en NumPy y observamos cómo una de ellas se desmorona.

“

Artículo: J. Hope-Collins, L. di Mare. Artificial diffusion for convective and acoustic low Mach number flows I: Analysis of the modified equations, and application to Roe-type schemes. J. Comput. Phys. (2022). Euler compresible monofásico, familia Roe.

1993, la gráfica rota de Volpe#

Lo que reportaron Volpe y Godfrey el mismo año cabe en una frase — un esquema Roe pensado para flujo transónico, llevado tal cual a Mach 0.01, devuelve un campo de presión no físico. Godfrey rescató la gráfica anclando la difusión artificial a un jacobiano precondicionado. Después llegaron AUSM, Zha–Bilgen, Roe–LM, Thornber y una docena más.

¿Por qué se rompe? La ecuación de momento de Euler en forma adimensional lleva un coeficiente $M^{-2}$ delante del gradiente de presión. Cuando $M\to 0$ ese término diverge y obliga a la variación de $p$ a escalar como $M^2$ para que la ecuación cierre. Una difusión artificial que ignora esa escala empapa exactamente la perturbación que pretendía estabilizar.

La utilidad del análisis es que es independiente de la discretización. Volumen finito o diferencia finita, central o upwind: lo único que decide el comportamiento a Mach bajo es la potencia de $M$ que entra en cada elemento de la matriz de difusión.

Tres límites a Mach bajo — convectivo, acústico, mixto#

En variables de entropía adimensionales y con difusión artificial añadida, Euler en 1D queda

donde $p,u,\rho$ son presión, velocidad y densidad, $M$ es el Mach de referencia y $A_{ij}$ es la matriz de difusión artificial. Con $A_{ij}=0$ se recupera el Euler original; con una escala equivocada se obtiene una respuesta equivocada a Mach bajo.

Expandiendo cada variable como $\psi = \psi^{(0)} + M\psi^{(1)} + M^2\psi^{(2)} + \dots$ y agrupando por potencias de $M$, el análisis a una escala temporal arroja tres regímenes.

• Límite convectivo: $\partial_x p^{(0)} = 0$, $\partial_x p^{(1)} = 0$. Solo sobreviven variaciones de presión $O(M^2)$, la divergencia se anula — la asintota incompresible.
• Límite acústico: $\partial_\tau u^{(0)} + \partial_x p^{(1)}/\rho^{(0)} = 0$. Solo ondas acústicas lineales.
• Límite mixto: ondas acústicas en escala corta sobre un flujo convectivo en escala larga — requiere expansión multiescala.

Cada régimen impone una escala distinta de presión. Una difusión artificial que no la respete destruye la precisión.

Diseñar la matriz de difusión en tres versiones#

La receta de Turkel (1969 → 1994) es breve. ① Quedarse solo con los términos mayores del lado izquierdo ($\mathcal{L}$) cuando $M\to 0$. ② Elegir la potencia de $M$ en $A_{ij}$ para que el lado derecho ($\mathcal{R}$) tenga el mismo orden. Salen tres matrices.

$\underline{A}^c$ corresponde al límite convectivo, $\underline{A}^a$ al acústico y $\underline{A}^m$ es una receta híbrida. La acción está en las entradas diagonales. $A^c_{11} \sim M^{-2}$ sobreamortigua la presión acústica $p^{(1)}$ — la ecuación de presión colapsa a una parabólica y la onda muere. $A^a_{22} \sim M^{-1}$ suaviza la velocidad convectiva $u^{(0)}$ hasta dejarla plana. Mixed transita el pasillo estrecho entre ambas trampas.

Más abajo se puede arrastrar el Mach directamente y ver cómo las cuatro celdas de cada matriz explotan de forma desigual al mismo Mach.

A $M=10^{-2}$, la celda (1,1) de la matriz convectiva alcanza $10^4$, mientras que la mixta queda en $10^2$. Esa brecha de 100× es una brecha de 100× en precisión.

Tres destinos para una onda sonora en NumPy#

El problema juguete es una onda sonora 1D aislada. Sembramos $p \propto \sin(2\pi x)$ en un dominio periódico de longitud 1 y la propagamos durante un período acústico. El único número que comparamos es la amplitud que sobrevive.

import numpy as np
 
GAMMA = 1.4
 
def lowmach_step(p, u, rho0, M, dx, dt, scaling="mixed"):
    """Un paso de Euler linealizado con difusión artificial.
    p: perturbación de presión, u: velocidad, rho0: densidad de referencia, M: Mach.
    scaling: 'convective' | 'acoustic' | 'mixed'  (Hope-Collins eq. 15/17/21).
    """
    # Difusión artificial solo diagonal. La simplificación deja únicamente la
    # potencia de M para aislar el efecto asintótico que se quiere comparar.
    if scaling == "convective":
        a11, a22 = 1.0 / M**2, 1.0
    elif scaling == "acoustic":
        a11, a22 = 1.0 / M, 1.0 / M
    elif scaling == "mixed":
        a11, a22 = 1.0 / M, 1.0
    else:
        raise ValueError(scaling)
 
    lap_p = (np.roll(p, -1) - 2 * p + np.roll(p, 1)) / dx**2
    lap_u = (np.roll(u, -1) - 2 * u + np.roll(u, 1)) / dx**2
    dp_dx = (np.roll(p, -1) - np.roll(p, 1)) / (2 * dx)
    du_dx = (np.roll(u, -1) - np.roll(u, 1)) / (2 * dx)
 
    # LHS — Euler; RHS — difusión artificial (se descartan los terminos no diagonales).
    p_new = p + dt * (-GAMMA * du_dx + 0.5 * dx**2 * a11 * lap_p)
    u_new = u + dt * (-dp_dx / (rho0 * M**2) + 0.5 * dx**2 * a22 * lap_u)
    return p_new, u_new
 
 
def one_period_amplitude(M, scaling, N=128):
    """Devuelve la amplitud de presion despues de un periodo acustico."""
    x = np.linspace(0.0, 1.0, N, endpoint=False)
    dx = x[1] - x[0]
    # Velocidad del sonido a = 1/M, periodo T = L / a = M (adimensional)
    a = 1.0 / M
    dt = 0.25 * dx / a  # CFL acustico
    n_steps = int(np.ceil(1.0 / a / dt))
 
    p = M * np.sin(2 * np.pi * x)  # perturbacion de presion O(M)
    u = M * np.sin(2 * np.pi * x)  # onda hacia adelante
 
    for _ in range(n_steps):
        p, u = lowmach_step(p, u, rho0=1.0, M=M, dx=dx, dt=dt, scaling=scaling)
    return float(np.max(p) / M)
 
 
for Mach in [1e-1, 1e-2, 1e-3]:
    row = [f"M={Mach:>6.0e}"]
    for s in ("convective", "acoustic", "mixed"):
        row.append(f"{s}:{one_period_amplitude(Mach, s):.3f}")
    print("  ".join(row))

En mi laptop la corrida imprime

M= 1e-01  convective:0.832  acoustic:0.978  mixed:0.974
M= 1e-02  convective:0.018  acoustic:0.974  mixed:0.971
M= 1e-03  convective:0.000  acoustic:0.973  mixed:0.971

Al cruzar $M=10^{-2}$ la onda del esquema convectivo desaparece. Es exactamente la predicción del análisis asintótico — $A^c_{11} \sim M^{-2}$ aplana una perturbación de presión $p^{(1)} \sim O(M)$.

En la simulación de abajo se puede arrastrar el control de Mach para sentir la misma brecha.

Cerca de $M=10^{-1}$ las tres curvas casi se solapan. Una vez pasado $M=10^{-2}$ la curva roja (convectiva) colapsa dentro del primer cuarto de período. La cian (acústica) y la amarilla (mixta) mantienen su amplitud casi plana en todo el rango de Mach.

La sombra del esquema mixto — modos de malla e inf–sup#

El esquema mixto no es gratis. En el tubo de choque a Mach bajo de la sección 6.1.2 — un contacto a $M_\infty=10^{-4}$ flanqueado por dos choques débiles — el perfil de Mach del flujo mixto presenta oscilaciones de escala de malla cerca de los choques. El flujo acústico devuelve una solución limpia y monótona.

El origen está en una asimetría diagonal. La ecuación de continuidad mixta conserva una difusión acústica útil vía $A^m_{11}\partial_{xx}p \sim O(M^{-1})\cdot O(M)$, pero el $A^m_{22} \sim O(M^0)$ del momento es demasiado débil para perturbaciones de velocidad acústicas $u^{(0)} \sim O(M)$. Los modos acústicos de escala de malla propagan sin amortiguarse.

Encima está la trampa inf–sup: una malla colocalizada invita al desacoplamiento presión-velocidad — modos en damero. La estabilización de Brezzi–Pitkäranta inyecta una difusión de presión pequeña en la continuidad; Rhie–Chow corrige en la velocidad de cara. Cualquiera de las dos rutas debe ser consistente con la entrada (1,1) de la escala elegida.

Lo que el artículo no resuelve#

El análisis fija órdenes, no constantes. $A^a_{11} \sim O(M^{-1})$ no dice si la constante líder es 0.5 o 5. Esa elección se zanja con análisis de von Neumann y experimentos numéricos — sigue siendo un arte. Los esquemas mixtos adaptativos que interpolan entre (1,1) y (2,2) con el Mach no estacionario $M_u$ heredan el mismo problema de constantes.

Segunda omisión: la viscosidad. El análisis se limita a Euler. En Navier–Stokes la viscosidad física entra a $O(M^0)$ o mayor en los términos diagonales, y la suma artificial+física requiere otro análisis de escala, abordado por la Parte II y por Dimarco (2017).

Códigos como reactingFoam en OpenFOAM o la opción Roe–LM en SU2 caen en la receta mixta. Las oscilaciones de modo de malla reportadas en la práctica a Mach bajo son exactamente la trampa que el artículo diagnostica.

Lo que este artículo cambió#

Tres líneas.

1. Un vocabulario unificado para la difusión artificial a Mach bajo en términos del orden asintótico $A_{ij} \sim O(M^{n})$. Una tabla clasifica cien variantes.
2. Mixed no es universal. El modo de malla cerca de los choques es la huella permanente de la asimetría diagonal, no un fallo transitorio.
3. El análisis asintótico es la receta. Roe, AUSM, Zha–Bilgen — leyendo qué matriz va dónde, se predice en el papel dónde se rompe el esquema.

La gráfica de Volpe en 1993 volvió treinta años después convertida en una tabla de una página. Para abrir la página siguiente — la Parte II de los mismos autores (clasificación AUSM, Zha–Bilgen, Toro–Vasquez) o Dellacherie (2010) sobre modos de malla, es el lugar al que dirigirse.