这个例子采用 MRT-LBM 模拟矩形腔绕流--基于 MRT-LBM 的流场与声场仿真计算 --王富海2017

%这个例子采用 MRT-LBM 模拟矩形腔绕流
%基于 MRT-LBM 的流场与声场仿真计算 --王富海2017
%上边界速度边界，其它边界-上下非平衡反弹格式-无滑移壁面
%还是老规矩先宣传一下QQ群： 格子玻尔兹曼救星：293267908。
clc
clear
close all
 
%设置仿真参数
 uMax=0.1;%中间最大速度
 xLen=151; %水平方向格子数   
 yLen=101;%垂直方向格子数
 Re=100; %雷诺数
 nu=uMax*(xLen-1)/Re; %运动粘性
 Cs=sqrt(1/3);%格子声速
 tau=1/2+3*nu;%松弛时间
 omega=1/tau;%松弛频率
 step=1000;%最大迭代次数
 rhoo=1;%初始化密度
 checkStep=100;%收敛计算间隔
 saveStep=20;%保存结果间隔
 filePath=uigetdir('*.*','D:\MyStudy\MATLAB\YuBrian')%仿真中间过程图片的保存路径
 VSSum=[];%所有节点格子速度总和-每 saveStep 步记录一次
 VSSum2=[]; %所有节点格子速度总和-每 checkStep 步记录一次，监视收敛曲线
 %D2Q9 模型参数
 w=[1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36 4/9 ]%各个方向的权重
 cx=[ 1 0 -1 0 1 -1   -1   1 0];%各方向 x 速度分量
 cy=[ 0 1 0 -1 1 1 -1 -1 0];%各方向 y 速度分量
 M=[1 1 1 1 1 1 1 1 1;-4 -1 -1 -1 -1 2 2 2 2;4 -2 -2 -2 -2 1 1 1 1; ...
        0 1 0 -1 0 1 -1 -1 1;0 -2 0 2 0 1 -1 -1 1;0 0 1 0 -1 1 1 -1 -1; ...
        0 0 -2 0 2 1 1 -1 -1;0 1 -1 1 -1 0 0 0 0;0 0 0 0 0 1 -1 1 -1;]
 %  由于分布函数 0 索引被放置在 matlab  索引 9 的位置，所以将 M 第一列放到最后。
 
M=circshift(M,[-1 -1]);                                                                  
        s1(9)=0;
        s1(3)=s1(9);
          s1(5)=s1(9);
          s1(7)=omega;
          s1(8)=s1(7);
          s1(4)=1.2;
          s1(6)=s1(4);
          s1(1)=omega;
          s1(2)=s1(1)-0.1;
 s=diag(s1);%对角矩阵-松弛因子
 Minv=inv(M);
 Minv_s=Minv*s;
 
 
%  网格节点初始化，初始化各节点密度=1  ，初始化各节点速度为 0，初始化各节点分布函数=平衡分布函数
 
for i=1:xLen
      for j=1:yLen
          rho(i,j)=rhoo;