Por dónde se filtra la densidad negativa — un limitador de flujo que preserva la positividad para ondas de explosión

Se conserva la θ más grande que deja la densidad del paso siguiente sobre el piso; por encima de eso está la fuga.

Las 2:27 de la madrugada. Un tubo de choque Le Blanc escupió NaN a los 0.003 s. Con el CFL en 0.1 y la malla duplicada, las dos curvas se cayeron juntas. La causa era pequeña. En la cabeza de una rarefacción fuerte, la densidad de una celda cruzó por debajo de cero. Dos pasos después la velocidad del sonido era NaN. El paso siguiente se llevó el dominio. La mezcla θ de Hu, Adams y Shu (2013) hace retroceder el flujo de alto orden un solo escalón, y solo ahí.

Al bajar la densidad derecha ρ_R por debajo de 10⁻³, la curva permitida entra en la zona roja y la θ de esa cara cae bajo 0.2.

Per-face θ search · single rarefaction face · curve ρ_iⁿ⁺¹(θ) at fixed left face

log₁₀(ρ_R) = -2.00 (ρ_R = 1.00e-2)outflow velocity u_R = 4.00

The cyan curve is ρⁿ⁺¹ after one Euler step as a function of how much central-flux we mix into the right face. The red zone is where ρ has crossed the floor 10⁻³. The yellow dot is the largest θ that still lands above the floor — that is the θ the per-face limiter picks. Drop ρ_R below 10⁻³ or push u_R past 6 and θ collapses toward zero: the limiter falls back to first-order Lax–Friedrichs.

La zona roja es donde la densidad del paso siguiente rompe el piso. Si θ se detiene en ese borde, la cara queda en la práctica como Lax–Friedrichs de primer orden.

Por dónde se filtra la densidad negativa#

Un esquema de alto orden oscila entre estabilidad y precisión. El upwind de primer orden emborrona un choque a lo largo de unas cinco celdas. La diferencia central de segundo orden siembra un diente de sierra a ambos lados. Si la oscilación es leve, solo tiembla la presión. Cuando el cociente de densidad entre celdas vecinas pasa de 100:1 — Le Blanc, una explosión de Sedov, la estela de una tobera supersónica — cambia el modo de fallo.

La conservación de densidad en la celda $i$ queda así:

\rho_i^{n+1} = \rho_i^n - \frac{\Delta t}{\Delta x}\bigl( F_{i+1/2}^\rho - F_{i-1/2}^\rho \bigr)

$F^\rho = \rho u$ es el flujo de densidad. En una celda con $F_{i+1/2}^\rho \gg F_{i-1/2}^\rho$ , donde sale más masa de la que entra, $\rho_i^{n+1}$ puede volverse negativa. En Euler 1D esa celda es la primera de la expansión justo detrás de un choque fuerte.

La densidad negativa mata primero a la velocidad del sonido. Con $a = \sqrt{\gamma p / \rho}$ , si $\rho \to 0^-$ entonces $a^2 < 0$ . La siguiente evaluación del flujo devuelve NaN desde $\sqrt{\cdot}$ . El NaN de una celda alcanza a las dos vecinas. Dos pasos después el dominio termina.

Lax–Friedrichs — borroso pero nunca roto#

El flujo LF que Lax escribió en 1954 está fuera de esta trampa.

F_{i+1/2}^{LF} = \tfrac{1}{2}\bigl(F_L + F_R\bigr) - \tfrac{1}{2}\alpha\bigl(U_R - U_L\bigr)

$\alpha$ es la velocidad de señal máxima de las dos celdas, y $U = (\rho, \rho u, E)^T$ es el estado conservativo. El término que se resta es viscosidad artificial. Promedia a los vecinos y mata el diente de sierra. Es el mismo gesto que la viscosidad que von Neumann untó a propósito sobre un choque en 1950. Bajo CFL $\leq$ 1, Lax mostró que $\rho^{n+1} \geq 0$ , una forma reducida del teorema 1 de Hu et al. Choque rápido o débil, LF no se cae.

El costo es la precisión. Un choque de dos celdas pasa a cinco. Un contacto se difumina una celda más en cada paso. En producción nadie usa LF solo. El margen de seguridad sigue siendo el punto de partida del limitador.

La mezcla θ — retroceder una muesca cada vez#

Un flujo de alto orden $F_{i+1/2}^{HO}$ (WENO, MUSCL, central) es más nítido que LF y aun así puede crear densidad negativa. Hu, Adams y Shu atan los dos flujos en una línea.

F_{i+1/2} = F_{i+1/2}^{LF} + \theta_{i+1/2}\bigl( F_{i+1/2}^{HO} - F_{i+1/2}^{LF} \bigr)

$\theta \in [0, 1]$ es un peso por cara. $\theta=1$ es alto orden puro. $\theta=0$ es LF puro. El limitador hace una sola cosa: elegir la θ más grande que deje positiva la densidad siguiente de cada celda. No retrocede por todo el dominio. Un escalón.

El estado LF ya cumple $\rho^{LF}_i \geq \rho_{\min}$ , así que cada cara resuelve la θ más grande que mantiene positiva la actualización tras sumar $\Delta F = F^{HO} - F^{LF}$ . Una cara toca dos celdas. Ambas tienen que quedar sobre el piso, de modo que gana la cota menor.

Esa es la cota que bajó θ por debajo de 0.2 cuando bajó ρ_R. La cara se resuelve entonces como LF de primer orden.

Limitación en dos pasos — densidad primero, presión después#

Densidad positiva no basta. Presión negativa mata el solver del mismo modo, porque la velocidad del sonido se vuelve imaginaria. Hu et al. parten el limitador en dos pasadas.

Pasada de densidad. La fórmula de arriba fija $\theta^\rho$ en cada cara para que toda densidad siguiente quede en o sobre $\rho_{\min}$ .

Pasada de presión. Con el flujo de densidad congelado por $\theta^\rho$ , los flujos de momento y energía reciben una $\theta^p \leq \theta^\rho$ . La restricción de presión

p = (\gamma - 1)\bigl(E - \tfrac{1}{2}\rho u^2\bigr) \geq p_{\min}

es no lineal, cuadrática en general. Una elección conservadora, tomar la raíz menor para $\theta^p$ , cabe en una línea.

Un solo peso. $\theta = \min(\theta^\rho, \theta^p)$ mezcla todos los flujos conservativos. Una cara, una θ, y las variables conservadas siguen consistentes entre sí.

50 líneas de NumPy — θ automática en un tubo Le Blanc#

Le Blanc aquí es un tubo fuerte, $\rho_L=1$ , $p_L=100$ frente a $\rho_R=0.02$ , $p_R=0.05$ . Una rarefacción intensa viaja a la izquierda y la densidad de la primera celda se vuelve negativa si nada la guarda.

Al dividir los dos estados, el cociente de densidad es $1/0.02 = 50$ y el de presión es $100/0.05 = 2000$ . La línea de 100:1 de arriba era una regla práctica, y el salto de densidad de este par queda por debajo. La violencia está en la presión. Con $\gamma = 1.4$ , $a = \sqrt{\gamma p / \rho}$ da $\sqrt{140} \approx 11.83$ a la izquierda y $\sqrt{3.5} \approx 1.871$ a la derecha. El cociente es exactamente $\sqrt{40} \approx 6.32$ .

Los números de la malla dicen que el paso de tiempo no fue el culpable. $N = 200$ , así que $\Delta x = 1/200 = 0.005$ , y $\Delta t = 9.0 \times 10^{-5}$ . Entonces $\Delta t/\Delta x = 0.018$ . Un número de Courant acústico solo con el estado izquierdo es $a_L \Delta t/\Delta x \approx 0.213$ . Encaja con una corrida que ya murió con CFL 0.1. Desde $\rho^{n+1} = \rho^n - (\Delta t/\Delta x)\Delta F$ , si la diferencia de flujo supera el margen de densidad, un paso más corto igual perfora el piso.

En el estado derecho $\rho = 0.02$ , el margen hasta el piso $10^{-3}$ es $0.019$ . Para que θ se quede en 1, la corrección de densidad $|d|$ no puede pasar de $0.019/0.018 \approx 1.06$ . Con $|d| = 2$ la cota es $0.019/(0.018 \times 2) \approx 0.53$ . Si el margen se encoge a $5 \times 10^{-4}$ , el mismo $|d| = 1$ baja la cota a $0.0005/0.018 \approx 0.028$ . Esa cara es casi LF puro.

En pesos, la mezcla es $F = (1-\theta)F^{LF} + \theta F^{HO}$ . Con $\theta = 0.53$ la parte de alto orden es 53%. Con $\theta = 0.028$ queda en 2.8% y el resto es LF. La cota del listado es la misma desigualdad de la actualización: la corrección mueve la densidad en $(\Delta t/\Delta x)\,\theta\,|d|$ , y eso no puede pasar el margen $0.019$ , así que $\theta \le 0.019/(0.018\,|d|)$ . La igualdad en $\theta = 1$ cae cerca de $|d| \approx 1.06$ .

La fórmula de la presión se lee igual. $p = (\gamma - 1)(E - \tfrac{1}{2}\rho u^2)$ y $\gamma - 1 = 0.4$ . Con $E$ y $\rho$ fijos, una unidad de energía cinética recorta 0.4 de presión. El estado derecho $p = 0.05$ toca cero con una energía cinética de $0.05/0.4 = 0.125$ . El izquierdo $p = 100$ aguanta hasta $250$ . Cuidar solo el piso de densidad desplaza la fuga hacia una densidad positiva con velocidad del sonido imaginaria.

import numpy as np
 
GAMMA = 1.4
RHO_FLOOR = 1e-3
N, NSTEP = 200, 220
dx = 1.0 / N
dt = 9.0e-5
 
def cons_to_prim(U):
    rho, mu, E = U[0], U[1], U[2]   # conservativo -> primitivo
    u = mu / rho
    p = (GAMMA - 1.0) * (E - 0.5 * rho * u * u)
    return u, p
 
def euler_flux(U):
    u, p = cons_to_prim(U)
    return np.array([U[1], U[1] * u + p, u * (U[2] + p)])
 
def alpha_face(UL, UR):
    uL, pL = cons_to_prim(UL); uR, pR = cons_to_prim(UR)
    aL = np.sqrt(GAMMA * max(pL, 0.0) / UL[0])
    aR = np.sqrt(GAMMA * max(pR, 0.0) / UR[0])
    return max(abs(uL) + aL, abs(uR) + aR)
 
def density_theta(U, dF, dt_dx):
    """θ por cara: máximo peso que mantiene ρ por encima del piso."""
    theta = np.ones(N + 1)
    for f in range(1, N):
        d = dF[0, f]                          # anti-difusión de densidad
        if abs(d) < 1e-14: continue
        rL_room = U[0, f - 1] - RHO_FLOOR     # margen celda izquierda
        rR_room = U[0, f]     - RHO_FLOOR     # margen celda derecha
        cap = (rL_room if d > 0 else rR_room) / (dt_dx * abs(d))
        theta[f] = max(0.0, min(theta[f], cap))
    return theta

El controlador hace una pasada LF de primer orden, una pasada de antidifusión central y una pasada para θ. 220 pasos, 200 celdas. Con el limitador apagado, NaN cerca del paso 30. Encendido, la corrida termina.

Con el interruptor del limitador apagado, la curva roja arrastra la densidad por debajo de cero cerca del paso 30.

1D strong shock tube · ρ_L=1, p_L=100 / ρ_R=0.02, p_R=0.05 · central + θ·anti-diffusion · 200 cells

t-index = 80 / 220show limited

Red is the raw centered scheme — its rarefaction tail dives below zero around t ≈ 0.006 and the solver blows up. Cyan is the same scheme blended with Lax–Friedrichs by a per-face θ ≤ θ_max so that ρ stays above the floor 10⁻³. Drag the time slider to watch the two solutions separate.

Ese cruce es el colapso. La curva cian se aplana cerca de $\rho \approx 0.02$ en el mismo instante. La diferencia de precisión queda bajo el 5%. Una corrida vive. La otra muere.

Lo que cuesta el limitador#

Una cara con $\theta = 0$ se resuelve por LF de primer orden. Los frentes de choque y las cabezas de rarefacción se embotan ahí. Esa zona ocupa menos del 5% del dominio, así que el error global apenas se mueve. La θ media queda entre 0.95 y 1.0.

Fijar $\rho_{\min}$ en algo como $10^{-12}$ hace que el limitador despierte tarde. El mismo NaN regresa. En la práctica el piso anda cerca de $10^{-4}$ de la densidad media. Un Runge–Kutta de varias etapas tiene que repetir el limitador en cada etapa. Proteger solo la primera etapa deja que la siguiente diverja igual.

Un piso de 1e-12 despierta al limitador después del NaN#

$\rho_{\min} = 10^{-12}$ es adorno, no un piso. Frente a una densidad derecha de 0.02 casi todo el margen sigue libre, y la cota de arriba se queda cerca de 1. El límite se aprieta solo después de que la densidad ya cayó bajo $10^{-8}$ . Para entonces $\gamma p / \rho$ cambió de signo. El limitador retrocede a LF cuando la celda ya está muerta.

La θ de la cara vigila el promedio de la celda. Si ese promedio sigue positivo mientras un punto reconstruido dentro de la celda cruza a negativo, esta línea no ve el punto. Queda abierto si una reconstrucción positiva en el promedio y negativa en un punto necesita todavía un arreglo propio.

Referencia: X. Y. Hu, N. A. Adams, C.-W. Shu. Positivity-preserving method for high-order conservative schemes solving compressible Euler equations. J. Comput. Phys. 242 (2013) 169–180.