%红色代码为增加代码；绿底为屏蔽的代码
%初始化
clc; clear; close;
lx=25; ly=9;
Re=1.25;
yu=0.1;
fx0=Re*8*yu^2/(ly-1)^3;
W=[4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36];
ex=[0, 1, 0, -1, 0 1,-1, -1, 1];
ey=[0, 0, 1, 0, -1, 1, 1, -1,-1];
cs=1/sqrt(3);
opposite=[1,4,5,2,3,8,9,6,7];
rho=ones(lx,ly);
u=zeros(lx,ly);
v=zeros(lx,ly);
fx=zeros(lx,ly)+fx0;
fy=zeros(lx,ly);
M1=W(1)*rho; M2=W(2)*rho; M3=W(6)*rho;
F= cat(3,M1,M2,M2,M2,M2,M3,M3,M3,M3);
Feq=F;
%迁移下标
D2x=[lx,1:lx-1]; D2y=[1:ly];
D3x=[1:lx]; D3y=[ly,1:ly-1];
D4x=[2:lx,1]; D4y=[1:ly];
D5x=[1:lx]; D5y=[2:ly,1];
D6x=[lx,1:lx-1]; D6y=[ly,1:ly-1];
D7x=[2:lx,1]; D7y=[ly,1:ly-1];
D8x=[2:lx,1]; D8y=[2:ly,1];
D9x=[lx,1:lx-1]; D9y=[2:ly,1];
%% 主循环===========================================
for tStep=1:20000
% 迁移
F0(:,:,1)=F(:,:,1);
F0(:,:,2)=F(D2x,D2y,2);
F0(:,:,3)=F(D3x,D3y,3);
F0(:,:,4)=F(D4x,D4y,4);
F0(:,:,5)=F(D5x,D5y,5);
F0(:,:,6)=F(D6x,D6y,6);
F0(:,:,7)=F(D7x,D7y,7);
F0(:,:,8)=F(D8x,D8y,8);
F0(:,:,9)=F(D9x,D9y,9);
F=F0;
%计算宏观量
rho=sum(F,3);
u=(sum(F(:,:,[2,6,9]),3)-sum(F(:,:,[4,7,8]),3))./rho+fx/2./rho;
v=(sum(F(:,:,[3,6,7]),3)-sum(F(:,:,[5,8,9]),3))./rho+fy/2./rho;
%设置宏观边界
rho(:,1)=rho(:,2);
rho(:,ly)=rho(:,ly-1);
u(:,1)=0; u(:,ly)=0;
v(:,1)=0; v(:,ly)=0;
%计算平衡态函数
for i=1:9
M1=(ex(i)*u+ey(i)*v)/cs^2;
M2=M1.^2/2;
M3=(u.^2+v.^2)/2/cs^2;
Feq(:,:,i)=rho.*W(i).*(1+M1+M2-M3);
end
%碰撞
for i=1:9
M1=3*(ex(i)-u).*fx+(ey(i)-v).*fy;
M2=9*(ex(i)*u+ey(i)*v).*(ex(i).*fx+ey(i).*fy);
force(:,:,i)=W(i)*(1-0.5/(0.5+3*yu))*(M1+M2);
end
F=F-1/(0.5+3*yu)*(F-Feq)+force;
%微观边界设置
F(:,[1,ly],:)=F(:,[1,ly],opposite);
% F(:,1,:)=Feq(:,1,:)+(F(:,2,:)-Feq(:,2,:));
% F(:,ly,:)=Feq(:,ly,:)+(F(:,ly-1,:)-Feq(:,ly-1,:));
%后处理
if mod(tStep,500)==0
clc;clf
tStep
U=(u.^2+v.^2).^0.5;
subplot(1,2,1),imagesc(U’);
colormap jet
axis equal
cy=(ly-1)/2;
y0=fx0/2/yu*(cy.^2-([0:ly-1]-cy).^2);
Ni=floor(lx/2);
subplot(1,2,2),hold on,plot(y0,’k’),plot(u(Ni,:),’r’); hold off
drawnow
end
end % 主循环结束============================