Métodos de Boltzmann reticular

Los métodos de Boltzmann reticular, también llamado métodos de Boltzmann en red, y conocido por LBM por las siglas del título en inglés, originados a partir del método de los autómatas de gas en red (LGA) (modelos de Hardy-Pomeau-Pazzis y Uriel Frisch-Hasslacher-Pomeau), son una clase de métodos de dinámica de fluidos computacional (CFD) para la simulación de fluidos. En lugar de resolver directamente las ecuaciones de Navier-Stokes, se simula la densidad de un fluido en una red mediante procesos de flujo y colisión (relajación),^[1]​ El método es versátil^[1]​ ya que el fluido modelo puede imitar fácilmente el comportamiento habitual de los fluidos, como la coexistencia de vapor y líquido, lo que permite simular sistemas fluidos como las gotas de líquido. Además, se pueden simular fácilmente fluidos en entornos complejos, como medios porosos, mientras que con otros métodos CFD puede resultar difícil trabajar con límites complejos.

A diferencia de los métodos CFD, que resuelven numéricamente las ecuaciones de conservación de las propiedades macroscópicas (es decir, masa, momento y energía), el LBM modela el fluido como partículas ficticias, y dichas partículas realizan procesos consecutivos de propagación y colisión sobre una red discreta. Debido a su naturaleza particulada y a su dinámica local, el LBM presenta varias ventajas con respecto a otros métodos CFD convencionales, especialmente en el tratamiento de límites complejos, la incorporación de interacciones microscópicas y la paralelización del algoritmo. ^[2]​ Una interpretación diferente de la ecuación de Boltzmann en red es la de una ecuación de Boltzmann de velocidad discreta. Los métodos numéricos de solución del sistema de ecuaciones diferenciales parciales dan lugar a un mapa discreto, que puede interpretarse como la propagación y colisión de partículas ficticias.

En un algoritmo, hay pasos de colisión y streaming. Estos hacen evolucionar la densidad del fluido ${\displaystyle \rho ({\vec {x}},t)}$, siendo ${\displaystyle {\vec {x}}}$ la posición y ${\displaystyle t}$ el tiempo. Como el fluido se encuentra en una red, la densidad tiene varios componentes ${\displaystyle f_{i},i=0,\ldots ,a}$ iguales al número de vectores de red conectados a cada punto de la red. A modo de ejemplo, se muestran aquí los vectores de red para una red simple utilizada en simulaciones en dos dimensiones. Esta red se suele denominar D2Q9, por sus dos dimensiones y nueve vectores: cuatro vectores a lo largo del norte, este, sur y oeste, más cuatro vectores a las esquinas de un cuadrado unitario, más un vector con ambos componentes iguales a cero. Entonces, por ejemplo, el vector ${\displaystyle {\vec {e}}_{4}=(0,-1)}$, es decir, apunta hacia el sur y, por lo tanto, no tiene componente ${\displaystyle x}$, pero tiene una componente ${\displaystyle y}$ de ${\displaystyle -1}$. Así, una de las nueve componentes de la densidad total en el punto central de la red, ${\displaystyle f_{4}({\vec {x}},t)}$, es la parte del fluido en el punto ${\displaystyle {\vec {x}}}$ que se mueve hacia el sur, a una velocidad en unidades de red de uno.

A continuación, los pasos que hacen evolucionar el fluido en el tiempo son:^[1]​

El paso de colisión

${\displaystyle f_{i}^{\ast }({\vec {x}},t)=f_{i}({\vec {x}},t)+{\frac {f_{i}^{eq}({\vec {x}},t)-f_{i}({\vec {x}},t)}{\tau _{f}}}\,\!}$

que es el Bhatnagar Gross y Krook (BGK),^[3]​ modelo de relajación hasta el equilibrio mediante colisiones entre las moléculas de un fluido. ${\displaystyle f_{i}^{eq}({\vec {x}},t)}$ es la densidad de equilibrio a lo largo de la dirección «i» en la densidad actual en ese punto, lo que puede expresarse mediante una aproximación de Taylor (véase más abajo, en Ecuaciones matemáticas para simulaciones):

${\displaystyle f_{i}^{eq}=\omega _{i}\rho \left(1+{\frac {3{\vec {e}}_{i}{\vec {u}}}{c^{2}}}+{\frac {9({\vec {e}}_{i}{\vec {u}})^{2}}{2c^{4}}}-{\frac {3({\vec {u}})^{2}}{2c^{2}}}\right)}$

El modelo asume que el fluido se relaja localmente hasta alcanzar el equilibrio en un intervalo de tiempo característico ${\displaystyle \tau _{f}}$. Este intervalo de tiempo determina la viscosidad cinemática: cuanto mayor es el intervalo de tiempo, mayor es la viscosidad cinemática.

El paso de streaming

${\displaystyle f_{i}({\vec {x}}+{\vec {e}}_{i},t+\delta _{t})=f_{i}^{\ast }({\vec {x}},t)\,\!}$

Como ${\displaystyle f_{i}^{\ast }({\vec {x}},t)}$ es, por definición, la densidad del fluido en el punto ${\displaystyle {\vec {x}}}$ en el momento ${\displaystyle t}$ que se mueve a una velocidad de ${\displaystyle {\vec {e}}_{i}}$ por periodo de tiempo, entonces en el siguiente periodo de tiempo ${\displaystyle t+\delta _{t}}$ habrá fluido hasta el punto ${\displaystyle {\vec {x}}+{\vec {e}}_{i}}$.

• El LBM se diseñó desde cero para funcionar de manera eficiente en arquitecturas masivamente paralelas, que van desde FPGA y DSP integrados y económicos hasta GPU y clústeres heterogéneos y superordenadores (incluso con una red de interconexión lenta). Permite realizar cálculos físicos complejos y algoritmos sofisticados. La eficiencia conduce a un nuevo nivel cualitativo de comprensión, ya que permite resolver problemas que antes no se podían abordar (o solo con una precisión insuficiente).
• El método tiene su origen en una descripción molecular de un fluido y puede incorporar directamente términos físicos derivados del conocimiento de la interacción entre moléculas. Por lo tanto, es un instrumento indispensable en la investigación fundamental, ya que acorta el ciclo entre la elaboración de una teoría y la formulación de un modelo numérico correspondiente.
• Preprocesamiento automatizado de datos y generación de redes en un tiempo que representa una pequeña fracción del total de la simulación.
• Análisis, posprocesamiento y evaluación de datos en paralelo.
• Flujo multifásico totalmente resuelto con pequeñas gotas y burbujas.
• Flujo totalmente resuelto a través de geometrías complejas y medios porosos.
• Flujo complejo y acoplado con transferencia de calor y reacciones químicas.

Al igual que con el CFD basado en Navier-Stokes, los métodos LBM se han combinado con éxito con soluciones térmicas específicas para permitir la simulación de la transferencia de calor (conducción, convección y radiación basadas en sólidos). En los modelos multifásicos/multicomponentes, el espesor de la interfaz suele ser grande y la relación de densidad a través de la interfaz es pequeña en comparación con los fluidos reales. Recientemente, este problema ha sido resuelto por Yuan y Schaefer, quienes mejoraron los modelos de Shan y Chen, Swift, y He, Chen y Zhang. Consiguieron alcanzar relaciones de densidad de 1000:1 simplemente cambiando la ecuación de estado. Se ha propuesto aplicar la transformación galileana para superar la limitación de la modelización de flujos de fluidos a alta velocidad. ^[4]​ Los rápidos avances de este método también han permitido simular con éxito la microfluídica,^[5]​ Sin embargo, por ahora, el LBM sigue teniendo limitaciones a la hora de simular flujos con un número de Knudsen elevado, en los que se utilizan en su lugar los métodos de Montecarlo, y los flujos con un número de Mach elevado en aerodinámica siguen siendo difíciles para el LBM, además de que no existe un esquema termohidrodinámico coherente.^[6]​

El LBM se originó a partir del método autómata de gas en celosía (LGA por las siglas en inglés), que puede considerarse un modelo simplificado de dinámica molecular ficticia en el que el espacio, el tiempo y las velocidades de las partículas son discretos. Por ejemplo, en el Modelo FHP bidimensional, cada nodo de la red está conectado a sus vecinos por seis velocidades reticulares en una red triangular; puede haber 0 o 1 partículas en un nodo de la red moviéndose con una velocidad reticular determinada. Tras un intervalo de tiempo, cada partícula se moverá hacia el nodo vecino en su dirección; este proceso se denomina paso de propagación o streaming. Cuando más de una partícula llega al mismo nodo desde diferentes direcciones, chocan y cambian sus velocidades según un conjunto de reglas de colisión. Los pasos de flujo y los pasos de colisión se alternan. Las reglas de colisión adecuadas deben conservar el número de partículas (masa), el momento y la energía antes y después de la colisión. El LGA adolece de varios defectos innatos para su uso en simulaciones hidrodinámicas: falta de invarianza galileana para flujos rápidos, ruido estadístico y mala escala del número de Reynolds con el tamaño de la red. Sin embargo, el LGA es muy adecuado para simplificar y ampliar el alcance de los modelos de difusión de reacciones y dinámica molecular.

La principal motivación para la transición de LGA a LBM fue el deseo de eliminar el ruido estadístico sustituyendo el número booleano de partículas en una dirección de la red por su media ensamblada, la denominada función de distribución de densidad. Junto con esta sustitución, la regla de colisión discreta también se sustituye por una función continua conocida como operador de colisión. En el desarrollo de LBM, una simplificación importante es aproximar el operador de colisión con el término de relajación Bhatnagar-Gross-Krook (BGK). Este modelo BGK en red (LBGK) hace que las simulaciones sean más eficientes y permite flexibilidad en los coeficientes de transporte. Por otro lado, se ha demostrado que el esquema LBM también puede considerarse una forma discretizada especial de la ecuación continua de Boltzmann. A partir de la teoría de Chapman-Enskog, se pueden recuperar las ecuaciones de continuidad y Navier-Stokes que rigen el algoritmo LBM.

Los modelos de Boltzmann en red pueden funcionar en diferentes redes, tanto cúbicas como triangulares, y con o sin partículas de reposo en la función de distribución discreta.

Una forma popular de clasificar los diferentes métodos según la red es el esquema DnQm. Aquí, «Dn» significa «n dimensiones», mientras que «Qm» significa «m velocidades». Por ejemplo, D3Q15 es un modelo de Boltzmann en red tridimensional sobre una cuadrícula cúbica, con partículas en reposo presentes. Cada nodo tiene forma de cristal y puede entregar partículas a 15 nodos: cada uno de los 6 nodos vecinos que comparten una superficie, los 8 nodos vecinos que comparten una esquina y él mismo.^[7]​ (El modelo D3Q15 no contiene partículas que se mueven a los 12 nodos vecinos que comparten un borde; si se añadieran, se crearía un modelo «D3Q27»).

Las cantidades reales, como el espacio y el tiempo, deben convertirse a unidades de red antes de la simulación. Las cantidades adimensionales, como el número de Reynolds, permanecen iguales.

En la mayoría de las simulaciones de Boltzmann en celosía, ${\displaystyle \delta _{x}\,\!}$ es la unidad básica para el espaciado de la celosía, por lo que si el dominio de longitud ${\displaystyle L\,\!}$ tiene ${\displaystyle N\,\!}$ unidades de celosía a lo largo de toda su longitud, la unidad de espacio se define simplemente como ${\displaystyle \delta _{x}=L/N\,\!}$. Las velocidades en las simulaciones de Boltzmann en red se suelen expresar en términos de la velocidad del sonido. Por lo tanto, la unidad de tiempo discreta se puede expresar como ${\displaystyle \delta _{t}={\frac {\delta _{x}}{C_{s}}}\,\!}$, donde el denominador ${\displaystyle C_{s}}$ es la velocidad física del sonido.^[8]​

Para flujos a pequeña escala (como los que se ven en la mecánica de medios porosos), operar con la velocidad real del sonido puede dar lugar a pasos de tiempo inaceptablemente cortos. Por lo tanto, es habitual elevar el número de Mach de la red a un valor mucho mayor que el número de Mach real y compensarlo aumentando también la viscosidad para preservar el número de Reynolds.^[9]​

La simulación de flujos multifásicos/multicomponentes siempre ha sido un reto para la CFD convencional debido a las interfaces móviles y deformables. Más fundamentalmente, las interfaces entre diferentes fases (líquido y vapor) o componentes (por ejemplo, aceite y agua) se originan a partir de las interacciones específicas entre las moléculas de los fluidos. Por lo tanto, es difícil implementar estas interacciones microscópicas en la ecuación macroscópica de Navier-Stokes. Sin embargo, en el LBM, la cinética de partículas proporciona una forma relativamente fácil y consistente de incorporar las interacciones microscópicas subyacentes mediante la modificación del operador de colisión. Se han desarrollado varios modelos LBM multifásicos/multicomponentes. Aquí, las separaciones de fases se generan automáticamente a partir de la dinámica de las partículas y no se necesita ningún tratamiento especial para manipular las interfaces, como en los métodos CFD tradicionales. Se pueden encontrar aplicaciones exitosas de los modelos LBM multifásicos/multicomponentes en diversos sistemas fluidos complejos, incluyendo la inestabilidad de interfaces, la dinámica de burbujas/gotas, la humectación en superficies sólidas, el deslizamiento interfacial y las deformaciones electrohidrodinámicas de gotas.

Recientemente se ha propuesto un modelo de Boltzmann en red para la simulación de la combustión de mezclas de gases capaz de adaptarse a variaciones significativas de densidad en regímenes de bajo número de Mach. ^[10]​

A este respecto, cabe señalar que, dado que el LBM se ocupa de un conjunto más amplio de campos (en comparación con el CFD convencional), la simulación de mezclas de gases reactivos presenta algunos retos adicionales en términos de demanda de memoria en lo que se refiere a mecanismos de combustión grandes y detallados. No obstante, estas cuestiones pueden abordarse recurriendo a técnicas sistemáticas de reducción de modelos.^[11]​^[12]​^[13]​

Actualmente (2009), el método térmico de Boltzmann en red (TLBM) se divide en tres categorías: el enfoque multivelocidad,^[14]​ el enfoque escalar pasivo,^[15]​ y la distribución de la energía térmica.^[16]​

Partiendo de la ecuación discreta de Boltzmann en red (también denominada ecuación LBGK debido al operador de colisión utilizado). Primero hacemos una expansión en serie de Taylor de segundo orden sobre el lado izquierdo de la LBE. Se elige esta expansión en lugar de una expansión de Taylor de primer orden más simple, ya que no se puede recuperar la LBE discreta. Al hacer la expansión en serie de Taylor de segundo orden, el término derivado cero y el primer término de la derecha se cancelan, dejando solo los términos de la primera y segunda derivada de la expansión de Taylor y el operador de colisión:

${\displaystyle f_{i}({\vec {x}}+{\vec {e}}_{i}\delta _{t},t+\delta _{t})=f_{i}({\vec {x}},t)+{\frac {\delta _{t}}{\tau _{f}}}(f_{i}^{eq}-f_{i}).}$

Para simplificar, escriba ${\displaystyle f_{i}({\vec {x}},t)}$ como ${\displaystyle f_{i}}$. La expansión en serie de Taylor ligeramente simplificada es entonces la siguiente, donde ":" es el producto colono entre díadas:

${\displaystyle {\frac {\partial f_{i}}{\partial t}}+{\vec {e}}_{i}\cdot \nabla f_{i}+\left({\frac {1}{2}}{\vec {e}}_{i}{\vec {e}}_{i}:\nabla \nabla f_{i}+{\vec {e}}_{i}\cdot \nabla {\frac {\partial f_{i}}{\partial t}}+{\frac {1}{2}}{\frac {\partial ^{2}f_{i}}{\partial t^{2}}}\right)={\frac {1}{\tau }}(f_{i}^{eq}-f_{i}).}$

Al expandir la función de distribución de partículas en componentes de equilibrio y de no equilibrio y utilizando la expansión de Chapman-Enskog, donde ${\displaystyle K}$ es el número de Knudsen, la LBE expandida de Taylor se puede descomponer en diferentes magnitudes de orden para el número de Knudsen con el fin de obtener las ecuaciones continuas adecuadas:

${\displaystyle f_{i}=f_{i}^{\text{eq}}+Kf_{i}^{\text{neq}},}$

${\displaystyle f_{i}^{\text{neq}}=f_{i}^{(1)}+Kf_{i}^{(2)}+O(K^{2}).}$

Las distribuciones en equilibrio y en desequilibrio satisfacen las siguientes relaciones con sus variables macroscópicas (estas se utilizarán más adelante, una vez que las distribuciones de partículas estén en la «forma correcta» para poder escalar desde el nivel de partículas al macroscópico):

${\displaystyle \rho =\sum _{i}f_{i}^{\text{eq}},}$

${\displaystyle \rho {\vec {u}}=\sum _{i}f_{i}^{\text{eq}}{\vec {e}}_{i},}$

${\displaystyle 0=\sum _{i}f_{i}^{(k)}\qquad {\text{for }}k=1,2,}$

${\displaystyle 0=\sum _{i}f_{i}^{(k)}{\vec {e}}_{i}.}$

La expansión de Chapman-Enskog es entonces:

${\displaystyle {\frac {\partial }{\partial t}}=K{\frac {\partial }{\partial t_{1}}}+K^{2}{\frac {\partial }{\partial t_{2}}}\qquad {\text{for }}t_{2}({\text{escala de tiempo difusiva}})\ll t_{1}({\text{escala de tiempo convectiva}}),}$

${\displaystyle {\frac {\partial }{\partial x}}=K{\frac {\partial }{\partial x_{1}}}.}$

Sustituyendo el equilibrio ampliado y el desequilibrio en la expansión de Taylor y separando en diferentes órdenes de ${\displaystyle K}$, se obtienen casi las ecuaciones continuas.

Para el orden ${\displaystyle K^{0}}$:

${\displaystyle {\frac {\partial f_{i}^{\text{eq}}}{\partial t_{1}}}+{\vec {e}}_{i}\nabla _{1}f_{i}^{\text{eq}}=-{\frac {f_{i}^{(1)}}{\tau }}.}$

For order ${\displaystyle K^{1}}$:

${\displaystyle {\frac {\partial f_{i}^{(1)}}{\partial t_{1}}}+{\frac {\partial f_{i}^{\text{eq}}}{\partial t_{2}}}+{\vec {e}}_{i}\nabla f_{i}^{(1)}+{\frac {1}{2}}{\vec {e}}_{i}{\vec {e}}_{i}:\nabla \nabla f_{i}^{\text{eq}}+{\vec {e}}_{i}\cdot \nabla {\frac {\partial f_{i}^{\text{eq}}}{\partial t_{1}}}+{\frac {1}{2}}{\frac {\partial ^{2}f_{i}^{\text{eq}}}{\partial t_{1}^{2}}}=-{\frac {f_{i}^{(2)}}{\tau }}.}$

A continuación, la segunda ecuación se puede simplificar con un poco de álgebra y la primera ecuación queda así:

${\displaystyle {\frac {\partial f_{i}^{\text{eq}}}{\partial t_{2}}}+\left(1-{\frac {1}{2\tau }}\right)\left[{\frac {\partial f_{i}^{(1)}}{\partial t_{1}}}+{\vec {e}}_{i}\nabla _{1}f_{i}^{(1)}\right]=-{\frac {f_{i}^{(2)}}{\tau }}.}$

Aplicando las relaciones entre las funciones de distribución de partículas y las propiedades macroscópicas anteriores, se obtienen las ecuaciones de masa y momento:

${\displaystyle {\frac {\partial \rho }{\partial t}}+\nabla \cdot \rho {\vec {u}}=0,}$

${\displaystyle {\frac {\partial \rho {\vec {u}}}{\partial t}}+\nabla \cdot \Pi =0.}$

El tensor de flujo de impulso ${\displaystyle \Pi }$ tiene entonces la siguiente forma:

${\displaystyle \Pi _{xy}=\sum _{i}{\vec {e}}_{ix}{\vec {e}}_{iy}\left[f_{i}^{eq}+\left(1-{\frac {1}{2\tau }}\right)f_{i}^{(1)}\right],}$

donde ${\displaystyle {\vec {e}}_{ix}{\vec {e}}_{iy}}$ es una abreviatura para el cuadrado de la suma de todos los componentes de ${\displaystyle {\vec {e}}_{i}}$ (es decir, ${\displaystyle \textstyle \left(\sum _{x}{\vec {e}}_{ix}\right)^{2}=\sum _{x}\sum _{y}{\vec {e}}_{ix}{\vec {e}}_{iy}}$), y la distribución de partículas en equilibrio con segundo orden comparable a la ecuación de Navier-Stokes es:

${\displaystyle f_{i}^{\text{eq}}=\omega _{i}\rho \left(1+{\frac {{\vec {e}}_{i}{\vec {u}}}{c_{s}^{2}}}+{\frac {({\vec {e}}_{i}{\vec {u}})^{2}}{2c_{s}^{4}}}-{\frac {{\vec {u}}^{2}}{2c_{s}^{2}}}\right).}$

La distribución de equilibrio solo es válida para velocidades pequeñas o números de Mach pequeños. Al insertar la distribución de equilibrio de nuevo en el tensor de flujo, se obtiene:

${\displaystyle \Pi _{xy}^{(0)}=\sum _{i}{\vec {e}}_{ix}{\vec {e}}_{iy}f_{i}^{eq}=p\delta _{xy}+\rho u_{x}u_{y},}$

${\displaystyle \Pi _{xy}^{(1)}=\left(1-{\frac {1}{2\tau }}\right)\sum _{i}{\vec {e}}_{ix}{\vec {e}}_{iy}f_{i}^{(1)}=\nu \left(\nabla _{x}\left(\rho {\vec {u}}_{y}\right)+\nabla _{y}\left(\rho {\vec {u}}_{x}\right)\right).}$

Finalmente, se recuperan las Ecuaciones de Navier-Stokes bajo el supuesto de que la variación de densidad es pequeña:

${\displaystyle \rho \left({\frac {\partial {\vec {u}}_{x}}{\partial t}}+\nabla _{y}\cdot {\vec {u}}_{x}{\vec {u}}_{y}\right)=-\nabla _{x}p+\nu \nabla _{y}\cdot \left(\nabla _{x}\left(\rho {\vec {u}}_{y}\right)+\nabla _{y}\left(\rho {\vec {u}}_{x}\right)\right).}$

Esta derivación sigue el trabajo de Chen y Doolen.^[17]​

La ecuación continua de Boltzmann es una ecuación de evolución para una función de distribución de probabilidad de una sola partícula ${\displaystyle f({\vec {x}},{\vec {e}}_{i},t)}$ y la función de distribución de densidad de energía interna ${\displaystyle g({\vec {x}},{\vec {e}}_{i},t)}$ (He et al.) son, respectivamente:

${\displaystyle \partial _{t}f+({\vec {e}}\cdot \nabla )f+F\partial _{v}f=\Omega (f),}$

${\displaystyle \partial _{t}g+({\vec {e}}\cdot \nabla )g+G\partial _{v}f=\Omega (g),}$

donde ${\displaystyle g({\vec {x}},{\vec {e}}_{i},t)}$ está relacionado con ${\displaystyle f({\vec {x}},{\vec {e}}_{i},t)}$ mediante

${\displaystyle g({\vec {x}},{\vec {e}}_{i},t)={\frac {({\vec {e}}-{\vec {u}})^{2}}{2}}f({\vec {x}},{\vec {e}}_{i},t),}$

${\displaystyle F}$ es una fuerza externa, ${\displaystyle \Omega }$ es una integral de colisión y ${\displaystyle {\vec {e}}}$ (también etiquetada como ${\displaystyle {\vec {\xi }}}$ en la literatura) es la velocidad microscópica. La fuerza externa ${\displaystyle F}$ está relacionada con la fuerza externa de temperatura ${\displaystyle G}$ mediante la relación siguiente. Una prueba típica para un modelo es la convección de Rayleigh-Bénard para ${\displaystyle G}$.

${\displaystyle F={\frac {{\vec {G}}\cdot ({\vec {e}}-{\vec {u}})}{RT}}f^{\text{eq}},}$

${\displaystyle {\vec {G}}=\beta g_{0}(T-T_{avg}){\vec {k}}.}$

Las variables macroscópicas, como la densidad ${\displaystyle \rho }$, la velocidad ${\displaystyle {\vec {u}}}$ y la temperatura ${\displaystyle T}$, pueden calcularse como los momentos de la función de distribución de la densidad:

${\displaystyle \rho =\int f\,d{\vec {e}},}$

${\displaystyle \rho {\vec {u}}=\int {\vec {e}}f\,d{\vec {e}},}$

${\displaystyle {\frac {\rho DRT}{2}}=\rho \epsilon =\int g\,d{\vec {e}}.}$

El método de Boltzmann en red discretiza esta ecuación limitando el espacio a una red y el espacio de velocidad a un conjunto discreto de velocidades microscópicas (es decir, ${\displaystyle {\vec {e}}_{i}=({\vec {e}}_{ix},{\vec {e}}_{iy})}$). Las velocidades microscópicas en D2Q9, D3Q15 y D3Q19, por ejemplo, se dan como:

${\displaystyle {\vec {e}}_{i}=c\times {\begin{cases}(0,0)&i=0\\(1,0),(0,1),(-1,0),(0,-1)&i=1,2,3,4\\(1,1),(-1,1),(-1,-1),(1,-1)&i=5,6,7,8\\\end{cases}}}$

${\displaystyle {\vec {e}}_{i}=c\times {\begin{cases}(0,0,0)&i=0\\(\pm 1,0,0),(0,\pm 1,0),(0,0,\pm 1)&i=1,2,...,5,6\\(\pm 1,\pm 1,\pm 1)&i=7,8,...,13,14\\\end{cases}}}$

${\displaystyle {\vec {e}}_{i}=c\times {\begin{cases}(0,0,0)&i=0\\(\pm 1,0,0),(0,\pm 1,0),(0,0,\pm 1)&i=1,2,...,5,6\\(\pm 1,\pm 1,0),(\pm 1,0,\pm 1),(0,\pm 1,\pm 1)&i=7,8,...,17,18\\\end{cases}}}$

La ecuación de Boltzmann discretizada en una sola fase para la densidad de masa y la densidad de energía interna son:

${\displaystyle f_{i}({\vec {x}}+{\vec {e}}_{i}\delta _{t},t+\delta _{t})-f_{i}({\vec {x}},t)+F_{i}=\Omega (f),}$

${\displaystyle g_{i}({\vec {x}}+{\vec {e}}_{i}\delta _{t},t+\delta _{t})-g_{i}({\vec {x}},t)+G_{i}=\Omega (g).}$

El operador de colisión se aproxima a menudo mediante un operador de colisión BGK, siempre que también satisfaga las leyes de conservación:

${\displaystyle \Omega (f)={\frac {1}{\tau _{f}}}(f_{i}^{\text{eq}}-f_{i}),}$

${\displaystyle \Omega (g)={\frac {1}{\tau _{g}}}(g_{i}^{\text{eq}}-g_{i}).}$

En el operador de colisión, ${\displaystyle f_{i}^{\text{eq}}}$ es la función de distribución de probabilidad de partículas en equilibrio discreta, equilibrium particle probability distribution function. En D2Q9 y D3Q19, se muestra a continuación para un flujo incompresible en forma continua y discreta, donde «D», «R» y «T» son la dimensión, la constante universal de los gases y la temperatura absoluta, respectivamente. La derivación parcial de la forma continua a la discreta se proporciona mediante una derivación simple con precisión de segundo orden.

${\displaystyle f^{\text{eq}}={\frac {\rho }{(2\pi RT)^{D/2}}}e^{-{\frac {({\vec {e}}-{\vec {u}})^{2}}{2RT}}}}$

${\displaystyle ={\frac {\rho }{(2\pi RT)^{D/2}}}e^{-{\frac {({\vec {e}})^{2}}{2RT}}}e^{{\frac {{\vec {e}}{\vec {u}}}{RT}}-{\frac {{\vec {u}}^{2}}{2RT}}}}$

${\displaystyle ={\frac {\rho }{(2\pi RT)^{D/2}}}e^{-{\frac {({\vec {e}})^{2}}{2RT}}}\left(1+{\frac {{\vec {e}}{\vec {u}}}{RT}}+{\frac {({\vec {e}}{\vec {u}})^{2}}{2(RT)^{2}}}-{\frac {{\vec {u}}^{2}}{2RT}}+...\right)}$

Al sustituir ${\displaystyle c={\sqrt {3RT}}}$, se obtiene el resultado final:

${\displaystyle f_{i}^{eq}=\omega _{i}\rho \left(1+{\frac {3{\vec {e}}_{i}{\vec {u}}}{c^{2}}}+{\frac {9({\vec {e}}_{i}{\vec {u}})^{2}}{2c^{4}}}-{\frac {3({\vec {u}})^{2}}{2c^{2}}}\right)}$

${\displaystyle g^{eq}={\frac {\rho ({\vec {e}}-{\vec {u}})^{2}}{2(2\pi RT)^{D/2}}}e^{-{\frac {({\vec {e}}-{\vec {u}})^{2}}{2RT}}}}$

${\displaystyle \omega _{i}={\begin{cases}4/9&i=0\\1/9&i=1,2,3,4\\1/36&i=5,6,7,8\\\end{cases}}}$

${\displaystyle \omega _{i}={\begin{cases}1/3&i=0\\1/18&i=1,2,...,5,6\\1/36&i=7,8,...,17,18\\\end{cases}}}$

Dado que ya se ha trabajado mucho en el flujo de un solo componente, se analizará el TLBM siguiente. El TLBM multicomponente/multifásico también es más interesante y útil que el de un solo componente. Para estar en consonancia con la investigación actual, definamos el conjunto de todos los componentes del sistema (es decir, paredes de medios porosos, múltiples fluidos/gases, etc.) ${\displaystyle \Psi }$ con elementos ${\displaystyle \sigma _{j}}$.

${\displaystyle f_{i}^{\sigma }({\vec {x}}+{\vec {e}}_{i}\delta _{t},t+\delta _{t})-f_{i}^{\sigma }({\vec {x}},t)+F_{i}={\frac {1}{\tau _{f}^{\sigma }}}(f_{i}^{\sigma ,eq}(\rho ^{\sigma },v^{\sigma })-f_{i}^{\sigma })}$

El parámetro de relajación,${\displaystyle \tau _{f}^{\sigma _{j}}\,\!}$, está relacionado con la viscosidad cinemática,${\displaystyle \nu _{f}^{\sigma _{j}}\,\!}$, mediante la siguiente relación:

${\displaystyle \nu _{f}^{\sigma _{j}}=(\tau _{f}^{\sigma _{j}}-0.5)c_{s}^{2}\delta _{t}.}$

Los momentos de ${\displaystyle f_{i}\,\!}$ dan las cantidades conservadas locales. La densidad viene dada por

${\displaystyle \rho =\sum _{\sigma }\sum _{i}f_{i}\,\!}$

${\displaystyle \rho \epsilon =\sum _{i}g_{i}\,\!}$

${\displaystyle \rho ^{\sigma }=\sum _{i}f_{i}^{\sigma }\,\!}$

y la velocidad media ponderada, ${\displaystyle {\vec {u'}}\,\!}$, y el momento local vienen dados por

${\displaystyle {\vec {u'}}=\left(\sum _{\sigma }{\frac {\rho ^{\sigma }{\vec {u^{\sigma }}}}{\tau _{f}^{\sigma }}}\right)/\left(\sum _{\sigma }{\frac {\rho ^{\sigma }}{\tau _{f}^{\sigma }}}\right)}$

${\displaystyle \rho ^{\sigma }{\vec {u^{\sigma }}}=\sum _{i}f_{i}^{\sigma }{\vec {e}}_{i}.}$

${\displaystyle v^{\sigma }={\vec {u'}}+{\frac {\tau _{f}^{\sigma }}{\rho ^{\sigma }}}{\vec {F}}^{\sigma }}$

En la ecuación anterior para la velocidad de equilibrio ${\displaystyle v^{\sigma }\,\!}$, el término ${\displaystyle {\vec {F}}^{\sigma }\,\!}$ es la fuerza de interacción entre un componente y los demás componentes. Sigue siendo objeto de mucho debate, ya que suele ser un parámetro de ajuste que determina cómo interactúan los fluidos entre sí, los fluidos con los gases, etc. Frank et al. enumeran los modelos actuales para este término de fuerza. Las derivaciones más utilizadas son el modelo cromodinámico de Gunstensen, el enfoque basado en la energía libre de Swift para sistemas líquido/vapor y fluidos binarios, el modelo basado en la interacción intermolecular de He, el enfoque de Inamuro y el enfoque de Lee y Lin.^[18]​

La siguiente es la descripción general de ${\displaystyle {\vec {F}}^{\sigma }\,\!}$ dada por varios autores.^[19]​^[20]​

${\displaystyle {\vec {F}}^{\sigma }=-\psi ^{\sigma }({\vec {x}})\sum _{\sigma _{j}}H^{\sigma \sigma _{j}}({\vec {x}},{\vec {x}}')\sum _{i}\psi ^{\sigma _{j}}({\vec {x}}+{\vec {e}}_{i}){\vec {e}}_{i}\,\!}$

${\displaystyle \psi ({\vec {x}})\,\!}$ es la masa efectiva y ${\displaystyle H({\vec {x}},{\vec {x}}')\,\!}$ es la función de Green que representa la interacción entre partículas con ${\displaystyle {\vec {x}}'\,\!}$ como sitio vecino. Satisfaciendo ${\displaystyle H({\vec {x}},{\vec {x}}')=H({\vec {x}}',{\vec {x}})\,\!}$ y donde ${\displaystyle H({\vec {x}},{\vec {x}}')>0\,\!}$ representa fuerzas repulsivas. Para D2Q9 y D3Q19, esto conduce a

${\displaystyle H^{\sigma \sigma _{j}}({\vec {x}},{\vec {x}}')={\begin{cases}h^{\sigma \sigma _{j}}&\left|{\vec {x}}-{\vec {x}}'\right|\leq c\\0&\left|{\vec {x}}-{\vec {x}}'\right|>c\\\end{cases}}}$ ${\displaystyle H^{\sigma \sigma _{j}}({\vec {x}},{\vec {x}}')={\begin{cases}h^{\sigma \sigma _{j}}&\left|{\vec {x}}-{\vec {x}}'\right|=c\\h^{\sigma \sigma _{j}}/2&\left|{\vec {x}}-{\vec {x}}'\right|={\sqrt {2c}}\\0&{\text{otherwise}}\\\end{cases}}}$

La masa efectiva propuesta por Shan y Chen utiliza la siguiente masa efectiva para un «sistema monocomponente y multifásico». La ecuación de estado también se da en condiciones de un solo componente y multifásico.

${\displaystyle \psi ({\vec {x}})=\psi (\rho ^{\sigma })=\rho _{0}^{\sigma }\left[1-e^{(-\rho ^{\sigma }/\rho _{0}^{\sigma })}\right]\,\!}$

${\displaystyle p=c_{s}^{2}\rho +c_{0}h[\psi ({\vec {x}})]^{2}\,\!}$

Hasta ahora, parece que ${\displaystyle \rho _{0}^{\sigma }\,\!}$ y ${\displaystyle h^{\sigma \sigma _{j}}\,\!}$ son constantes libres que se pueden ajustar, pero una vez introducidas en la ecuación de estado (EOS) del sistema, deben satisfacer las relaciones termodinámicas en el punto crítico, de modo que ${\displaystyle (\partial P/\partial {\rho })_{T}=(\partial ^{2}P/\partial {\rho ^{2}})_{T}=0\,\!}$ y ${\displaystyle p=p_{c}\,\!}$. Para la EOS, ${\displaystyle c_{0}\,\!}$ es 3,0 para D2Q9 y D3Q19, mientras que es igual a 10,0 para D3Q15.^[21]​

Más tarde, Yuan y Schaefer demostraron^[22]​ que la densidad de masa efectiva debe modificarse para simular con mayor precisión el flujo multifásico. Compararon las EOS de Shan y Chen (SC), Carnahan-Starling (C–S), van der Waals (vdW), Redlich–Kwong (R–K), Redlich–Kwong Soave (RKS) y Peng–Robinson (P–R). Sus resultados revelaron que la EOS SC era insuficiente y que las EOS C–S, P–R, R–K y RKS son todas más precisas para modelar el flujo multifásico de un solo componente.

Para los populares métodos isotérmicos de Lattice Boltzmann, estas son las únicas cantidades conservadas. Los modelos térmicos también conservan la energía y, por lo tanto, tienen una cantidad conservada adicional:

${\displaystyle \rho \theta +\rho uu=\sum _{i}f_{i}{\vec {e}}_{i}{\vec {e}}_{i}.}$

Normalmente, los métodos de Boltzmann en red se implementan en mallas regulares. Sin embargo, el uso de mallas no estructuradas puede ayudar a resolver límites complejos. Las mallas no estructuradas están formadas por triángulos o tetraedros con variaciones.

Suponiendo que ${\displaystyle \Omega ^{j}}$ es un volumen formado por todos los baricentros de los tetraedros, caras y aristas conectados al vértice ${\displaystyle {\boldsymbol {v}}^{j}}$, la función de densidad de velocidad discreta:

${\displaystyle f_{i}({\boldsymbol {v}}^{j},t+\delta t)=f_{i}({\boldsymbol {v}}^{j},t)-\delta t\sum _{k}S_{i}^{jk}f_{i}({\boldsymbol {v}}^{k},t)-{\delta t \over \tau }\sum _{k}C^{jk}(f_{i}({\boldsymbol {v}}^{k},t)-f_{i}^{eq}({\boldsymbol {v}}^{k}))}$

donde ${\displaystyle {\boldsymbol {v}}^{k}}$ son la posición de un vértice y sus vecinos, y:

${\displaystyle C^{jk}={1 \over V^{j}}\int _{\Omega ^{j}}w_{k}({\boldsymbol {x}})d\Omega }$

${\displaystyle S_{i}^{jk}={1 \over V^{j}}\oint _{\partial \Omega ^{j}}({\vec {e_{i}}}{\vec {n}})w_{k}({\boldsymbol {x}})d\Omega }$

donde ${\displaystyle w_{k}({\boldsymbol {x}})}$ son los pesos de una interpolación lineal de ${\displaystyle {\boldsymbol {x}}}$ por los vértices del triángulo o tetraedro en el que se encuentra ${\displaystyle {\boldsymbol {x}}}$.^[23]​

Durante los últimos años, el LBM ha demostrado ser una herramienta poderosa para resolver problemas a diferentes escalas de longitud y tiempo. Algunas de las aplicaciones del LBM incluyen:

• Flujos en medios porosos ^[24]​
• Flujos biomédicos
• Ciencias de la Tierra (filtración del suelo).
• Ciencias energéticas (pilas de combustible^[25]​).

Esta es una implementación básica de LBM en una cuadrícula de 100x100, utilizando Python:

#Este es un simulador de fluidos que utiliza el método de Boltzmann en red.
#Se utilizó D2Q9 y límites periódicos, sin utilizar ninguna biblioteca externa.
#Genera dos ondas a 50,50 y 50,40.
#Reference: Erlend Magnus Viggen's Master thesis, "The Lattice Boltzmann Method with Applications in Acoustics".
#For Wikipedia under CC-BY-SA license.
import math
#Define some utilities
def sum(a):
    s=0
    for e in a:
        s=s+e
    return s
#Weights in D2Q9
Weights=[1/36,1/9,1/36,
   1/9, 4/9,1/9,
   1/36,1/9,1/36]
#Discrete velocity vectors
DiscreteVelocityVectors=[[-1,1],[0,1],[1,1],
     [-1,0],[0,0],[1,0],
     [-1,-1],[0,-1],[1,-1]
]
#A Field2D class
class Field2D():
    def __init__(self,res : int):
        self.field=[]
        for b in range(res):
            fm=[]
            for a in range(res):
                fm.append([0,0,0,
                           0,1,0,
                           0,0,0])
            self.field.append(fm[:])
        self.res = res
    #This visualize the simulation, can only be used in a terminal
    @staticmethod
    def VisualizeField(a,sc,res):
        stringr=""
        for u in range(res):
            row=""
            for v in range(res):
                n=int(u*a.res/res)
                x=int(v*a.res/res)
                flowmomentem=a.Momentum(n,x)
                col="\033[38;2;{0};{1};{2}m██".format(int(127+sc*flowmomentem[0]),int(127+sc*flowmomentem[1]),0)
                row=row+col
            print(row)
            stringr=stringr+row+"\n"
        return stringr
    #Momentum of the field
    def Momentum(self,x,y):
        return velocityField[y][x][0]*sum(self.field[y][x]),velocityField[y][x][1]*sum(self.field[y][x])
#Resolution of the simulation
res=100
a=Field2D(res)
#The velocity field
velocityField=[]
for DummyVariable in range(res):
    DummyList=[]
    for DummyVariable2 in range(res):
        DummyList.append([0,0])
    velocityField.append(DummyList[:])
#The density field
DensityField=[]
for DummyVariable in range(res):
    DummyList=[]
    for DummyVariable2 in range(res):
        DummyList.append(1)
    DensityField.append(DummyList[:])
#Set initial condition
DensityField[50][50]=2
DensityField[40][50]=2
#Maximum solving steps
MaxSteps = 120
#The speed of sound, specifically 1/sqrt(3) ~ 0.57
SpeedOfSound=1/math.sqrt(3)
#time relaxation constant
TimeRelaxationConstant=0.5
#Solve
for s in range(MaxSteps):
    #Collision Step
    df=Field2D(res)
    for y in range(res):
        for x in range(res):
            for v in range(9):
                Velocity=a.field[y][x][v]
                FirstTerm=Velocity
                #The Flow Velocity
                FlowVelocity=velocityField[y][x]
                Dotted=FlowVelocity[0]*DiscreteVelocityVectors[v][0]+FlowVelocity[1]*DiscreteVelocityVectors[v][1]
                # #The taylor expainsion of equilibrium term
                taylor=1+((Dotted)/(SpeedOfSound**2))+((Dotted**2)/(2*SpeedOfSound**4))-((FlowVelocity[0]**2+FlowVelocity[1]**2)/(2*SpeedOfSound**2))
                #The current density
                density=DensityField[y][x]
                #The equilibrium
                equilibrium=density*taylor*Weights[v]
                SecondTerm=(equilibrium-Velocity)/TimeRelaxationConstant
                df.field[y][x][v]=FirstTerm+SecondTerm
    #Streaming Step
    for y in range(0,res):
        for x in range(0,res):
            for v in range(9):
                #Target, the lattice point this iteration is solving
                TargetY=y+DiscreteVelocityVectors[v][1]
                TargetX=x+DiscreteVelocityVectors[v][0]
                # Peiodic Boundary
                if TargetY == res and TargetX == res:
                    a.field[TargetY-res][TargetX-res][v]=df.field[y][x][v]
                elif TargetX == res:
                    a.field[TargetY][TargetX-res][v]=df.field[y][x][v]
                elif TargetY == res:
                    a.field[TargetY-res][TargetX][v]=df.field[y][x][v]
                elif TargetY == -1 and TargetX == -1:
                    a.field[TargetY+res][TargetX+res][v]=df.field[y][x][v]   
                elif TargetX == -1:
                    a.field[TargetY][TargetX+res][v]=df.field[y][x][v]
                elif TargetY == -1:
                    a.field[TargetY+res][TargetX][v]=df.field[y][x][v]
                else:
                    a.field[TargetY][TargetX][v]=df.field[y][x][v]
    #Calculate macroscopic variables
    for y in range(res):
        for x in range(res):
            #Recompute Density Field
            DensityField[y][x]=sum(a.field[y][x])
            #Recompute Flow Velocity
            FlowVelocity=[0,0]
            for DummyVariable in range(9):
                FlowVelocity[0]=FlowVelocity[0]+DiscreteVelocityVectors[DummyVariable][0]*a.field[y][x][DummyVariable]
            for DummyVariable in range(9):
                FlowVelocity[1]=FlowVelocity[1]+DiscreteVelocityVectors[DummyVariable][1]*a.field[y][x][DummyVariable]
            FlowVelocity[0]=FlowVelocity[0]/DensityField[y][x]
            FlowVelocity[1]=FlowVelocity[1]/DensityField[y][x]
            #Insert to Velocity Field
            velocityField[y][x]=FlowVelocity
    #Visualize
    Field2D.VisualizeField(a,128,100)

• Deutsch, Andreas; Sabine Dormann (2004). Cellular Automaton Modeling of Biological Pattern Formation. Birkhäuser Verlag. ISBN 978-0-8176-4281-5. 
• Succi, Sauro (2001). The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. Oxford University Press. ISBN 978-0-19-850398-9. 
• Wolf-Gladrow, Dieter (2000). Lattice-Gas Cellular Automata and Lattice Boltzmann Models. Springer Verlag. ISBN 978-3-540-66973-9. 
• Sukop, Michael C.; Daniel T. Thorne, Jr. (2007). Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers. Springer. ISBN 978-3-540-27981-5. 
• Jian Guo Zhou (2004). Lattice Boltzmann Methods for Shallow Water Flows. Springer. ISBN 978-3-540-40746-1. 
• He, X., Chen, S., Doolen, G. (1998). A Novel Thermal Model for the Lattice Boltzmann Method in Incompressible Limit. Academic Press. 
• Guo, Z. L.; Shu, C (2013). Lattice Boltzmann Method and Its Applications in Engineering. World Scientific Publishing. 
• Huang, H.; M.C. Sukop; X-Y. Lu (2015). Multiphase Lattice Boltzmann Methods: Theory and Application. Wiley-Blackwell. ISBN 978-1-118-97133-8. 
• Krüger, T.; Kusumaatmaja, H.; Kuzmin, A.; Shardt, O.; Silva, G.; Viggen, E. M. (2017). The Lattice Boltzmann Method: Principles and Practice. Springer Verlag. ISBN 978-3-319-44647-9. 

• LBM Method
• Entropic Lattice Boltzmann Method (ELBM)
• dsfd.org: Website of the annual DSFD conference series (1986 -- now) where advances in theory and application of the lattice Boltzmann method are discussed
• Website of the annual ICMMES conference on lattice Boltzmann methods and their applications