%%主程序%%%%%%%%%%%%%%%%%%%%%%%
clc; clear; close;
lx=1680; ly=140;
lyp=0.005;%m
Lr=lyp/ly;
cs=1/3^0.5; c_squ=cs^2;
t1=4/9; t2=1/9; t3=1/36;
Re=400;
Uin=0.02; D=fix(0.022/Lr);
nu=D*Uin/Re;
tau=3*nu+0.5;
omega=1/tau;
t= [1/9,1/36,1/9,1/36,1/9,1/36,1/9,1/36,4/9];
ex=[1, 1, 0, -1, -1, -1, 0, 1, 0];
ey=[0, 1, 1, 1, 0, -1, -1, -1, 0];
oppo=[5,6,7,8,1,2,3,4,9];
rho=ones(ly,lx);
UX=zeros(ly,lx);
UX=Uin;
UY=zeros(ly,lx);
F=zeros(ly,lx,9);
FEQ=FEQf(UX,UY,rho,t1,t2,t3,c_squ);
F=FEQ; Fb=F;
[y,x]=meshgrid(1:ly,1:lx); %生成序号网格
bound =[((x’-[lx/4]).^2 + (y’-[D/2]+510).^2 )<=(D/2)^2]+…
[y’==1]+[y’==ly];
obs=find(bound); %找到壁面网格的全局位置
m=length(obs);
b2b=zeros(9*m,1); bid=zeros(9*m,1);
for i=0:8
bid(i*m+1:(i+1)*m)=(oppo(i+1)-1)*lx*ly+obs;
b2b(i*m+1:(i+1)*m)=i*lx*ly+obs;
end
for i=1:80000
F=stream(F,ly,lx);
rho = sum(F,3); %计算密度
UX=(sum(F(:,:,[1 2 8]),3)-sum(F(:,:,[4 5 6]),3))./rho; %计算x速度
UY=(sum(F(:,:,[2 3 4]),3)-sum(F(:,:,[6 7 8]),3))./rho; %计算y速度
UX(obs)=0; UY(obs)=0;
UX(2:ly-1,1)=Uin*(2:ly-1)/(ly-1).*(1-(2:ly-1)/(ly-1))*4; UY(2:ly-1,1)=0;
rho (2:ly-1,1)=rho (2:ly-1,2); %密度采用外推
UX(2:ly-1,lx)= UX(2:ly-1,lx-1); %右端出口，速度采用外推
UY(2:ly-1,lx)= UY(2:ly-1,lx-1);
rho (2:ly-1,lx)=1; %密度指定
P=cs^2*rho;
FEQ=FEQf(UX,UY,rho,t1,t2,t3,c_squ);
BOUNCEBACK=F(b2b); %获得反弹量
F=omega*FEQ+(1-omega)*F; %碰撞
F(bid)=BOUNCEBACK; %将反弹量赋予相应位置
F(2:ly-1,1,:)=FEQ(2:ly-1,1,:)+((F(2:ly-1,2,:)-FEQ(2:ly-1,2,:))); %左边入口
F(2:ly-1,lx,:)=FEQ(2:ly-1,lx,:)+((F(2:ly-1,lx-1,:)-FEQ(2:ly-1,lx-1,:))); %右边出口
%结果显示
OutputRes;
end
%%平衡态分布函数%%%%%%%%%%%%%%%%%%%%%%%
function FEQ=FEQf(UX,UY,rho,t1,t2,t3,c_squ)
U_SQU=UX.^2+UY.^2; U_C2=UX+UY;
U_C4=-UX+UY; U_C6=-U_C2;
U_C8=-U_C4;
FEQ(:,:,9)=t1*rho.*(1-U_SQU/(2*c_squ));
FEQ(:,:,1)=t2*rho.*(1+UX/c_squ+0.5*(UX/c_squ).^2-U_SQU/(2*c_squ));
FEQ(:,:,3)=t2*rho.*(1+UY/c_squ+0.5*(UY/c_squ).^2-U_SQU/(2*c_squ));
FEQ(:,:,5)=t2*rho.*(1-UX/c_squ+0.5*(UX/c_squ).^2-U_SQU/(2*c_squ));
FEQ(:,:,7)=t2*rho.*(1-UY/c_squ+0.5*(UY/c_squ).^2-U_SQU/(2*c_squ));
FEQ(:,:,2)=t3*rho.*(1+U_C2/c_squ+0.5*(U_C2/c_squ).^2-U_SQU/(2*c_squ));
FEQ(:,:,4)=t3*rho.*(1+U_C4/c_squ+0.5*(U_C4/c_squ).^2-U_SQU/(2*c_squ));
FEQ(:,:,6)=t3*rho.*(1+U_C6/c_squ+0.5*(U_C6/c_squ).^2-U_SQU/(2*c_squ));
FEQ(:,:,8)=t3*rho.*(1+U_C8/c_squ+0.5*(U_C8/c_squ).^2-U_SQU/(2*c_squ));
%%迁移%%%%%%%%%%%%%%%%%%%%%%%
function F=stream(F,ly,lx)
F(:,:,1)=F(:,[lx 1:lx-1],1);
F(:,:,2)=F([ly 1:ly-1], [lx 1:lx-1],2);
F(:,:,3)=F([ly 1:ly-1],:,3);
F(:,:,4)=F( [ly 1:ly-1], [2:lx 1],4);
F(:,:,5)=F(:,[2:lx 1],5);
F(:,:,6)=F([2:ly 1], [2:lx 1],6);
F(:,:,7)=F([2:ly 1],:,7);
F(:,:,8)=F([2:ly 1],[lx 1:lx-1],8);
%%结果显示%%%%%%%%%%%%%%%%%%%%%%%
if mod(i,300)==0
clc ; clf; i
U = sqrt(UX.^2+UY.^2);
hold on
DUX=diff(UX);
DUY=diff(UY);
Vor=DUY(1:ly-1,1:lx-1)-DUX(1:ly-1,1:lx-1);
contour(Vor,100);
imagesc(0,300,U,[0,1.5*Uin])
XL=(2:2:ly-1)+150;
streamline(x’,y’+150,UX,UY,ones(1,length(XL)),XL,1)
colormap(‘jet’);axis equal; drawnow;
hold off
end