力驱动泊萧叶流，基于 MRT-LBM 的流场与声场仿真计算 --王富海2017

%这个例子采用 MRT-LBM 模拟体积力驱动的二维泊肃叶流，
%左右周期边界格式，上下非平衡反弹格式.
%基于 MRT-LBM 的流场与声场仿真计算 --王富海2017

%还是老规矩先宣传一下QQ群： 格子玻尔兹曼救星：293267908。
clc
clear
close all
 
%设置仿真参数
xLen=60; %水平方向格子数
yLen=30; %垂直方向格子数
Re=10;%雷诺数
uMax=0.1;%中间最大速度
uAve=2/3*uMax; %泊肃叶流解析解的平均流速
nu=uAve*(yLen-1)/Re;%运动粘性
Cs=sqrt(1/3);   %格子声速
tau=1/2+3*nu; %松弛时间
omega=1/tau;%松弛频率
step=10000; %最大迭代次数
y1=1;
y2=yLen;
GG=uMax/((y1-y2)^2/4);
f_u=GG*(2*nu);%常外力项,每个节点的力都是一样的
f_v=0; %常外力项,每个节点的力都是一样的
rhoo=1; %初始化密度
checkStep=100;%收敛计算间隔
saveStep=20; %保存结果间隔
filePath=uigetdir('*.*','D:\MyStudy\MATLAB\YuBrian');%仿真中间过程图片的保存路径
VSSum=[];   %所有节点格子速度总和-每 save Step 步记录一次
VSSum2=[]; %所有节点格子速度总和-每 check Step 步记录一次，监视收敛曲线
%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