Cuando la presión oscila en tablero de ajedrez — SIMPLE y Rhie-Chow

El campo de velocidad se veía bien. El de presión, no. Cada celda saltaba contra su vecina. Un tablero de ajedrez (checkerboard). Los residuos bajaban. La presión seguía en diente de sierra. No era un error de código. Quedó sellado al guardar presión y velocidad en los mismos puntos.

Al cambiar el control entre Naive y Rhie-Chow, el diente de sierra sembrado o se queda quieto o se asienta en una curva suave.

En Naive la amplitud no baja. En Rhie-Chow el mismo tablero se asienta. Reseed checkerboard vuelve a sembrar el diente.

Por qué el diente sobrevive, y cómo la continuidad se vuelve una ecuación de presión, viene después. El estiramiento compresible del mismo solver también queda más adelante.

La continuidad no ve la presión de al lado#

En una malla colocada (velocidad y presión almacenadas en el mismo centro de celda), la ecuación de cantidad de movimiento semidiscreta se escribe

donde $a_P$ es el coeficiente diagonal, $\mathbf{H}(\mathbf{u})$ reúne las contribuciones vecinas, $V_P$ es el volumen de celda y $(\nabla p)_P$ es el gradiente de presión en el centro de celda.

El problema es ese gradiente en el centro. Con diferencias centradas, $(\partial p/\partial x)_P = (p_E - p_W)/2\Delta x$. Nótese que $p_P$ no aparece. Una celda nunca ve su propia presión: solo la de los vecinos a dos celdas de distancia. Así, las celdas pares e impares se desacoplan (desacoplamiento par-impar). La presión en diente de sierra $(+,-,+,-,\dots)$ queda en el espacio nulo de las ecuaciones discretas y no produce residuo alguno.

Rhie–Chow: construir la velocidad de cara desde la cantidad de movimiento#

La cura es dejar de promediar velocidades de celda hacia la cara. Se reaplica la ecuación de cantidad de movimiento en la propia cara, de modo que la diferencia de presión vecina a través de esa cara entre directamente.

La barra superior es interpolación lineal. El término de atrás es el gradiente de presión recalculado en la cara. Ahora la velocidad de cara queda ligada a $p_E - p_P$, la diferencia entre celdas adyacentes. El modo diente de sierra ya no tiene dónde esconderse.

Mientras el salto de cara $p_E-p_P$ no sea cero, el diente queda atrapado en ese término. Otra forma de sacar la presión de la cantidad de movimiento es el método de proyección de Chorin.

SIMPLE: convertir la continuidad en una ecuación de presión#

La dificultad de fondo del flujo incompresible es que no existe una ecuación que resuelva la presión de forma directa. SIMPLE (Semi-Implicit Method for Pressure Linked Equations) convierte la continuidad $\nabla\cdot\mathbf{u}=0$ en una. Al sustituir la velocidad de cara de Rhie-Chow en la continuidad se obtiene una ecuación elíptica (de Poisson) para la presión.

El lado izquierdo es un laplaciano (forma de difusión); el derecho es la divergencia de la velocidad predicha. El procedimiento es predictor-corrector.

1. Resolver cantidad de movimiento con un $p^*$ supuesto para predecir $\mathbf{u}^*$.
2. Resolver la ecuación de presión anterior para actualizar $p$.
3. Corregir las velocidades de cara y de celda con el nuevo gradiente de presión para que se cumpla la continuidad.
4. Repetir 1–3 hasta converger.

PISO añade dos o tres pasos correctores más por iteración y se usa en cálculos transitorios. PIMPLE anida un bucle SIMPLE externo y un bucle PISO interno para pasos de tiempo grandes.

Sin subrelajación, revienta#

Hay una trampa. Si se reemplaza la presión por completo en cada paso ($p \leftarrow p^*$), el acoplamiento sobrecorrige y diverge. SIMPLE solo deja pasar una fracción de la actualización.

$\alpha_p$ es el factor de subrelajación de presión (de 0 a 1). La velocidad tiene su propio $\alpha_u$. Demasiado pequeño y la convergencia se arrastra; demasiado grande y oscila hasta estallar. Una regla práctica habitual es $\alpha_p \approx 0.3$ y $\alpha_u \approx 0.7$, elegidos para que ambos sumen alrededor de uno.

En el laboratorio de abajo conviene cambiar el factor de relajación por cuenta propia. Resuelve la ecuación de corrección de presión con un barrido de Gauss-Seidel relajado.

Al bajar el factor a 0.3, la curva de residuo desciende con lentitud. Cerca de 1.5 cae más rápido. Al acercarlo a 2, la sobrecorrección se acumula en cada barrido y el residuo vuelve a saltar hacia arriba. El $\alpha_p$ de SIMPLE vive justo sobre este mismo equilibrio.

Un solo solver para toda velocidad — el alcance compresible#

El verdadero atractivo del método basado en presión es que no le importa el número de Mach. En flujo compresible la continuidad se vuelve una ecuación de densidad, y una ecuación de estado (EOS) liga la densidad a la presión. Al escribir $\psi \equiv \partial\rho/\partial p$ (la compresibilidad, $=1/c^2$), aparece un término temporal en la ecuación de presión.

Esta ecuación tiene carácter convectivo y difusivo a la vez. Cuando $M\to 0$, $\psi$ crece y el término temporal deja de dominar, de modo que la presión se resuelve de forma elíptica (acoplamiento global e instantáneo). Con $M$ grande se vuelve hiperbólica (las ondas acústicas viajan a velocidad finita). Los solvers basados en densidad se vuelven rígidos a bajo Mach porque el acoplamiento densidad-presión se debilita, pero el método basado en presión mantiene ese acoplamiento explícito mediante la EOS. Así un solo código abarca desde subsónico hasta supersónico.

Revivir y borrar el tablero en Python#

Las afirmaciones no bastan. Conviene ejecutar ambos esténciles en una malla periódica 1D y medir si el tablero sobrevive o desaparece.

import numpy as np
 
def poisson_sweep(p, f, naive):
    """Un barrido de Gauss-Seidel de -p'' = f en una malla periodica 1D."""
    N = len(p)
    for i in range(N):
        if naive:                       # interp lineal naive: esténcil desacoplado que salta una celda
            p[i] = 0.5 * (p[(i - 2) % N] + p[(i + 2) % N] + f[i])
        else:                           # Rhie-Chow: esténcil compacto de 3 puntos que acopla vecinos
            p[i] = 0.5 * (p[(i - 1) % N] + p[(i + 1) % N] + f[i])
    p -= p.mean()                       # la presion es libre salvo una constante -> se fija la media
    return p
 
def checkerboard_metric(p):
    """Tamano de la componente (+,-,+,-,...). Cero significa sin diente de sierra."""
    signs = (-1.0) ** np.arange(len(p))
    return abs(np.dot(p, signs)) / len(p)
 
N = 48
x = 2 * np.pi * np.arange(N) / N
f = np.sin(x) + 0.4 * np.sin(2 * x)
f -= f.mean()
 
for naive in (True, False):
    p = 0.8 * (-1.0) ** np.arange(N)    # inicializar como tablero de ajedrez
    for _ in range(4000):
        p = poisson_sweep(p, f, naive)
    tag = "naive linear" if naive else "Rhie-Chow  "
    print(f"{tag}  checkerboard = {checkerboard_metric(p):.2e}")
 
# naive linear   checkerboard = 8.00e-01   <- el diente de sierra permanece
# Rhie-Chow      checkerboard = 3.1e-16    <- desaparece hasta la precision de maquina

La misma fuente, el mismo valor inicial, el mismo número de iteraciones. Solo cambió el esténcil. Naive no quita el diente de sierra. Rhie-Chow lo borra. Es el mismo sitio que la figura de arriba.

La presión inicial es $0.8\times(-1)^i$. Si la métrica es el valor absoluto del producto interno de $p$ con esa fila de signos, dividido por el número de celdas, arranca en 0.8. Tras 4000 barridos en 48 celdas, naive sigue en $8.00\times 10^{-1}$. No se movió ni un dígito. Rhie-Chow queda en $3.1\times 10^{-16}$. Precisión de máquina.

El esténcil es la razón de que naive conserve 0.8. En un diente de sierra, el índice $i$ y $i\pm 2$ comparten signo. $p_i=0.5(p_{i-2}+p_{i+2}+f_i)$ no promedia ese modo hasta cero. El gradiente en el centro $(p_E-p_W)/2\Delta x$ hace lo mismo. $E$ y $W$ son ambos opuestos a $P$, así que $p_E-p_W=0$. Un diente de amplitud 0.8 tiene gradiente nulo en el centro. No aparece en el residuo. El salto de cara $p_E-p_P$ cambia de signo, y su tamaño es 1.6. Ahí es donde Rhie-Chow agarra. La parte suave de $f=\sin x+0.4\sin 2x$ la resuelven ambos esténciles. La única diferencia que queda es el diente.

Con alpha_p en 1 el corrector revienta primero#

Reemplazar la presión entera, $p\leftarrow p^*$, sobrecorrige el acoplamiento en un solo paso. La actualización solo deja pasar la fracción $\alpha_p$. El par habitual $\alpha_p\approx 0.3$, $\alpha_u\approx 0.7$ suma 1. En la curva del laboratorio, 0.3 se arrastra, cerca de 1.5 cae más rápido, y cerca de 2 el residuo vuelve a saltar en cada barrido. SIMPLE deja $\alpha_p$ entre 0 y 1, en el lado lento de esa curva.

Donde cambiar la relajación cambia la solución convergida es la velocidad de cara que recuerda la subrelajación. En una malla colocada, promediar velocidades de celda hacia la cara pierde el salto de presión vecino, y el diente vive fuera del residuo. Revivir el término temporal con $\psi=\partial\rho/\partial p$ viene después. Añadir compresibilidad a una ecuación de presión que todavía carga el diente solo acelera el sonido. El tablero se queda.