这个例子采用 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;
u(i,j)=0;
v(i,j)=0;
t1=u(i,j)^2+v(i,j)^2; %速度 U*V 内积,t1 写到外面，节省计算量
for k=1:9
t2=u(i,j)*cx(k)+v(i,j)*cy(k);%U*C 内积
f(k,i,j)=rho(i,j)*w(k)*(1+3*t2+4.5*t2^2-1.5*t1);% 计算初始态下的 分布函数=平衡分布函数
end
end
end
%主循环
tstart = tic;%计算时钟
[meshX,meshY]=meshgrid(1:xLen,yLen:-1:1);%产生网格数据。
figure
pc=pcolor(meshX,meshY,ones(yLen,xLen));%绘制伪彩图
shading interp%伪彩色图
axis equal
ylim([1 yLen])
%set(gca,'clim',[0 0.25])
for kk=1:step
step2=kk%显示程序的步进
mm=max(max(rho))%观察最大密度
if mm>100%判断仿真是否发散
break;
end
%1 碰撞过程
for i=1:xLen
for j=1:yLen %
%速度空间到动量空间
m=M*f(1:9,i,j);
%计算平衡动量
m_eq(9,1)=rho(i,j);
m_eq(1,1)=rho(i,j)*(-2+3*(u(i,j)^2+v(i,j)^2));
m_eq(2,1)=rho(i,j)*(1-3*(u(i,j)^2+v(i,j)^2));
m_eq(3,1)=rho(i,j)*u(i,j);
m_eq(4,1)=-rho(i,j)*u(i,j);
m_eq(5,1)=rho(i,j)*v(i,j);
m_eq(6,1)=-rho(i,j)*v(i,j);
m_eq(7,1)=rho(i,j)*(u(i,j)^2-v(i,j)^2);
m_eq(8,1)=rho(i,j)*u(i,j)*v(i,j);
%在动量空间松弛
m_temp=Minv_s*(m-m_eq);
for k=1:9
f(k,i,j)=f(k,i,j)-m_temp(k);
end
end
end
% 2 迁移过程
f2=f;%将碰撞后分布函数保存给变量 f2
for i=1:9
f(i,:,: ) = circshift(f2(i,:,: ), [0,cx(i),cy(i)]);
end
%施加上下边界条件-非平衡反弹模式
u(1:xLen,1)=0;
v(1:xLen,1)=0;
u(1:xLen,yLen)=uMax;
v(1:xLen,yLen)=0;
for i=2:xLen-1
%下边界
rho(i,1)=(f(9,i,1)+f(1,i,1)+f(3,i,1)+2*(f(4,i,1)+f(7,i,1)+f(8,i,1)))/(1-v(i,1));
f(2,i,1)=f(4,i,1)+2/3*rho(i,1)*v(i,1);
f(5,i,1)=f(7,i,1)-1/2*(f(1,i,1)-f(3,i,1))+1/6*rho(i,1)*v(i,1)+1/2*rho(i,1)*u(i,1);
f(6,i,1)=f(8,i,1)+1/2*(f(1,i,1)-f(3,i,1))+1/6*rho(i,1)*v(i,1)-1/2*rho(i,1)*u(i,1);
%上边界
rho(i,yLen)=(f(9,i,yLen)+f(1,i,yLen)+f(3,i,yLen)+2*(f(2,i,yLen)+f(6,i,yLen)+f(5,i,yLen)))/(1+v(i,yLen));
f(4,i,yLen)=f(2,i,yLen)-2/3*rho(i,yLen)*v(i,yLen);
f(7,i,yLen)=f(5,i,yLen)+1/2*(f(1,i,yLen)-f(3,i,yLen))-1/6*rho(i,yLen)*v(i,yLen)-1/2*rho(i,yLen)*u(i,yLen);
f(8,i,yLen)=f(6,i,yLen)-1/2*(f(1,i,yLen)-f(3,i,yLen))-1/6*rho(i,yLen)*v(i,yLen)+1/2*rho(i,yLen)*u(i,yLen);
end
%施加左右边界条件-非平衡反弹模式
u(1,:)=0;
v(1,:)=0;
u(xLen,:)=0;
v(xLen,:)=0;
for j=2:yLen-1
%左边界-速度边界
rho(1,j)=(f(9,1,j)+f(2,1,j)+f(4,1,j)+2*(f(3,1,j)+f(6,1,j)+f(7,1,j)))/(1-u(1,j));
%rho(1,j)=1;
f(1,1,j)=f(3,1,j)+2/3*rho(1,j)*u(1,j);
f(5,1,j)= f(7,1,j)-1/2*( f(2,1,j)-f(4,1,j))+1/6*rho(1,j)*u(1,j)+1/2*rho(1,j)*v(1,j);
f(8,1,j)= f(6,1,j)+1/2*( f(2,1,j)-f(4,1,j))+1/6*rho(1,j)*u(1,j)-1/2*rho(1,j)*v(1,j);
%右边界-速度边界
rho(xLen,j)=(f(9,xLen,j)+f(2,xLen,j)+f(4,xLen,j)+2*(f(1,xLen,j)+f(5,xLen,j)+f(8,xLen,j)))/(1+u(xLen,j));
f(3,xLen,j)=f(1,xLen,j)-2/3*rho(xLen,j)*u(xLen,j);
f(7,xLen,j)= f(5,xLen,j)+1/2*( f(2,xLen,j)-f(4,xLen,j))-1/6*rho(xLen,j)*u(xLen,j)-1/2*rho(xLen,j)*v(xLen,j);
f(6,xLen,j)= f(8,xLen,j)-1/2*( f(2,xLen,j)-f(4,xLen,j))-1/6*rho(xLen,j)*u(xLen,j)-1/2*rho(xLen,j)*v(xLen,j);
end
%四个角落的处理
%左上角
f(1,1,yLen)=f(3,1,yLen);
f(4,1,yLen)=f(2,1,yLen);
f(8,1,yLen)=f(6,1,yLen);
f(5,1,yLen)=0.5*(rho(1,yLen-1)-f(9,1,yLen)-f(1,1,yLen)-f(2,1,yLen)-f(3,1,yLen)-f(4,1,yLen)-f(6,1,yLen)-f(8,1,yLen));
f(7,1,yLen)=f(5,1,yLen);
%左下角
f(1,1,1)=f(3,1,1);
f(2,1,1)=f(4,1,1);
f(5,1,1)=f(7,1,1);
f(6,1,1)=0.5*(rho(1,2)-f(9,1,1)-f(1,1,1)-f(2,1,1)-f(3,1,1)-f(4,1,1)-f(5,1,1)-f(7,1,1));
f(8,1,1)=f(6,1,1);
%右上角
f(3,xLen,yLen)=f(1,xLen,yLen);
f(4,xLen,yLen)=f(2,xLen,yLen);
f(7,xLen,yLen)=f(5,xLen,yLen);
f(6,xLen,yLen)=0.5*(rho(xLen,yLen-1)-f(9,xLen,yLen)-f(1,xLen,yLen)-f(2,xLen,yLen) ...
-f(3,xLen,yLen)-f(4,xLen,yLen)-f(5,xLen,yLen)-f(7,xLen,yLen));
f(8,xLen,yLen)=f(6,xLen,yLen);
%右下角
f(3,xLen,1)=f(1,xLen,1);
f(2,xLen,1)=f(4,xLen,1);
f(6,xLen,1)=f(8,xLen,1);
f(5,xLen,1)=0.5*(rho(xLen,2)-f(9,xLen,1)-f(1,xLen,1)-f(2,xLen,1) -f(3,xLen,1)-f(4,xLen,1)-f(6,xLen,1)-f(6,xLen,1));
f(7,xLen,1)=f(5,xLen,1);
%计算宏观量
for i=1:xLen
for j=1:yLen
rho(i,j)=sum(f(:,i,j));%计算每一点的密度
end
end
for i=1:xLen
for j=1:yLen
uSum=0;
vSum=0;
for k=1:9
uSum=uSum+f(k,i,j)*cx(k);
vSum=vSum+f(k,i,j)*cy(k);
end
u(i,j)=uSum/rho(i,j);%计算速度场
v(i,j)=vSum/rho(i,j);
end
end
%每 saveStep 步保存下误差，检查收敛性
if mod(kk,saveStep)==0
temp=0;
for i=1:xLen
for j=1:yLen
temp=temp+sqrt(u(i,j)^2+v(i,j)^2); %节点的速度和
end
end
VSSum2=[VSSum2;temp];
u2=rot90(u);
v2=rot90(v);
for i=1:xLen
for j=1:yLen
speed(j,i)=sqrt(u2(j,i)^2+v2(j,i)^2);%节点的速度
end
end
set(pc,'cdata',speed);
drawnow
saveas(gcf,['D:\MyStudy\MATLAB\YuBrian','\',num2str(kk/saveStep),'.tif'])
picLength=kk/saveStep; %保存的图片数量
%每 checkStep 次保存下误差
if mod(kk,checkStep)==0
VSSum=[VSSum;temp];
%计算节点速度和变化（网格点总数*误差），判断是否跳出循环
if length(VSSum)>=2 && abs(VSSum(end)-VSSum(end-1))/VSSum(end-1)<=(1e-6)
fprintf('迭代步数:%d \n',kk);
fprintf('迭代误差:%d\n',abs(VSSum(end)-VSSum(end-1))/(xLen*yLen));
break;
end
end
end
end
runtime = toc(tstart);%仿真用的时间
%结果后处理
u2=rot90(u);
v2=rot90(v);
for i=1:xLen
for j=1:yLen
speed(j,i)=sqrt(u2(j,i)^2+v2(j,i)^2);
end
end
meshX2=(meshX-1)/(xLen-1);
meshY2=(meshY-1)/(xLen-1);
%速度场彩图 1
figure
pcolor(meshX2,meshY2,speed);
shading interp
colorbar
%axis equal tight
xlim([0 1])
ylim([0 (yLen-1)/(xLen-1)])
xlabel('x/L','fontsize',16)
ylabel('y/L','fontsize',16)
axis equal tight
set(gcf, 'Color', 'w');
set(gcf,'position',[440 447 560 351])
print(gcf,'-djpeg','-r600',[pwd '\Re' num2str(Re) '-11.jpg'])
%绘制流线图 2
figure
streamslice(meshX2,meshY2,u2,v2,4,'nearest ' )
axis equal tight
xlim([0 1])
ylim([0 (yLen-1)/(xLen-1)])
xlabel('x/L','fontsize',16)
ylabel('y/L','fontsize',16)
set(gcf, 'Color', 'w');
set(gcf,'position',[440 447 560 351])
print(gcf,'-djpeg','-r600',[pwd '\Re' num2str(Re) '-22.jpg'])