una calle de vórtices de Kármán (los remolinos alternos que se forman detrás de cualquier objeto contundente en un flujo rápido). Has visto este patrón detrás de los pilares de puentes, alas de aviones y chimeneas. Es uno de los fenómenos más reconocibles en dinámica de fluidos.
Lo generé sin resolver una sola ecuación de fluido.
En cambio, esta simulación del Método Lattice Boltzmann (LBM) trata el fluido como lo que realmente es a nivel microscópico: una distribución estadística de partículas que chocan entre sí. El comportamiento macroscópico (vórtices, gradientes de presión, amortiguación viscosa) es una consecuencia, no un insumo.
LBM ocupa una posición realmente inusual en la física computacional, a medio camino entre una simulación de dinámica molecular y un solucionador de Navier-Stokes, que opera a escala mesoscópica. En este artículo, lo derivaré de los primeros principios, lo implementaré en ~200 líneas de C++ y lo ejecutaré en la supercomputadora MareNostrum 5. El código fuente completo y una versión simplificada de Python/NumPy están en GitHub.
[insertar]https://www.youtube.com/watch?v=TGJ6W_Jrg2w[/embed]
Resumen rápido: LBM reemplaza el intratable PDE de Navier-Stokes con una regla de actualización mucho más simple en una red de velocidad discreta. Las ecuaciones NS se recuperan en el límite correcto, no se asumen. Esto hace que LBM sea teóricamente elegante, prácticamente rápido y vergonzosamente paralelizable. Funciona mejor en el régimen subsónico de bajo Mach y se vuelve numéricamente inestable a viscosidades muy bajas.
Por qué el enfoque obvio es doloroso
La ruta estándar hacia la simulación de fluidos pasa por las ecuaciones de Navier-Stokes, las leyes de conservación macroscópicas que gobiernan prácticamente todo lo que nos importa en dinámica de fluidos. Son analíticamente intratables excepto en un puñado de geometrías de juguete, por lo que en la práctica se discretizan en una malla y se resuelven numéricamente: métodos de volumen finito, diferencias finitas o elementos finitos.
Estos funcionan bien, pero conllevan un costo oculto que se complica en geometrías complejas. Para un límite curvo (un cilindro, la carrocería de un automóvil, un medio poroso), se necesita una malla ajustada a la carrocería que es costosa de generar y difícil de automatizar, o correcciones de plantilla especializadas cerca de la pared. Cada vez que cambia la geometría, es necesario regenerar la malla. Es un problema de ingeniería superpuesto al problema de física.
LBM evita esto casi por completo. La razón queda clara una vez que comprendes dónde comienza.
Un punto de partida diferente: la mecánica estadística
En lugar de preguntar “¿cuál es el campo de velocidades en este punto?”, LBM pregunta “¿cuántas partículas se mueven en cada dirección en este punto?”
El objeto central es la función de distribución de una sola partícula f(x, v, t). Le indica el número esperado de partículas cerca de la posición x que se mueven con velocidad v en el momento t. Los campos macroscópicos que realmente nos interesan son sólo momentos de f:
[
rho(mathbf{x},t) = mint f, d^3v qquad rho,mathbf{u}(mathbf{x},t) = mint mathbf{v},f, d^3v
]
En equilibrio global (densidad uniforme, sin flujo aparente), f toma la forma familiar de Maxwell-Boltzmann:
[
f^{eq}(mathbf{v}) = frac{rho}{(2pi RT)^{D/2}} exp!left(-frac{|mathbf{v} – mathbf{u}|^2}{2RT}right)
]
Lejos del equilibrio, f toma otra forma. La ecuación que rige su evolución es la ecuación de Boltzmann:
[
frac{partial f}{partial t} + mathbf{v}cdotnabla_mathbf{x} f = Omega[f]
]
El lado izquierdo es transporte puro: las partículas llevan su distribución de velocidades a medida que fluyen por el espacio. El lado derecho, Ω[f], es el operador de colisión que hace que f regrese a f^eq redistribuyendo partículas entre clases de velocidad. El problema es que Ω[f] es una integral de cinco dimensiones sobre todos los posibles compañeros de colisión y ángulos de dispersión, que codifica la física de colisión microscópica completa. Para cualquier simulación práctica, es completamente intratable.
El milagro de una sola línea: BGK
En 1954, Bhatnagar, Gross y Krook propusieron una simplificación brillantemente contundente. En lugar de modelar cada colisión en detalle, simplemente capture el efecto neto: las colisiones llevan a f hacia el equilibrio a una tasa proporcional a qué tan lejos del equilibrio se encuentra actualmente.
[
Omega[f] approx -frac{1}{tau}bigl(f – f^{eq}bigr)
]
Esta es la aproximación BGK y reduce una integral de cinco dimensiones intratable a una multiplicación escalar. El parámetro τ es el tiempo de relajación, es decir, el tiempo promedio que tarda el sistema en regresar al equilibrio local después de una perturbación. τ grande significa relajación lenta; pequeña τ significa rápido.
Lo que hace que esto sea más que un truco conveniente es lo que τ significa físicamente. Mediante la expansión de Chapman-Enskog (un análisis perturbativo sistemático) se puede demostrar que en el límite del número de Knudsen pequeño (muchas colisiones por camino libre medio), la ecuación de BGK recupera exactamente las ecuaciones de Navier-Stokes, con viscosidad cinemática:
[
nu = c_s^2!left(tau – frac{1}{2}right)
]
Tenga en cuenta que esta relación opera en unidades de red adimensionales; la conversión a unidades físicas requiere una calibración de longitud y escala de tiempo por separado. Pero dentro de la simulación, τ es toda la perilla física: gírela hacia arriba para obtener un flujo espeso y viscoso; empújelo hacia 0,5 para un flujo con un número de Reynolds alto (aunque los valores por debajo de ~0,55 corren el riesgo de inestabilidad numérica).
Fundamentalmente, esta es la razón por la que LBM es algo más que Navier-Stokes disfrazado. La conexión con NS se deriva, no se asume: estás resolviendo una ecuación cinética que probablemente converge con el comportamiento de NS en el límite apropiado. Este es un tema muy común en toda la mecánica estadística.
De continuo a discreto: la red D2Q9
La ecuación BGK todavía vive en un espacio de velocidad continuo. Para simularlo necesitamos un conjunto finito de velocidades. El enfoque es una expansión de Taylor de f^eq en potencias de la velocidad total u, válida cuando u es pequeña en comparación con la velocidad reticular del sonido (el supuesto de bajo Mach):
[
f_i^{eq} = rho, w_i!left[1 + 3(mathbf{c}_icdotmathbf{u}) + frac{9}{2}(mathbf{c}_icdotmathbf{u})^2 – frac{3}{2}|mathbf{u}|^2right]
]
La velocidad continua c se reemplaza por un conjunto finito de velocidades discretas c_i, cada una de las cuales lleva un peso de cuadratura w_i elegido para recuperar correctamente los momentos de velocidad. Para simulaciones 2D, la opción estándar es D2Q9: 2 dimensiones, 9 direcciones de velocidad: reposo, cuatro direcciones alineadas con los ejes y cuatro diagonales.
Los pesos se fijan por isotropía: 4/9 para el resto de partículas, 1/9 para las direcciones de los ejes, 1/36 para las diagonales. La velocidad reticular del sonido es cs = 1/√3. Todo lo demás en el algoritmo se deriva de estos nueve números.
El algoritmo: dos pasos, repetidos hasta el infinito
Con el equilibrio discreto en la mano, el algoritmo LBM completo se reduce a cuatro operaciones por paso de tiempo y se ejecuta en un bucle hasta la convergencia:
Recuperación macroscópica: calcule ρ y u en cada nodo a partir de las poblaciones actuales f_i Colisión: relaje cada f_i hacia f_i^eq a una velocidad de 1/τ (puramente local, sin comunicación con vecinos) Streaming: propague cada población al nodo vecino en su dirección de viaje Condiciones de contorno: reconstruya las poblaciones faltantes en los bordes del dominio y las paredes sólidas
El paso de colisión es donde ocurre toda la física. El paso del streaming es una contabilidad puramente cinemática. Su alternancia es la que produce un comportamiento fluido.
colisión vacía() { const doble omega = 1.0 / tau; #pragma omp paralelo para colapsar(2) programar(estático) for (int x = 0; x < Nx; x++) { for (int y = 0; y < Ny; y++) { if (solid[x][y]) continuar; doble feq[9]; calcularEquilibrio(rho[x][y], ux[x][y], uy[x][y], feq); for (int i = 0; i < 9; i++) f[x][y][i] -= omega * (f[x][y][i] – feq[i]); } } }
La implementación del streaming resulta de gran importancia para el rendimiento. El ingenuo patrón de inserción (donde cada subproceso en (x,y) escribe poblaciones salientes en nodos vecinos) provoca un intercambio falso: los subprocesos escriben en regiones de memoria adyacentes, martillando las mismas líneas de caché y forzando una invalidación constante. El siguiente patrón de extracción (reunión) soluciona este problema, porque cada subproceso lee de sus vecinos ascendentes y solo escribe en su propio nodo.
void stream() { #pragma omp paralelo para colapso(2) horario(estático) para (int x = 0; x < Nx; x++) { para (int y = 0; y < Ny; y++) { if (solid[x][y]) continuar; for (int i = 0; i < 9; i++) { const int xsrc = x – EIX[i]; const int ysrc = y – EIY[i]; const bool oob = (xsrc < 0 || xsrc >= Nx || ysrc < 0 || ysrc >= Ny); f_new[x][y][i] = (oob || solid[xsrc][ysrc])? f[x][y][OPP[i]] // rebote desde la pared: f[xsrc][ysrc][i]; // flujo desde el vecino ascendente } } } }
Debido a que la colisión es completamente local y la transmisión solo requiere lecturas del vecino más cercano, todo el algoritmo es vergonzosamente paralelo: cada nodo se actualiza de forma independiente, sin soluciones globales, sin sistemas implícitos y sin inversiones de matrices. Esto se asigna claramente a CPU, GPU y clústeres distribuidos de múltiples núcleos; Las implementaciones de LBM en CUDA logran de forma rutinaria un escalamiento casi lineal con un número de subprocesos hasta el límite del ancho de banda de la memoria.
Condiciones de contorno: geometría compleja gratis
Después de cada paso de transmisión, a los nodos adyacentes a las paredes sólidas les faltan algunas poblaciones entrantes, ya que habrían llegado desde el interior del sólido, que no tiene fluido para transmitir. La forma en que se llenan esas poblaciones faltantes define la condición de contorno.
Asas de rebote para paredes antideslizantes. Cualquier población que fluya hacia un nodo sólido se refleja en la dirección opuesta: f_ī = f*_i. La imagen física es una partícula que invierte su velocidad en la pared, y la condición de no deslizamiento (velocidad cero en el límite) surge automáticamente, sin necesidad de una restricción explícita de velocidad. La pared física se encuentra a medio camino entre el último nodo fluido y el primer nodo sólido, lo que le da al esquema una precisión espacial de segundo orden.
La elegante consecuencia es que cualquier forma de obstáculo es sólo una bandera booleana. Para agregar un cilindro, marque los nodos relevantes como sólidos. La regla de rebote maneja la física de límites en cada nodo marcado automáticamente, sin malla, sin correcciones de plantilla y sin cambios de código entre geometrías.
En la entrada utilizamos el esquema Zou-He: se impone una velocidad fija uw, y las tres poblaciones faltantes que apuntan hacia la derecha (f_1, f_5, f_8) se reconstruyen consistentemente a partir de las poblaciones conocidas restantes y la velocidad impuesta, recuperando la densidad de entrada correcta sin imponerla de forma independiente.
En la salida, una condición de gradiente cero copia todas las poblaciones de la penúltima columna, lo que permite que el fluido abandone el dominio libremente sin reflejar ondas espurias aguas arriba.
Resultados: dos regímenes
El dominio de simulación es una red de 400 × 400 con un obstáculo circular de 40 unidades de red de radio, centrado en un cuarto de la longitud del dominio y ligeramente desplazado verticalmente para romper la simetría. La velocidad de entrada es uw = 0,1 unidades de red, muy por debajo de cs ≈ 0,577 (supuesto de bajo Mach satisfecho en todo momento). Dos valores de τ revelan regímenes de flujo cualitativamente diferentes.
En τ = 0,75 (Re ≈ 96) la estela es constante y simétrica. La función de distribución se ha relajado cerca del equilibrio local en todas partes, ya que las colisiones son lo suficientemente vigorosas como para amortiguar cualquier perturbación antes de que pueda crecer. El flujo es viscoso y laminar.
En τ = 0,55 (Re ≈ 480) la imagen cambia completamente.
Los vórtices brotan alternativamente desde las superficies superior e inferior, formando la calle Kármán. Desde la perspectiva de la mecánica estadística, esta es la idea clave: la formación del vórtice es impulsada por la desviación de f de su forma de equilibrio de Maxwell-Boltzmann. Cuanto más pequeño es τ, las colisiones relajan f hacia el equilibrio más lentamente, lo que permite que el componente de no equilibrio f − f^eq crezca lo suficiente como para retroalimentar el campo de impulso macroscópico y sostener la oscilación. Un fluido en perfecto equilibrio no mostraría tal estructura.
El enfoque de bandera booleana significa que cambiar la geometría no cuesta nada. Aquí está el mismo solucionador, sin modificar, aplicado a tres obstáculos cuadrados seguidos:
Las estelas se acoplan, las calles de vórtice interactúan y emerge una rica estructura de flujo de múltiples cuerpos. El único cambio en el código fue marcar diferentes nodos como sólidos.
Escalando a MareNostrum 5
Con #pragma omp paralelo tanto en el bucle de colisión como en el de transmisión, el solucionador escala a múltiples núcleos con un mínimo esfuerzo. La elección importante es el patrón de transmisión: la transmisión pull garantiza que cada subproceso sea propietario exclusivo de su memoria de salida, eliminando el intercambio falso que hace que la implementación ingenua de inserción se estabilice temprano.
Realizamos una evaluación comparativa de un nodo completo de MareNostrum 5 en el Centro de Supercomputación de Barcelona: 2× Intel Xeon Platinum 8480+ (Sapphire Rapids), 112 núcleos en dos dominios NUMA.
Tres regímenes emergen claramente. De 1 a 32 subprocesos, la eficiencia se mantiene por encima del 75%, alcanzando una aceleración de 24,6x, lo que es una escala realmente fuerte para un código de plantilla con mucha memoria. Con 64 subprocesos, la eficiencia cae al 42 % a medida que el controlador de memoria se satura (los ~30 MB de datos de la función de distribución se leen y escriben más rápido de lo que el bus puede soportar). Con 112 subprocesos, la aceleración en realidad cae a 23,5 veces, porque abarcar ambos sockets físicos significa que la mitad de todos los accesos a la memoria cruzan la interconexión entre sockets, lo que agrega latencia al ancho de banda ya saturado. El punto óptimo práctico para una cuadrícula de 400 × 400 es 32 núcleos en un solo enchufe.
Un número en el que vale la pena detenerse: la línea de base de un solo subproceso después de la reescritura de la transmisión pull-stream es de 31,7 segundos, frente a los 113 segundos con la implementación push original. Esa mejora de 3,6 veces provino únicamente de un cambio de código, en un solo núcleo, antes de tocar cualquier paralelismo. El patrón de acceso a la memoria importa mucho antes que el hardware.
Cuándo usar LBM (y cuándo no)
LBM es la herramienta adecuada cuando la geometría es compleja o cambia con el tiempo, cuando necesita ejecutarse en hardware masivamente paralelo o cuando está modelando fenómenos de mesoescala como el deslizamiento de velocidad cerca de paredes para los cuales los solucionadores continuos necesitan correcciones empíricas. También es una base excelente para simulaciones multifásicas, el tema del próximo artículo de esta serie.
No es la herramienta adecuada para flujos de alto Mach (la expansión de Taylor de bajo orden se descompone por encima de Mach ~0,3), números de Reynolds muy altos (el requisito de estabilidad τ > 0,5 limita cómo puede funcionar la baja viscosidad sin operadores de colisión más avanzados) o situaciones en las que la calibración precisa de la unidad física es fundamental desde el principio.
Qué sigue: flujos multifásicos con Shan-Chen
El solucionador que se describe aquí maneja un fluido monofásico de un solo componente. Muchos de los problemas físicamente más interesantes involucran dos fluidos inmiscibles (aceite y agua, burbujas ascendentes, dinámica de las gotas en el flujo de aire).
El modelo Shan-Chen extiende LBM para manejar estos casos agregando fuerzas entre partículas de corto alcance que producen una ecuación de estado efectiva con una transición de fase. La tensión interfacial y la separación de fases surgen de las estadísticas de colisión en lugar de imponerse como condiciones límite. En el próximo artículo, derivaré el modelo Shan-Chen, lo implementaré como una modesta extensión de este solucionador y mostraré cómo se ve la separación de fases espontánea cuando se deja que dos fluidos compitan por la misma red.
Pensamiento final
Lo que encuentro más satisfactorio del LBM es lo que revela sobre la relación entre escalas. Las ecuaciones de Navier-Stokes parecen fundamentales: son ecuaciones de dinámica de fluidos, escritas en la escala que podemos ver y tocar. Pero no son fundamentales en absoluto. Son una consecuencia estadística de la colisión de muchas partículas, recuperables a partir de una descripción cinética en el límite apropiado.
LBM hace que esa relación sea concreta y computacional. La viscosidad no es un parámetro que usted establezca; es una escala de tiempo de colisión. La calle del vórtice de Kármán no es una inestabilidad fluida que programas; es lo que surge cuando las colisiones no pueden seguir el ritmo de la inercia. Cada fenómeno macroscópico en la simulación tiene una interpretación microscópica directa en la mecánica estadística subyacente.
Ese puente entre escalas, hecho tangible en ~200 líneas de C++, es el punto.
Referencias:
Para aquellos que quieran profundizar más, las siguientes son las referencias canónicas en orden creciente de dificultad técnica.
Mohamed, AA (2019). Método Lattice Boltzmann: fundamentos y aplicaciones de ingeniería con códigos informáticos. Saltador. El punto de entrada más accesible; Ejemplos trabajados y código en todas partes.
Krüger, T. y col. (2017). El método Lattice Boltzmann: principios y práctica. Saltador. La referencia integral; Se cubren derivaciones, condiciones de contorno, modelos multifase e implementación de GPU.
Bhatnagar, PL, Gross, EP y Krook, M. (1954). Un modelo para procesos de colisión en gases. Revisión física, 94(3), 511–525. El documento BGK original.
Chapman, S. y Cowling, TG (1970). La teoría matemática de los gases no uniformes. Prensa de la Universidad de Cambridge. Por la ampliación Chapman-Enskog con todo su rigor.
El código fuente completo de C++, los scripts de referencia de OpenMP y la implementación de referencia de Python/NumPy para este artículo están disponibles en GitHub.