El hombre que borró la presión — el método de proyección de Chorin y el paso fraccionado

En 1967, Alexandre Chorin borró la presión de las ecuaciones, al menos por un instante. La presión es el término más incómodo al resolver las ecuaciones de Navier–Stokes incompresibles. No tiene derivada temporal. Como la densidad es constante, tampoco se puede recuperar la presión con una ecuación de estado. La respuesta de Chorin fue audaz: avanzar la velocidad ignorando la presión y luego forzar el resultado de vuelta a un estado de divergencia nula (divergence-free). Este artículo recorre ese paso fraccionado —conocido normalmente como método de proyección— partiendo de un teorema de una sola línea, la descomposición de Helmholtz–Hodge, hasta llegar a un solver 2D que de verdad corre. Al final quedará claro por qué la mayor parte del tiempo de ejecución de un código incompresible se lo come la ecuación de Poisson de presión.

La presión no tiene derivada temporal#

Al escribir las ecuaciones de Navier–Stokes incompresibles, el problema salta a la vista.

\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u}, \qquad \nabla\cdot\mathbf{u} = 0

Aquí $\mathbf{u}$ es la velocidad, $p$ la presión y $\nu$ la viscosidad cinemática (la difusividad del momento). La ecuación de momento indica cómo evoluciona $\mathbf{u}$ en el tiempo. Pero la segunda ecuación, $\nabla\cdot\mathbf{u}=0$ , no es una ecuación de evolución. Es una restricción que debe cumplirse en todo instante.

En flujo compresible, la ecuación de continuidad evoluciona la densidad, y la densidad fija la presión a través de una ecuación de estado. En el límite incompresible esa cadena se rompe. La presión no "evoluciona" en el tiempo. Es un multiplicador de Lagrange (una incógnita que impone la restricción) cuya única tarea es mantener la velocidad con divergencia nula en cada instante. Por eso intentar avanzar la presión en el tiempo es un esfuerzo inútil.

Todo campo de velocidad se parte en dos — Helmholtz–Hodge#

La salida está en un viejo teorema del cálculo vectorial. En un dominio razonable, cualquier campo vectorial $\mathbf{w}$ se descompone de forma única en una parte con divergencia nula más una parte de gradiente.

\mathbf{w} = \mathbf{u} + \nabla\phi, \qquad \nabla\cdot\mathbf{u} = 0

$\mathbf{u}$ es la componente rotacional (solenoidal) y $\nabla\phi$ la compresible (de gradiente). La operación que extrae solo $\mathbf{u}$ es el operador de proyección $\mathbb{P}$ . La receta es simple. Al tomar la divergencia de la ecuación anterior, como $\nabla\cdot\mathbf{u}=0$ , queda

\nabla^2\phi = \nabla\cdot\mathbf{w}

Se resuelve esta ecuación de Poisson para $\phi$ , se resta $\nabla\phi$ del campo original y solo sobrevive la parte de divergencia nula.

\mathbf{u} = \mathbb{P}\mathbf{w} = \mathbf{w} - \nabla\phi

Conviene experimentar con los dos paneles de abajo. El de la izquierda muestra un remolino más una fuente radial; el de la derecha muestra el mismo campo tras la proyección.

source strength6.0max|∇·u*| = 0.000 → 0.0000

Red = positive divergence (source), blue = negative (sink), dark = zero. Raise the source and the left panel lights up; the right panel stays dark — projection strips the compressible part and keeps only the swirl.

Al subir el source strength de 0 a 14, el panel izquierdo se tiñe de rojo y azul (su divergencia crece), mientras el derecho permanece oscuro. La proyección ha eliminado toda la componente compresible y ha conservado solo el remolino.

Predecir y luego proyectar — el paso fraccionado#

El algoritmo de Chorin coloca esta descomposición directamente en la integración temporal. Un paso se divide en dos.

Predictor — avanzar la velocidad sin presión, obteniendo una velocidad intermedia $\mathbf{u}^*$ .

\frac{\mathbf{u}^* - \mathbf{u}^n}{\Delta t} = -(\mathbf{u}^n\cdot\nabla)\mathbf{u}^n + \nu\nabla^2\mathbf{u}^n

$\mathbf{u}^*$ no tiene divergencia nula; no sorprende, porque se ignoró la presión.

Corrector — proyectar $\mathbf{u}^*$ de vuelta a divergencia nula. Se resuelve la ecuación de Poisson de presión,

\nabla^2\phi = \frac{1}{\Delta t}\nabla\cdot\mathbf{u}^*

y se resta el gradiente para construir la velocidad del paso siguiente.

\mathbf{u}^{n+1} = \mathbf{u}^* - \Delta t\,\nabla\phi

Aquí $\phi$ hace el papel de la presión ( $p^{n+1}\approx\rho\phi$ ). La presión no es algo que se evoluciona, sino un valor que se vuelve a resolver en cada paso para satisfacer la restricción. Esas dos líneas son todo el método de proyección.

La ecuación de Poisson de presión y sus trampas de contorno#

Aquí empieza el dolor práctico. El predictor es barato. El corrector, en cambio, resuelve una ecuación de Poisson en cada paso. En una malla grande, ese solver elíptico se lleva entre el 70% y el 90% del coste total. Por eso GMRES y multigrid son el corazón de un código incompresible.

Las condiciones de contorno también están llenas de trampas. La condición sobre $\phi$ es de Neumann en las paredes ( $\partial\phi/\partial n = 0$ ). Un problema de Neumann queda definido salvo una constante, y para que exista solución debe cumplir una condición de compatibilidad: $\int \nabla\cdot\mathbf{u}^*\,dV = 0$ en todo el dominio. Si los caudales de entrada y salida no coinciden, la ecuación de Poisson simplemente no tiene solución. Cuando aparezca un NaN, conviene sospechar primero del caudal de contorno.

En una malla colocalizada (collocated, no escalonada), poner presión y velocidad en el mismo punto provoca oscilaciones de tablero de ajedrez. Hay que usar una malla escalonada o frenarlas con la interpolación de Rhie–Chow. Como nota al margen: los métodos tipo SIMPLE repiten esta misma idea de proyección varias veces por paso para converger al estado estacionario, mientras que el paso fraccionado de aquí sigue un flujo no estacionario con una sola proyección por paso, conservando precisión temporal.

Python: enrollando una capa de cizalla doble con una FFT#

En contornos periódicos, la Poisson de presión se resuelve de un tirón con una FFT, porque en el espacio de Fourier el laplaciano se convierte en una multiplicación por $-|\mathbf{k}|^2$ . Con una capa de cizalla doble (dos chorros entrelazados) como condición inicial, el campo se enrolla en vórtices de Kelvin–Helmholtz, el problema clásico de verificación del método de proyección.

import numpy as np
 
N = 128                # malla (por lado)
h = 1.0 / N            # espaciado de malla
nu = 5e-3              # viscosidad cinematica -> Re = U L / nu ~ 200
dt = 2e-3              # dentro del rango estable de difusion/adveccion
steps = 3000
 
x = (np.arange(N) + 0.5) * h
X, Y = np.meshgrid(x, x, indexing='ij')
 
# capa de cizalla doble + pequena perturbacion
rho, delta = 1.0 / 30.0, 0.05
u = np.where(Y <= 0.5, np.tanh((Y - 0.25) / rho), np.tanh((0.75 - Y) / rho))
v = delta * np.sin(2 * np.pi * X)
 
# numeros de onda para la Poisson por FFT
k = 2 * np.pi * np.fft.fftfreq(N, d=h)
KX, KY = np.meshgrid(k, k, indexing='ij')
K2 = KX**2 + KY**2
K2[0, 0] = 1.0         # evitar dividir el modo medio por cero
 
def divergence(a, b):
    dadx = (np.roll(a, -1, 0) - np.roll(a, 1, 0)) / (2 * h)
    dbdy = (np.roll(b, -1, 1) - np.roll(b, 1, 1)) / (2 * h)
    return dadx + dbdy
 
def projection_correct(a, b):
    # Laplacian(phi) = div  ->  a <- a - grad(phi)
    phi_hat = np.fft.fft2(divergence(a, b)) / (-K2)
    phi_hat[0, 0] = 0.0
    phi = np.real(np.fft.ifft2(phi_hat))
    dpx = (np.roll(phi, -1, 0) - np.roll(phi, 1, 0)) / (2 * h)
    dpy = (np.roll(phi, -1, 1) - np.roll(phi, 1, 1)) / (2 * h)
    return a - dpx, b - dpy
 
def predictor(a, b):
    # adveccion (diferencia central) + difusion, un paso de Euler explicito
    ax = (np.roll(a, -1, 0) - np.roll(a, 1, 0)) / (2 * h)
    ay = (np.roll(a, -1, 1) - np.roll(a, 1, 1)) / (2 * h)
    bx = (np.roll(b, -1, 0) - np.roll(b, 1, 0)) / (2 * h)
    by = (np.roll(b, -1, 1) - np.roll(b, 1, 1)) / (2 * h)
    lap = lambda f: (np.roll(f, -1, 0) + np.roll(f, 1, 0)
                     + np.roll(f, -1, 1) + np.roll(f, 1, 1) - 4 * f) / h**2
    astar = a + dt * (-(a * ax + b * ay) + nu * lap(a))
    bstar = b + dt * (-(a * bx + b * by) + nu * lap(b))
    return astar, bstar
 
for n in range(steps):
    us, vs = predictor(u, v)                 # predecir: velocidad intermedia u*
    d_before = np.abs(divergence(us, vs)).max()
    u, v = projection_correct(us, vs)        # proyectar: eliminar la divergencia
    if n % 500 == 0:
        d_after = np.abs(divergence(u, v)).max()
        print(f"step {n:4d}  |div u*|={d_before:.2e} -> |div u|={d_after:.2e}")

Al mirar la salida: en cada paso $|\nabla\cdot\mathbf{u}^*|$ es del orden de $10^{-2}$ , pero tras la proyección $|\nabla\cdot\mathbf{u}|$ cae hasta $10^{-13}$ . La proyección por FFT ha aniquilado la divergencia hasta la precisión de máquina.

Qué pasa cuando se apaga la proyección#

Sin el corrector, el campo de velocidad acumula un poco de divergencia en cada paso. La masa deja de conservarse, los vórtices se emborronan y el campo pronto se convierte en ruido que llena la malla. Conviene verlo en la simulación de abajo.

dt0.60max|∇·u| = 0.000

Double shear layer rolling up. Red/blue = vorticity sign. Turn projection OFF and max|∇·u| climbs while the vortices dissolve into noise.

Con Projection en ON, la capa de cizalla se enrolla en un par limpio de vórtices. En el instante en que se cambia a OFF, la lectura max|∇·u| de la esquina superior se dispara y el patrón de vórtices se derrumba. Al subir dt, la advección se vuelve más agresiva y el colapso llega antes. Ese único interruptor es la respuesta más corta a "¿por qué necesitamos la proyección?".

Un resumen para quien no volverá a leer esto#

En flujo incompresible, la presión no es una variable que se evoluciona, sino un multiplicador de Lagrange que impone $\nabla\cdot\mathbf{u}=0$ en cada paso.
El método de proyección tiene dos etapas: predecir sin presión ( $\mathbf{u}^*$ ) y luego resolver una ecuación de Poisson para $\phi$ y restar su gradiente.
La mayor parte del coste vive en el solver de Poisson de presión, y la condición de compatibilidad de Neumann y el tablero de ajedrez son sus trampas clásicas.