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

El gas de red murió en 1986. La ocupación era cero o uno, y el ruido no se iba con el promedio. En 1992 esos enteros pasaron a una distribución real $f_i$. El ruido desapareció. La pregunta se quedó. Cómo nueve números reales sustituyen a un Maxwell–Boltzmann de dimensión infinita.

Al mover la deriva $u$ y el orden de truncamiento $N$, la curva sólida de Maxwell–Boltzmann se apoya en el polinomio discontinuo o se separa de él.

Con $N=2$ y $u=0.4$, densidad, cantidad de movimiento y energía quedan superpuestas. Pasado $u=1$, solo la curva discontinua se despega. El techo de Mach $|\mathbf{u}|/c_s \lesssim 0.3$ no es ese despegue. Los números se parten más adelante.

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

El LGCA (Lattice Gas Cellular Automaton) era una malla cuadrada, seis direcciones, una regla de colisión. El promedio macroscópico volvía como Navier–Stokes. Las partículas booleanas (cero o uno) agrandaban el ruido estadístico. Alargar el promedio temporal mata el ruido y, con él, el transitorio que se quería ver.

McNamara y Zanetti colocaron una distribución real $f_i(\mathbf{x}, t)$ en cada nodo. $f_i$ es el recuento esperado de partículas en la posición $\mathbf{x}$ y el instante $t$ que se mueven a lo largo de $\mathbf{c}_i$. Un paso se parte en dos. Colisión, luego streaming.

$\tau$ es el tiempo de relajación (relaxation time). $f_i^{eq}$ es la distribución de equilibrio. La colisión recorta la parte de no equilibrio a escala $\tau$, y el streaming desplaza esa distribución una celda. Nueve números alcanzan solo por la forma de $f_i^{eq}$.

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, $T$ la temperatura, $R$ la constante de los gases. La bomba es que $\mathbf{c}$ es continua. Dejar solo nueve o diecinueve direcciones exige adelgazar esta distribución.

Sustituir la integral gaussiana por una suma de rectángulos falla. Los momentos que pide 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. Rota la conservación y el LBM se detiene ahí.

El camino que queda es un desarrollo. Se expande Maxwell–Boltzmann en polinomios ortogonales y se conservan solo los momentos necesarios. Esos polinomios son los de Hermite.

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

Los polinomios probabilísticos de Hermite (probabilist Hermite polynomial) $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$. Con el equilibrio unidimensional en isotermo, es decir $RT=1$, la función generatriz da lo siguiente.

$n=0$ es una gaussiana en reposo. $n=1$ es una deriva pequeña. $n=2$ corrige la energía cinética media. Los $n$ más altos son colas asimétricas. El LBM corta la suma en un orden finito.

En ese desarrollo isotermo $RT=1$, así que la velocidad térmica vale 1. La $u$ de la figura ya está dividida por esa velocidad. El techo escrito como 0.3 es $|u|$ en esta escala. $u=0.4$ queda 0.1 por encima del techo, y los tres primeros momentos con $N=2$ siguen coincidiendo. El despegue en $|u|>1$ ya pasó Mach 1. El 0.3 y el 1 no son el mismo fallo.

El término $n=3$ es $(u^3/6)\,He_3(c)$, con $He_3=c^3-3c$. El coeficiente solo vale $0.3^3/6=0.027/6=0.0045$ en $u=0.3$. En $u=0.4$ vale $0.4^3/6=0.064/6\approx 0.0107$, cerca de 2.37 veces el valor del techo. En $u=1$ vale $1/6\approx 0.1667$, cerca de 37 veces 0.0045. Esa razón es $1^3/0.3^3=1/0.027$. $N=2$ tira el término entero. Así, densidad, cantidad de movimiento y energía pueden coincidir en $u=0.4$ mientras el cúbico descartado ya supera el doble de su tamaño en el techo. El despegue en $u=1$ es donde ese coeficiente ha subido hasta cerca de 0.17. El techo de Mach en 0.3 no es el punto en que la curva se separa a simple vista. Es una línea trazada de antemano, en el instante en que el coeficiente cúbico pasa de 0.0045.

Qué compra el truncamiento de segundo orden#

Navier–Stokes viscoso de primer orden necesita tres capas. $\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 $u^2$, la distribución tiene que cargar un término en $u^2$.

La serie se corta en $N=2$.

Esa es la forma continua del equilibrio. $N=3$ y más alimentan el flujo de calor. El LBM isotermo estándar no los conserva. Los términos descartados dejan a la vez la hipótesis isoterma y el techo de Mach 0.3. Al devolver el cúbico, la hipótesis isoterma se rompe primero. El flujo de calor es el trozo que carga ese término, así que un código detenido en segundo orden no puede tomar un problema donde la temperatura fluye y seguir con este polinomio. El techo de Mach y la hipótesis isoterma son dos caras de un mismo corte. Aunque la velocidad se mantenga bajo 0.173, una pared caliente y una pared fría en la misma malla no reciben transporte de calor de este equilibrio.

Cómo se construyó D2Q9#

Pasar la integral $\int f^{eq}(c)\,\phi(c)\,dc$ a la suma reticular $\sum_i w_i f_i^{eq}\phi(\mathbf{c}_i)$ pide una cuadratura. Igualar hasta orden $N$ significa integrar con exactitud los polinomios de ese orden.

En dos dimensiones, la exactitud hasta $u^2$ exige polinomios de quinto orden exactos bajo el peso gaussiano. El tensor de cantidad de movimiento produce un término $c^2\times u^2$. La cuadratura de Gauss–Hermite de tres puntos hace esa integral.

Los nodos unidimensionales son $c\in\{-\sqrt{3},0,+\sqrt{3}\}$, con pesos $\{1/6,2/3,1/6\}$. Por qué esas tres raíces caen ahí lo cuenta la cuadratura de Hermite que plegó un momento de sexto orden sobre nueve velocidades. El producto tensorial bidimensional es $3\times 3=9$. Eso es D2Q9. En unidades de malla $c_s=1/\sqrt{3}$, y las velocidades caen en coordenadas enteras $(0,0),(\pm 1,0),(0,\pm 1),(\pm 1,\pm 1)$. Los pesos $w_i$ son los de Gauss–Hermite con el factor gaussiano absorbido.

El equilibrio sobre la malla es la forma continua con $c_s^2=1/3$ sustituido.

$1/c_s^2=3$, $1/(2c_s^4)=9/2$, $1/(2c_s^2)=3/2$. El 3, el $9/2$ y el $3/2$ del polinomio son esos tres números. Con $\mathbf{u}=0$ el corchete vale 1 y $f_i^{eq}=w_i\rho$. Los pesos suman $4/9+4\cdot(1/9)+4\cdot(1/36)=1$. Cuatro novenos de una masa en reposo quedan en la dirección 0. Los tres pesos de la tabla son productos del conjunto unidimensional $\{1/6,\ 2/3,\ 1/6\}$. El reposo es $(2/3)\times(2/3)=4/9$. Una dirección de eje es $(2/3)\times(1/6)=1/9$. Una diagonal es $(1/6)\times(1/6)=1/36$. Cuatro ejes y cuatro diagonales cierran otra vez la suma, $4/9+4\cdot(1/9)+4\cdot(1/36)=1$. El nueve y los tres denominadores salen del mismo producto tensorial.

Esa línea es el equilibrio de un código LBM isotermo. Un polinomio, con los tres primeros momentos de Maxwell–Boltzmann encerrados dentro.

Python — un paso de colisión BGK#

Un paso de colisión, un 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}")

El estado inicial parte de $f=w$, suma 0.05 a $f_1$ y resta 0.04 de $f_3$. La masa es $1+0.05-0.04=1.01$. La cantidad de movimiento en $x$ es $1\cdot 0.05+(-1)\cdot(-0.04)=0.09$. Así $u_x=0.09/1.01\approx 0.0891$ y $u_y=0$. Al dividir por $c_s=1/\sqrt{3}$ el Mach queda cerca de 0.154. La velocidad de malla del techo 0.3 es $0.3/\sqrt{3}\approx 0.173$, y 0.0891 queda por debajo. El residuo que la corrida lleva bajo $10^{-6}$ en unos veinte pasos es el acercamiento al equilibrio sobre el par conservado $(1.01,\ 0.0891)$. Cada colisión conserva solo $\rho$ y $\rho\mathbf{u}$.

Al bajar $\tau$ hacia 0.55 el residuo desaparece en uno o dos pasos. La viscosidad es $\nu=c_s^2(\tau-1/2)$. Con el $\tau=0.8$ del código, $\nu=(1/3)\cdot 0.3=0.1$. Con $\tau=0.55$, $\nu=(1/3)\cdot 0.05=1/60\approx 0.0167$, un sexto del valor del código. Si $\tau$ llega a $1/2$, la viscosidad se anula y la malla se vuelve inestable. Al subir $\tau$, la viscosidad crece y el arrastre hacia el equilibrio se vuelve lento. Un paso multiplica la parte de no equilibrio por $1-1/\tau$. Con el $\tau=0.8$ del código, $1/\tau=1.25$, así que el factor es $-0.25$. El signo se invierte. Eso es sobrerrelajación. Veinte factores de $0.25$ dan $0.25^{20}\approx 9.1\times 10^{-13}$. Una desviación de orden 1 ya queda muy por debajo de $10^{-6}$ a los veinte pasos. Tomar veinte pasos como marca de $10^{-6}$ es holgado con este factor solo. Con $\tau=1$ el factor es 0 y el estado cae en el equilibrio en un paso. El 0.8 del código está un puesto al lado de ese equilibrio de un solo paso.

Pasado Mach 0.3, este polinomio no se sostiene#

$N=2$ ajusta $\rho$, $\rho\mathbf{u}$ y el flujo de cantidad de movimiento. Lo que está por encima se descartó. Cuando la velocidad de malla pasa de unos 0.173, los términos descartados dejan de ser pequeños. Bajar $\tau$ para que la viscosidad caiga de $0.1$ a $1/60$ no cambia el polinomio de equilibrio. Cruzar el techo con la velocidad y bajar solo la viscosidad son fallos distintos.

Si una viscosidad ajustada en un marco en reposo sobrevive a un desplazamiento de ese marco no tiene respuesta en esta expresión. El término cúbico descartado guarda esa respuesta. Si el tercer momento ausente cambia la viscosidad al mover el marco sigue fuera de este polinomio.