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

Un desajuste entre el flujo que entra y el que sale deja la Poisson de Neumann sin solución. Sin proyección, esa divergencia se queda en cada paso.

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.

Al subir la intensidad de la fuente en el panel izquierdo de 0 a 14, los colores de la divergencia se encienden, y la proyección de la derecha conserva solo el remolino.

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.

Con la fuente en 0 los dos paneles casi coinciden. En 14 solo el izquierdo se llena de rojo y azul, y el derecho sigue oscuro. La proyección ha quitado la parte compresible.

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

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. Hace falta una malla escalonada, o frenarlas con la interpolación de Rhie–Chow. Si esa velocidad de cara recuerda el factor de subrelajación, la solución convergida misma cambia. 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.

Un lado de 128 deja un espaciado de 1/128. Viscosidad 5e-3 y escalas unidad dan Re = 200, como dice el comentario. Tres mil pasos de 2e-3 integran un tiempo 6. El espesor 1/30 ronda 0.033, poco más de cuatro celdas. La perturbación 0.05 es el 5 por ciento del chorro, y el enrollamiento arranca ahí. La divergencia $10^{-2}$ frente a $10^{-13}$ son once cifras. Fijar el modo medio de Fourier en 0 cierra la constante de Neumann. Un desbalance entre entrada y salida no se arregla con ese denominador. El 70 al 90 por ciento del coste es esta matriz. Al predictor le queda el otro 10 al 30 por ciento.

Borrar una divergencia de $10^{-2}$ con paso 2e-3 pide al laplaciano de $\phi$ un valor cercano a 5. Esa es la corrección del gradiente de presión en un paso. Sin proyección se queda como divergencia, y 3000 pasos bastan para llenar la malla.

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?".

Si el caudal se desfasa una gota, la Poisson no tiene respuesta#

En flujo incompresible, la presión no es una variable que se evoluciona. Es un multiplicador de Lagrange que impone $\nabla\cdot\mathbf{u}=0$ en cada paso. La proyección construye $\mathbf{u}^*$ sin presión, resuelve una Poisson para $\phi$ y resta el gradiente. Ese orden es corto. La rotura está en el contorno. Si la entrada y la salida no coinciden, no hay condición de compatibilidad, y la matriz se queda con la libertad de una constante y no devuelve solución.

Conservar la precisión temporal con una proyección por paso significa no volver a resolver la convección dentro de ese paso. Subir solo el orden viscoso y empujar la convección afuera es la bifurcación de dónde la extrapolación de la convección en Gear agrandó el error. Hacer que la convección del predictor parezca más exacta puede mover la divergencia nula que la proyección acaba de fijar.