%红色代码为增加代码；绿底为屏蔽的代码
%初始化
clc; clear; close;
% lx=25; ly=9;
% Re=1.25;
% yu=0.1;
% fx0=Re*8*yu^2/(ly-1)^3;
lx=400; ly=81;
%构建壁面及圆柱区域
Bound=zeros(lx,ly);
Bound(:,[1,ly])=1;
xc=lx/4; yc=ly/2; Rc=6;
for i=floor(xc)-7:floor(xc)+7
for j=floor(yc)-7:floor(yc)+7
if ((i-xc)^2+(j-yc)^2)^0.5<=Rc
Bound(i,j)=1;
end
end
end
Bin=[]; BFin0=[]; BFin=[];
opposite=[1,4,5,2,3,8,9,6,7];
for i=1:lx
for j=1:ly
if Bound(i,j)==1
Bin=[Bin,i+lx*(j-1)];
for k=1:9
BFin=[BFin,(i+lx*(j-1))+(opposite(k)-1)*lx*ly];
BFin0=[BFin0,(i+lx*(j-1))+(k-1)*lx*ly];
end
end
end
end
Re=100;
Uin=0.1;
yu=2*Rc*Uin/Re;
fx0=0;
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);
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:15000
% 迁移
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(Bin)=0; v(Bin)=0; %所有固体
rho(1,:)=rho(2,:); %入口
u(1,:)=Uin; v(1,:)=0; %
rho(lx,:)=rho(lx-1,:); %出口
u(lx,:)=u(lx-1,:); v(lx,:)=v(lx-1,:);
%计算平衡态函数
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(BFin0)=F(BFin); %半步反弹操作
F(1,:,:)=Feq(1,:,:)+(F(2,:,:)-Feq(2,:,:));
F(lx,:,:)=Feq(lx,:,:)+(F(lx-1,:,:)-Feq(lx-1,:,:));
% 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,100)==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 % 主循环结束============================