学习LBM的第一个代码

这是我学习LBM时接触到的第一个代码，是一个“祖传”的代码，在很多博客上都能找到，我按照自己初学时的理解加了注释。这个代码是关于D2Q9圆柱绕流模拟的。b站有个up主讲解过这个代码，这是网站链接：https://www.bilibili.com/video/BV1nZ4y1h7DF/?spm_id_from=333.999.0.0&vd_source=9ab5ddb7b0de0c8f2730f37d076c851d

%D2Q9圆柱绕流模拟 
clear all 
close all 
clc 

% 一些常数 --------------------------------------------------
lx = 400;                 % x方向格子数
ly = 100;                 % y方向格子数
obst_x = lx/5+1;          % 圆柱位置
obst_y = ly/2+3;                         
obst_r = ly/10+1;         % 圆柱半径
uMax = 0.1;               % Poiseuille流入的最大速度
Re = 100;                 % 雷诺数
nu = uMax * 2.*obst_r / Re; % 运动粘度
omega = 1. / (3*nu+1./2.);  % 松弛参数
maxT = 400000;              % 迭代数
tPlot = 50;                 % cycles
 
% =========================================================================
% D2Q9格子常数 
t     = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36]; %权重
cx    = [ 0, 1, 0, -1, 0, 1, -1, -1, 1]; %速度离散的9个方向
cy    = [ 0, 0, 1, 0, -1, 1, 1, -1, -1]; 
opp   = [ 1, 4, 5, 2, 3, 8, 9, 6, 7]; %反弹边界条件要用到，2、4互换，3、5互换等
col   = [2:(ly-1)]; 
in    = 1;                  % 流入位置
out   = lx;                 % 流出位置
[y,x] = meshgrid(1:ly,1:lx);    % 获取矩阵索引的坐标
obst  = (x-obst_x).^2 + (y-obst_y).^2 <= obst_r.^2;% 判断是否在圆柱内，圆柱内为1

obst(:,[1,ly]) = 1;         % 边界为1，标识边界
bbRegion = find(obst);      % obst在圆柱内的索引
 
 
% =========================================================================
% 初始条件INITIAL CONDITION: Poiseuille profile at equilibrium --------------------
L      = ly-2;  %98，长度99，物理边界从0.5，99.5开始
y_phys = y-1.5; %会生成一个矩阵，实际的物理边界
%初始速度，二次函数，大红书P185
ux = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys); 
uy = zeros(lx,ly); 
rho = 1; 
 
%宏观转介观
for i=1:9 
    cu = 3*(cx(i)*ux+cy(i)*uy);%方便写代码  
    %初始分布速度直接设置成平衡态，根据初始速度确定平衡分布函数，后面再根据速度改变并趋向最终的平衡
    %这样设置达到平衡快 
    fIn(i,:,:) = rho .* t(i) .*( 1 + cu + 1/2*(cu.*cu) - 3/2*(ux.^2+uy.^2) ); %书本3.4式
end 
 
% 主循环（时间循环）
for cycle = 1:maxT 
    %介观转宏观
    % 宏观量
    rho = sum(fIn); %密度
    ux = reshape ( (cx * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho; %速度，先变成二维，再reshape回去
    uy = reshape ( (cy * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho; 
 
    % 宏观 (DIRICHLET) 边界条件 ------------------------- 
    % Inlet: Poiseuille profile -------------------------------------------
    y_phys = col-1.5; 
    ux(:,in,col) = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys); %前面出现过
    uy(:,in,col) = 0; 
    rho(:,in,col) = 1 ./ (1-ux(:,in,col)) .* (sum(fIn([1,3,5],in,col)) + 2*sum(fIn([4,7,8],in,col))); 
 
    % Outlet: Constant pressure -------------------------------------------
    rho(:,out,col) = 1; 
    ux(:,out,col) = -1 + 1 ./ (rho(:,out,col)) .* (sum(fIn([1,3,5],out,col)) + 2*sum(fIn([2,6,9],out,col))); 
    uy(:,out,col) = 0; 
    % 微观边界条件:  (Zou/He BC 非平衡态反弹)，in代表入边界，out代表出边界。大红书198页
    %有了宏观条件才有微观条件
    fIn(2,in,col) = fIn(4,in,col) + 2/3*rho(:,in,col).*ux(:,in,col); 
    fIn(6,in,col) = fIn(8,in,col) + 1/2*(fIn(5,in,col)-fIn(3,in,col)) ...
                                  + 1/2*rho(:,in,col).*uy(:,in,col) ...
                                  + 1/6*rho(:,in,col).*ux(:,in,col); 
                               
    fIn(9,in,col) = fIn(7,in,col) + 1/2*(fIn(3,in,col)-fIn(5,in,col)) ...
                                  - 1/2*rho(:,in,col).*uy(:,in,col) ...
                                  + 1/6*rho(:,in,col).*ux(:,in,col); 
      
    % 微观边界条件: OUTLET (Zou/He BC)。大红书198页
    fIn(4,out,col) = fIn(2,out,col) - 2/3*rho(:,out,col).*ux(:,out,col); 
 
    fIn(8,out,col) = fIn(6,out,col) + 1/2*(fIn(3,out,col)-fIn(5,out,col)) ...
                                    - 1/2*rho(:,out,col).*uy(:,out,col) ...
                                    - 1/6*rho(:,out,col).*ux(:,out,col); 
 
    fIn(7,out,col) = fIn(9,out,col) + 1/2*(fIn(5,out,col)-fIn(3,out,col)) ...
                                    + 1/2*rho(:,out,col).*uy(:,out,col) ...
                                    - 1/6*rho(:,out,col).*ux(:,out,col); 
% 碰撞（内部碰撞） ------------------------------------------------------
for i=1:9 
    cu = 3*(cx(i)*ux+cy(i)*uy); 
    fEq(i,:,:) = rho .* t(i) .* ( 1 + cu + 1/2*(cu.*cu) - 3/2*(ux.^2+uy.^2) ); 
    %书本3.13式
    fOut(i,:,:) = fIn(i,:,:) - omega .* (fIn(i,:,:)-fEq(i,:,:)); %fOut对应fIn,fOut是新的fIn
end 
 
% =========================================================================
% OBSTACLE (BOUNCE-BACK)（边界反弹）
for i=1:9 
    fOut(i,bbRegion) = fIn(opp(i),bbRegion); 
end 
 
% =========================================================================
% 迁移，书本3.3.3.2节有写
for i=1:9 
    fIn(i,:,:) = circshift(fOut(i,:,:), [0,cx(i),cy(i)]); 
end 
 
% =========================================================================
% 可视化
if (mod(cycle,tPlot)==1) 
    u = reshape(sqrt(ux.^2+uy.^2),lx,ly); %变成二维的
    u(bbRegion) = nan; %标为1的边界及圆柱区域
    imagesc(u');
    axis equal off; drawnow
end
end