[passster password=”IB-LBM” area=”3212″]
%红色代码为增加代码；绿底为屏蔽的代码
%初始化
clc; clear; close;
lx=400; ly=81;
%构建壁面
Bound=zeros(lx,ly);
Bound(:,[1,ly])=1;
%构建圆柱边界
xc=lx/4; yc=ly/2; Rc=10;
Nc=floor(2*2*pi*Rc);
ceta=2*pi/Nc:2*pi/Nc:2*pi;
Xcy=xc+Rc*cos(ceta); Ycy=yc+Rc*sin(ceta);
ds=0.5;
Fx=zeros(1,Nc); Fy=zeros(1,Nc);
Uxi=zeros(1,Nc); Uyi=zeros(1,Nc);
Cd=[]; Cl=[];
% 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
for tStep=1:25000
% 迁移
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)=1;%rho(:,2); %下壁面
u(:,1)=0; v(:,1)=0; %
rho(:,ly)=1;%rho(:,ly-1); %上壁面
u(:,ly)=0; v(:,ly)=0; %
% u(Bin)=0; v(Bin)=0; %所有固体
rho(1,:)=rho(2,:); %入口
u(1,:)=Uin; v(1,:)=0; %
rho(lx,:)=1; %出口
u(lx,:)=u(lx-1,:); v(lx,:)=v(lx-1,:);
% 获得浸入边界上的速度
Ux=zeros(1,Nc); Uy=zeros(1,Nc);
for i=1:Nc
for j=floor(Xcy(i))-1:floor(Xcy(i))+2
for k=floor(Ycy(i))-1:floor(Ycy(i))+2
rx=j-Xcy(i); ry=k-Ycy(i);
Dx=1/8*(3-2*abs(rx)+sqrt(1+4*abs(rx)-4*rx^2))*(abs(rx)<=1);
Dx=Dx+1/8*(5-2*abs(rx)-sqrt(-7+12*abs(rx)-4*rx^2))*(abs(rx)<=2)*(abs(rx)>1);
Dy=1/8*(3-2*abs(ry)+sqrt(1+4*abs(ry)-4*ry^2))*(abs(ry)<=1);
Dy=Dy+1/8*(5-2*abs(ry)-sqrt(-7+12*abs(ry)-4*ry^2))*(abs(ry)<=2)*(abs(ry)>1);
Ux(i)=Ux(i)+u(j,k)*Dx*Dy; Uy(i)=Uy(i)+v(j,k)*Dx*Dy;
end
end
end
%计算总体积力
Uxi=Uxi+Ux; Uyi=Uyi+Uy;
Fx=-0.02*Uxi-2*Ux; Fy=-0.02*Uyi-2*Uy;
%浸入边界计算，将L点n时刻体积力spread到n+1时刻的Eular点上去。
fx=zeros(lx,ly); fy=zeros(lx,ly);
%清空体积力存放空间
for i=1:Ncfor i=1:Nc
for j=floor(Xcy(i))-1:floor(Xcy(i))+2
for k=floor(Ycy(i))-1:floor(Ycy(i))+2
rx=j-Xcy(i); ry=k-Ycy(i);
Dx=1/8*(3-2*abs(rx)+sqrt(1+4*abs(rx)-4*rx^2))*(abs(rx)<=1);
Dx=Dx+1/8*(5-2*abs(rx)-sqrt(-7+12*abs(rx)-4*rx^2))*(abs(rx)<=2)*(abs(rx)>1);
Dy=1/8*(3-2*abs(ry)+sqrt(1+4*abs(ry)-4*ry^2))*(abs(ry)<=1);
Dy=Dy+1/8*(5-2*abs(ry)-sqrt(-7+12*abs(ry)-4*ry^2))*(abs(ry)<=2)*(abs(ry)>1);
fx(j,k)=fx(j,k)+Fx(i)*ds*Dx*Dy; fy(j,k)=fy(j,k)+Fy(i)*ds*Dx*Dy;
end
end
end
%计算平衡态函数
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,:)=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(2,2,1),imagesc(U’);
colormap jet
box on
axis equal
cy=(ly-1)/2;
y0=fx0/2/yu*(cy.^2-([0:ly-1]-cy).^2);
Ni=floor(lx/2);
subplot(2,2,2),hold on,plot(y0,’k’),plot(u(Ni,:),’r’); hold off
if tStep>1000
Cd=[Cd,-sum(Fx*ds)/(0.5*Uin^2*2*Rc)]; Cl=[Cl,-sum(Fy*ds)/(0.5*Uin^2*2*Rc)];
end
subplot(2,2,3),plot(Cd,’b’)
subplot(2,2,4),plot(Cl,’b’)
drawnow
end
end % 主循环结束============================