不太一样的CFD方法：格子玻尔兹曼方法LBM实例
格子玻尔兹曼方法（Lattice Boltzmann Method，简称LBM）是一种用于模拟流体流动的计算方法。
和传统的流体力学不同的是，LBM基于玻尔兹曼方程，通过将连续的流体宏观行为离散化为微小的空间格子和时间步骤，从而简化了流体动力学问题的求解过程。
格子玻尔兹曼方法的基本思想和步骤：
格子表示空间： 想象整个流体区域被划分成许多小格子，就像在一张网格纸上一样。每个格子内都有一些微小的粒子，这些粒子代表了流体中的微观粒子。
时间步骤模拟： 时间被划分成小的步骤。在每个时间步骤内，粒子按照事先定义的规则在格子之间移动，这个规则基于碰撞和传播的原理。
碰撞模拟： 粒子在每个格子内相互碰撞，模拟了流体的宏观行为。这个步骤使用玻尔兹曼方程的一部分来更新粒子的速度和分布函数。
转移模拟： 粒子根据规定的速度在格子之间传播，这模拟了流体在空间中的运动。这一步使得流体在每个时间步骤中都有一定的流动。
边界处理： 处理流体与边界的交互，确保在边界上符合物理条件，比如反射、吸收或施加外部条件。
我们可以通过一个类比来更容易理解格子玻尔兹曼方法。
在一片水域中，水域被划分成了很多小方块，每个小方块内都有一群小鱼，它们代表了流体中微观的粒子。我们想要模拟这片水域中的水流动态，但是一下子考虑所有的水分子可能会非常复杂，类比到LBM方法：
小方块表示空间： 把整个水域划分成小方块，每个小方块代表了一个微小的区域，就像在游戏棋盘上一样。
时间步骤模拟： 时间被分成小的步骤，就好像在游戏中每过一轮。在每一轮内，小鱼（粒子）按照一定的规则在小方块之间移动。
碰撞模拟： 在每一轮中，小鱼在小方块内相互碰撞，这样模拟了水的宏观行为，比如水分子之间的相互作用。这个碰撞的过程使得小鱼的速度和分布发生变化。
转移模拟： 小鱼根据预定的速度在小方块之间传播，就像水在空间中流动一样。这使得水在每一轮中都有了一定的流动。
边界处理： 处理水域与边界的交互，确保在水域边界上符合物理条件，比如小鱼不会跑到陆地上去。
通过不断重复这些步骤，我们就可以观察到整个水域中水流的模拟。格子玻尔兹曼方法就好比是通过模拟小方块内小鱼的行为，来近似描述整个水域的流动情况。
下面给出一个代码例子：充满空气的方形空腔，从左侧壁加热，从右侧壁冷却，通过自然对流的方式来驱动空气流动。
代码使用Matlab书写，复 制到Matlab中直接运行即可，代码是我自己写的，可以正常运行，存在可以优化的地方。
运行得到的温度等值线图：
clc;clear;%充满空气的方形空腔（Pr = 0.7），从左侧壁加热，从右侧壁冷却。其他边界底部和顶部是绝热的% 对于Ra = 10^5，确定等温线、努塞特数和流线。%分析：Ra数=10^5。在LBM中格子粘度<=0.1% 假设格子粘度=0.03，pr=v/α=0.7（air）,α=0.043，g_beta=Ra*α*v/(ΔT*H^3)% 格子数等于NX=NY=100% 边界条件：上下绝热；左温度=1，右温度=0；%modelLX=0.2;LY=0.2;%meshNX=100;NY=100;dx=LX/NX;dy=LY/NY;y=1:NY+1;[X,Y]=meshgrid(1:NX+1,1:NY+1);%lattice constantw=[4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36];c=[0,0;1,0;0,1;-1,0;0,-1;1,1;-1,1;-1,-1;1,-1];%alpha和U的取值会影响收敛性alpha_v=0.03;%格子粘度U=0.1;Ra=10.e5;%瑞利数alpha_T=0.043;%格子热扩散系数T_l=1.0;T_r=0.;Pr=alpha_v/alpha_T;%格子普朗特数g_beta=Ra*alpha_T*alpha_v/((T_l-T_r)*NY^3);omega_v=1/(3*alpha_v+0.5);%速度松弛时间倒数omega_T=1/(3*alpha_T+0.5);%温度松弛时间倒数cs2=1/3;%声明数组变量f=zeros(NY+1,NX+1,9);%动量f_eq=zeros(NY+1,NX+1,9);rho_f=ones(NY+1,NX+1);u=zeros(NY+1,NX+1,2);u(NY+1,:,1)=U;g=zeros(NY+1,NX+1,9);%温度g_eq=zeros(NY+1,NX+1,9);rho_g=ones(NY+1,NX+1);T=zeros(NY+1,NX+1);%% 初始化t1=u(:,:,1).*u(:,:,1)+u(:,:,2).*u(:,:,2);%u*ufor j=1:9force=3*w(j)*g_beta*T(:,:)*c(j,2);t2=c(j,1)*u(:,:,1)+c(j,2)*u(:,:,2);%ck*uf_eq(:,:,j)=w(j)*rho_f(:,:).*(1+t2/cs2+0.5*t2.*t2/cs2^2-0.5*t1/cs2);f(:,:,j)=f(:,:,j)*(1-omega_v)+omega_v*f_eq(:,:,j)+force;g_eq(:,:,j)=w(j)*T(:,:).*(1+t2/cs2);enditer=50000;for i=1:iter%% collsiont1=u(:,:,1).*u(:,:,1)+u(:,:,2).*u(:,:,2);%u*ufor j=1:9force=3*w(j)*g_beta*T(:,:)*c(j,2);t2=c(j,1)*u(:,:,1)+c(j,2)*u(:,:,2);%ck*uf_eq(:,:,j)=w(j)*rho_f(:,:).*(1+t2/cs2+0.5*t2.*t2/cs2^2-0.5*t1/cs2);f(:,:,j)=f(:,:,j)*(1-omega_v)+omega_v*f_eq(:,:,j)+force;g_eq(:,:,j)=w(j)*T(:,:).*(1+t2/cs2);g(:,:,j)=g(:,:,j)*(1-omega_T)+omega_T*g_eq(:,:,j);end%% streamingf(:,2:NX+1,2)=f(:,1:NX,2);f(:,1:NX,4)=f(:,2:NX+1,4);f(2:NY+1,:,3)=f(1:NY,:,3);f(1:NY,:,5)=f(2:NY+1,:,5);f(2:NY+1,2:NX+1,6)=f(1:NY,1:NX,6);f(2:NY+1,1:NX,7)=f(1:NY,2:NX+1,7);f(1:NY,1:NX,8)=f(2:NY+1,2:NX+1,8);f(1:NY,2:NX+1,9)=f(2:NY+1,1:NX,9);g(:,2:NX+1,2)=g(:,1:NX,2);g(:,1:NX,4)=g(:,2:NX+1,4);g(2:NY+1,:,3)=g(1:NY,:,3);g(1:NY,:,5)=g(2:NY+1,:,5);g(2:NY+1,2:NX+1,6)=g(1:NY,1:NX,6);g(2:NY+1,1:NX,7)=g(1:NY,2:NX+1,7);g(1:NY,1:NX,8)=g(2:NY+1,2:NX+1,8);g(1:NY,2:NX+1,9)=g(2:NY+1,1:NX,9);%% boundary of v%top为wall，反弹边界条件f(NY+1,:,5)=f(NY+1,:,3);f(NY+1,:,8)=f(NY+1,:,6);f(NY+1,:,9)=f(NY+1,:,7);%bottom为wall，反弹边界条件f(1,:,3)=f(1,:,5);f(1,:,6)=f(1,:,8);f(1,:,7)=f(1,:,9);%right为wall，反弹边界条件f(:,NX+1,4)=f(:,NX+1,2);f(:,NX+1,7)=f(:,NX+1,9);f(:,NX+1,8)=f(:,NX+1,6);%left为为wall，反弹边界条件f(:,1,2)=f(:,1,4);f(:,1,6)=f(:,1,8);f(:,1,9)=f(:,1,7);%% boundary condition of T%top 反弹格式g(NY+1,:,5)=g(NY+1,:,3);g(NY+1,:,8)=g(NY+1,:,6);g(NY+1,:,9)=g(NY+1,:,7);%leftg(:,1,2)=w(2)*T_l+w(4)*T_l-g(:,1,4);g(:,1,6)=w(6)*T_l+w(8)*T_l-g(:,1,8);g(:,1,9)=w(9)*T_l+w(7)*T_l-g(:,1,7);%rightg(:,NX+1,4)=w(4)*T_r+w(2)*T_r-g(:,NX+1,2);g(:,NX+1,7)=w(7)*T_r+w(9)*T_r-g(:,NX+1,9);g(:,NX+1,8)=w(8)*T_r+w(6)*T_r-g(:,NX+1,6);%bottomg(1,:,3)=g(1,:,5);g(1,:,7)=g(1,:,9);g(1,:,6)=g(1,:,8);T=sum(g,3);rho_f=sum(f,3);u=zeros(NY+1,NX+1,2);for m=1:9u(:,:,1)=u(:,:,1)+c(m,1)*f(:,:,m)./rho_f;u(:,:,2)=u(:,:,2)+c(m,2)*f(:,:,m)./rho_f;endu(NY+1,2:NX,1)=U;u(NY+1,2:NX,2)=0;end%plot(T(:,25),y,T(:,50),y,T(:,75),y);%legend("25","50","75");figurecontourf(X,Y,u(:,:,1).*u(:,:,1)+u(:,:,2).*u(:,:,2));figurecontour(T,30);
