1950 — Cuando von Neumann difuminó el choque a propósito: viscosidad artificial y desacople par-impar

Un diente de sierra de densidad 1 y 2 sigue dando fuerza nula: el estencil de presión a dos celdas no ve el salto.

Los Álamos, 1950. John von Neumann y Robert Richtmyer publicaron un artículo breve en el Journal of Applied Physics: "A Method for the Numerical Calculation of Hydrodynamic Shocks". Si el choque no cabe en una celda, se difumina antes. Esa era la línea.

El tubo es el de Sod, membrana retirada en $t = 0$ : $\rho = 1, p = 1$ a la izquierda, $\rho = 0.125, p = 0.1$ a la derecha. Al subir el deslizador ξ de 0 a 5, el diente de sierra tras el choque se apaga y el salto se reparte en unas pocas celdas.

ξ artificial viscosity0.0

Sod shock tube. ξ=0: post-shock oscillations. ξ≈2-3: shock smeared over a few cells, wiggles damped.

Con ξ = 0 la espalda del choque es una sierra, cerca de ξ = 2 las oscilaciones mueren y el salto queda en dos o tres celdas, y con ξ = 5 el plateau se ondula. El $\xi \approx 2$ - $3$ recomendado por von Neumann queda a la vista como el punto dulce del compromiso.

Por qué funciona la viscosidad falsa, cuánto espesor deja y el modo par-impar que no cura: ese es el alcance de la nota. Al final, NumPy atrapa un tubo a Mach 2.

En una malla, los choques oscilan#

Las ecuaciones de Euler ideal no llevan viscosidad, así que un choque real se comprime hasta unos pocos recorridos libres medios. Una celda CFD es aproximadamente un millón de veces más ancha. La densidad y la presión deben duplicarse dentro de una celda.

Si se pasa esa discontinuidad por casi cualquier esquema conservativo, aparece un fenómeno único: una oscilación tipo diente de sierra detrás del frente, las famosas wiggles. Von Neumann ya las conocía desde 1944 en el Proyecto Manhattan. Cuando la región post-choque oscilaba, el cálculo del rendimiento del TNT salía mal.

La causa es directa. El upwind de primer orden tiene una difusión numérica fuerte que atrapa el choque pero suaviza todo lo demás. Los esquemas de segundo orden son precisos pero suenan junto a las discontinuidades, un fenómeno de Gibbs numérico. La naturaleza fabrica la entropía faltante con viscosidad molecular dentro del choque. La malla no tiene ese canal. Resolver el salto dentro de una celda como problema de Riemann es el tubo de choque que Godunov pone en cada celda. La viscosidad artificial hace lo contrario: estira el salto unas tres celdas para que las oscilaciones mueran.

La única línea de 1950#

La respuesta de von Neumann y Richtmyer fue imitar a la naturaleza. La viscosidad molecular no se puede resolver, así que se añade una presión de masa falsa a escala de celda.

Q_i = \begin{cases} \xi^2 \, (u_{i+1} - u_{i-1})^2 \, \rho_i & \text{if } u_{i+1} \le u_{i-1} \\ 0 & \text{otherwise} \end{cases}

$Q_i$ es una presión de masa artificial en la celda $i$ . $\xi$ es un parámetro de ajuste que controla sobre cuántas celdas se extiende el choque (típicamente 2-3). $u_{i+1} - u_{i-1}$ es la diferencia de velocidad. $\rho_i$ es la densidad local. Sustituir cada $P$ en momento y energía por $P + Q$ . Listo.

$Q$ va como $(\Delta u)^2$ , así que prácticamente vale cero en flujo suave y crece solo cerca de los choques. La condición es $u_{i+1} \le u_{i-1}$ : se activa solo cuando la velocidad converge localmente. Los abanicos de expansión y las ondas acústicas quedan intactos. Multiplicar por $\rho_i$ da un empuje mayor a medios más densos.

La ξ escrita y el 0,25 del código no dan el mismo espesor#

La forma corta es $Q = \xi^2 (\Delta u)^2 \rho$ . Solo se enciende si $\Delta u < 0$ . El listado multiplica por 0.25 * xi**2 * du**2 * rho. Ese du salta una celda. La pendiente de una sola celda mete $2\Delta x$ en el denominador, y el cuadrado aporta $1/4$ . El 0,25 es esa diferencia centrada. La fórmula de arriba se queda con $\xi^2$ y pierde el factor.

Con $\xi = 2,5$ , $\xi^2 = 6,25$ . El código multiplica de hecho por $0,25 \times 6,25 = 1,5625$ . Sin el 0,25, el mismo $\xi$ deja $Q$ cuatro veces mayor. Como $Q$ escala con $\xi^2$ , un factor cuatro equivale a duplicar $\xi$ . Un choque que se calma en el deslizador con $\xi = 2$ corresponde, en una copia literal de la fórmula, a $\xi = 4$ . Con $\xi = 5$ el plateau ya se ondula. Meter $\xi = 2,5$ sin el 0,25 cae de ese lado romo.

El listado fija 200 celdas, un tiempo de 0,2 y un CFL de 0,45. La frase bajo el código solo dice que el pico de densidad queda dentro de unos pocos por ciento del valor 2,667 de Rankine–Hugoniot para Mach 2. Contrastar la fórmula sin el 0,25 con esa frase mide otro espesor de choque.

Siguiendo el algoritmo en NumPy#

Cincuenta líneas guardan el núcleo. Advección donor-cell más el término de viscosidad artificial, nada más.

import numpy as np
 
def vnr_shock_solver(N=200, T=0.2, xi=2.5, gamma=1.4):
    # Condicion inicial: tubo de choque de Sod
    dx = 1.0 / N
    rho = np.where(np.arange(N) < N // 2, 1.0, 0.125)
    mom = np.zeros(N)
    p   = np.where(np.arange(N) < N // 2, 1.0, 0.1)
    E   = p / (gamma - 1)
    t = 0.0
    while t < T:
        u = mom / rho
        c = np.sqrt(gamma * p / np.maximum(rho, 1e-9))
        dt = 0.45 * dx / np.max(np.abs(u) + c)
        # Q de von Neumann-Richtmyer
        du = np.zeros(N)
        du[1:-1] = u[2:] - u[:-2]
        Q = np.where(du < 0, 0.25 * xi**2 * du**2 * rho, 0.0)
        # Adveccion donor-cell de rho, mom, E
        uf = 0.5 * (u[:-1] + u[1:])
        for q_arr, q_new in [(rho, 'rho'), (mom, 'mom'), (E, 'E')]:
            flux = np.where(uf > 0, q_arr[:-1] * uf, q_arr[1:] * uf)
            q_arr[1:-1] -= dt / dx * (flux[1:] - flux[:-1])
        # Actualizacion de presion y fuerza de cuerpo P+Q
        p = (gamma - 1) * (E - 0.5 * rho * (mom / rho) ** 2)
        Ptot = p + Q
        mom[1:-1] -= dt / (2 * dx) * (Ptot[2:] - Ptot[:-2])
        E[1:-1]   -= dt / (2 * dx) * (Ptot[2:] * u[2:] - Ptot[:-2] * u[:-2])
        t += dt
    return rho, mom, E
 
rho, mom, E = vnr_shock_solver(xi=2.5)
print(f"shock peak rho ≈ {rho.max():.3f}   (Rankine-Hugoniot exact ≈ 2.667 for Mach 2)")

Al correrlo, el pico de densidad cae dentro del uno por ciento del valor de Rankine-Hugoniot. El mismo código con ξ = 0 deja que la oscilación empuje el máximo y el mínimo asimétricamente lejos del plateau real.

Desacople par-impar: la otra trampa#

Atrapar el choque no cierra la historia. La misma malla colocalizada (todas las variables en el centro de celda) esconde otra patología que von Neumann notó. Basta un fluido en reposo, con energía térmica específica constante 1 y densidad en diente de sierra.

\rho_i = \begin{cases} 1 & i = 1, 3, 5, \dots \\ 2 & i = 2, 4, 6, \dots \end{cases}

La presión $p_i = \rho_i$ (caso isotérmico), de modo que la celda $i = 2$ está al doble de presión que sus vecinas. Debería relajarse al instante. Pero la actualización del momento en la celda $i$ es

\frac{\partial (\rho u)_i}{\partial t} = -\frac{F_{i+1/2} - F_{i-1/2}}{\Delta x} = -\frac{p_{i+1} - p_{i-1}}{2\Delta x} = 0

porque $p_{i+1}$ y $p_{i-1}$ vienen de la misma sub-malla de paridad. La malla se ha partido en dos sub-mallas independientes que no se ven entre sí. El modo diente de sierra no crece ni decae.

Initial: ρ alternates 1↔2 (isothermal, p=ρ). Collocated: ∂p/∂x at cell center cancels — sawtooth survives. Staggered: face gradient sees the jump and damps it.

Arriba, el modo pasa de colocalizado a desplazado (staggered). Una malla desplazada pone el momento en las caras de celda. La diferencia de presión en la cara $i+1/2$ es $p_{i+1} - p_i$ : las paridades par e impar se encuentran a través de una sola celda. El diente de sierra cae dentro del alcance de la primera derivada y se amortigua de inmediato.

Qué se ganó y qué se perdió#

El truco de 1950 mostró dos cosas a la vez. Los límites de la malla se pueden remendar por imitación. Y cada imitación crea su propia patología de malla. Ambas lecciones siguen vivas en el CFD moderno.

Códigos de astrofísica como ZEUS, Athena y PLUTO siguen usando mallas desplazadas con viscosidad artificial. OpenFOAM y ANSYS Fluent toman mallas colocalizadas y les ponen encima la interpolación de Rhie-Chow para esquivar el mismo desacople. Dos rutas alrededor del mismo problema.

Un diente de sierra detrás de un choque no siempre es un bug. Es una patología de malla conocida desde 1944.

Quitar el 0,25 con ξ = 2,5 cae del lado ξ = 5 del deslizador#

Sin el 0,25, meter $\xi = 2,5$ en la fórmula escrita iguala la $Q$ que el listado produce con $\xi = 5$ . En el deslizador ese punto ya es un choque romo. El rango $\xi \approx 2$ - $3$ pertenece al código que todavía lleva el 0,25. Trasladar el número sin el factor cambia el relato del espesor.

Que las celdas pares e impares no se vean es una enfermedad distinta de la viscosidad artificial. Poner la velocidad en la cara hace que la presión mire al vecino de una celda. Esa disposición fallando por un motivo distinto del diente de sierra post-choque está en la restricción de rotacional rota en una malla desplazada.