Maxwell–Boltzmann en nueve flechas — por qué el equilibrio del LBM es un polinomio

En 1986, Frisch, Hasslacher y Pomeau colocaron partículas sobre una retícula y las hicieron imitar a Navier–Stokes. Aquello era el LGCA (Lattice Gas Cellular Automaton). Murió antes de cumplir cinco años — la condición entera de cero o una partícula por dirección generaba ruido imposible de borrar. En 1992, McNamara y Zanetti cambiaron enteros por reales. Llenaron la retícula con funciones de distribución $f_i$ en lugar de bits. El ruido desapareció, pero llegó otra pregunta. ¿Cómo pueden nueve números reales sustituir a la distribución infinito-dimensional de Maxwell–Boltzmann? Este artículo sigue la respuesta — un truncamiento Hermite de segundo orden hace de equilibrio en el LBM, y el mismo truncamiento lo encadena al régimen de bajo Mach.

De enteros a reales — el LGCA murió en 1986 y nació el LBM#

La idea del LGCA era elegante. Retícula cuadrada, seis direcciones de partícula, una sola regla de colisión. Promediada a escala macroscópica devolvía Navier–Stokes. Pero las partículas booleanas (cero o uno) amplificaban el ruido estadístico. Para matarlo hacían falta promedios temporales largos, y entonces el transitorio que se buscaba ver también desaparecía.

McNamara y Zanetti pusieron distribuciones reales $f_i(\mathbf{x}, t)$ en cada nodo. $f_i$ es el recuento esperado de partículas que en el instante $t$, en la posición $\mathbf{x}$, se mueven a lo largo de la dirección $\mathbf{c}_i$. El ruido se evapora por construcción. Cada paso se parte en dos — colisión y streaming.

$\tau$ es el tiempo de relajación y $f_i^{eq}$ es la distribución de equilibrio. La colisión pela la parte de no-equilibrio a escala $\tau$, y el streaming desplaza la distribución pelada una celda. Toda la cuestión profunda se concentra en $f_i^{eq}$ — ¿qué hace que basten nueve números reales?

Maxwell–Boltzmann carga una bomba de tiempo#

En la teoría cinética continua, el equilibrio es Maxwell–Boltzmann.

$\mathbf{c}$ es la velocidad molecular, $\mathbf{u}$ la velocidad macroscópica, $D$ la dimensión espacial, $T$ la temperatura, $R$ la constante de los gases. La bomba es que $\mathbf{c}$ es continua. Reducirla a nueve o diecinueve direcciones discretas exige cribar la distribución con cuidado.

Intento ingenuo — aproximar la integral gaussiana por una suma de Riemann en malla fija. Falla. Los momentos que necesita Navier–Stokes — $\rho$, $\rho \mathbf{u}$, el tensor de cantidad de movimiento $\rho \mathbf{u}\mathbf{u}$, la energía — no coinciden con la suma reticular. En cuanto se rompe la conservación, el LBM queda inservible.

La respuesta correcta es expansión funcional. Desarrollar Maxwell–Boltzmann en una serie de polinomios ortogonales, conservar los órdenes que llevan los momentos necesarios y dejar caer el resto. Esos polinomios son los de Hermite.

Los polinomios de Hermite cortan la distribución a tajadas#

Los polinomios probabilísticos de Hermite $He_n(c)$ son ortogonales bajo el peso gaussiano estándar $w(c) = (2\pi)^{-1/2} e^{-c^2/2}$.

Los primeros son $He_0 = 1$, $He_1 = c$, $He_2 = c^2 - 1$, $He_3 = c^3 - 3c$. Si se asume isotermia en 1D ($RT = 1$ tras adimensionalizar), la función generadora exponencial da de inmediato la identidad siguiente.

Cada término atrapa una rebanada del mundo. $n=0$ es el gaussiano estacionario; $n=1$ es una pequeña deriva; $n=2$ es la corrección de energía cinética media. Cuanto mayor el $n$, más asimétricos son los efectos de cola. Todo el juego del LBM consiste en cortar esta suma en un orden finito.

Toca la simulación de arriba. Sube la deriva $u$ a 0.4 y elige el orden de truncamiento $N=2$: los tres primeros momentos (densidad, cantidad de movimiento y energía) coinciden con Maxwell–Boltzmann al detalle. Empuja $u$ por encima de 1.0 y la curva discontinua se separa. Allí está, en una sola figura, el origen de la restricción de Mach del LBM, $|\mathbf{u}|/c_s \lesssim 0.3$.

Qué compra el truncamiento de segundo orden#

¿Qué momentos hacen falta para Navier–Stokes con esfuerzo viscoso de primer orden? Justo tres órdenes — $\rho$, $\rho \mathbf{u}$ y el flujo de cantidad de movimiento $\Pi_{\alpha\beta} = \rho u_\alpha u_\beta + p \delta_{\alpha\beta}$. Como $\rho \mathbf{u}\mathbf{u}$ es de orden $u^2$, la distribución debe acarrear un término en $u^2$.

Cortemos entonces la serie en $N=2$.

Esta es la forma continua del equilibrio LBM. $N=3$ y mayores sirven para momentos no viscosos más altos (flujo de calor), pero son redundantes en el LBM isotérmico estándar. Lo descartado explica los dos límites conocidos — la hipótesis isotérmica y el techo de Mach 0.3.

Cómo se construyó D2Q9#

Cambiar la integral continua $\int f^{eq}(c) \phi(c)\, dc$ por la suma reticular $\sum_i w_i f_i^{eq} \phi(\mathbf{c}_i)$ exige una cuadratura. Hecho clave — una cuadratura exacta hasta orden $N$ basta con que integre exactamente polinomios de ese orden.

En 2D, ser exactos hasta el término $u^2$ obliga a integrar exactamente polinomios de orden cinco bajo el peso gaussiano (el tensor de cantidad de movimiento produce un término $c^2 \times u^2$). La cuadratura Gauss–Hermite de tres puntos hace justamente ese trabajo.

En 1D, tres nodos $c \in \{-\sqrt{3}, 0, +\sqrt{3}\}$ con pesos $\{1/6, 2/3, 1/6\}$. En 2D el producto tensorial da $3 \times 3 = 9$ nodos — eso es D2Q9. Reescalando a unidades de retícula, $c_s = 1/\sqrt{3}$ y los vectores de velocidad reticular caen en coordenadas enteras $(0,0), (\pm 1, 0), (0, \pm 1), (\pm 1, \pm 1)$. Los pesos reticulares $w_i$ absorben el factor gaussiano en los pesos Gauss–Hermite.

Sustituyendo $c_s^2 = 1/3$ en la forma continua, queda el equilibrio sobre la retícula.

Esta línea es el corazón de todo código LBM isotérmico. Un polinomio modesto, pero encierra los tres primeros momentos de Maxwell–Boltzmann.

Python — un paso de colisión BGK#

Un paso de colisión sobre un solo nodo.

import numpy as np
 
cx = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
cy = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
 
def equilibrium(rho, ux, uy):
    """Equilibrio D2Q9 — nueve polinomios"""
    ciu = cx * ux + cy * uy            # producto escalar
    usq = ux * ux + uy * uy
    return w * rho * (1 + 3 * ciu + 4.5 * ciu**2 - 1.5 * usq)
 
def bgk_collision(f, tau):
    rho = f.sum()
    ux  = (cx * f).sum() / rho
    uy  = (cy * f).sum() / rho
    feq = equilibrium(rho, ux, uy)
    return f - (f - feq) / tau, rho, ux, uy
 
# distribución inicial fuera de equilibrio
f = w * 1.0
f[1] += 0.05; f[3] -= 0.04          # ligero sesgo hacia el este
 
for step in range(20):
    f, rho, ux, uy = bgk_collision(f, tau=0.8)
    err = np.linalg.norm(f - equilibrium(rho, ux, uy))
    print(f"step {step:2d}: rho={rho:.4f}, u=({ux:+.4f},{uy:+.4f}), ||delta||={err:.2e}")

En unos veinte pasos el residuo cae por debajo de $10^{-6}$ — los nueve $f_i$ se asientan sobre el polinomio de equilibrio. ¿Qué sobrevive en cada paso? La densidad $\rho$ y la cantidad de movimiento $\rho \mathbf{u}$. La colisión solo recorta los momentos de no-equilibrio.

Arrastra el deslizador de $\tau$ hacia 0.55 y el residuo se anula en uno o dos pasos. La viscosidad es $\nu = c_s^2 (\tau - 1/2)$, así que cuando $\tau \to 1/2$ la viscosidad desaparece y la retícula se vuelve inestable. Empujar $\tau$ hacia arriba aumenta la viscosidad pero ralentiza la relajación. Esta es la perilla que toda persona que trabaja con LBM termina ajustando.

Tres líneas para recordar#

• El equilibrio $f_i^{eq}$ del LBM es Maxwell–Boltzmann desarrollado en polinomios de Hermite y cortado en segundo orden — ese corte conserva exactamente los momentos $\rho$, $\rho \mathbf{u}$ y el flujo de cantidad de movimiento.
• D2Q9 se deduce como el producto tensorial 2D de la cuadratura Gauss–Hermite de tres puntos, la regla mínima que integra exactamente el equilibrio de segundo orden. El 9 no es arbitrario; es un subproducto de la integración ortogonal.
• El orden de truncamiento fija todos los límites. La hipótesis isotérmica y el techo $|\mathbf{u}|/c_s \lesssim 0.3$ se siguen automáticamente de una sola frase: "cortamos la serie en $N=2$".