Main.m
clc; clear; close;
colormap jet
lx=300; ly=60;
cs=1/3^0.5; c_squ=cs^2;
t1=4/9; t2=1/9; t3=1/36;
Re=100;
Uin=0.05;
nu=(ly-1)*Uin/Re;
tau=3*nu+0.5;
omega=1/tau;
oppo=[5,6,7,8,1,2,3,4,9];
t= [1/9,1/36,1/9,1/36,1/9,1/36,1/9,1/36,4/9];%D2Q9权重
ex=[1, 1, 0, -1, -1, -1, 0, 1, 0]; % D2Q9方向向量
ey=[0, 1, 1, 1, 0, -1, -1, -1, 0]; % D2Q9方向向量
rho=ones(ly,lx);
UX=zeros(ly,lx);
UX(:,1:lx/2)=0.1;
UY=zeros(ly,lx);
F=zeros(ly,lx,9);
F=FEQf(UX,UY,rho,t1,t2,t3,c_squ);
for i=1:10000
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速度
P=cs^2*rho;
UX([1,ly],:)=0;UY([1,ly],:)=0; %壁面
UX(:,1)=Uin; UY(:,1)=0; %左入口
rho(:,1)=rho(:,2);
UX(:,lx)=UX(:,lx-1); UY(:,lx)=0; %右出口
rho(:,1)=1;
FEQ=FEQf(UX,UY,rho,t1,t2,t3,c_squ);
Fb=F;
F=omega*FEQ+(1-omega)*F; %碰撞
F([1,ly],:,:)=Fb([1,ly],:,oppo);
F(:,1,:)=FEQ(:,1,:)+(F(:,2,:)-FEQ(:,2,:));
F(:,lx,:)=FEQ(:,lx,:)+(F(:,lx-1,:)-FEQ(:,lx-1,:));
if mod(i,20)==0
U= (UX.^2+ UY.^2).^0.5;
imagesc(U);
axis equal
drawnow
axis xy
end
end
stream.m
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);
FEQf.m
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(:,:,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));
FEQ(:,:,9)=t1*rho.*(1-U_SQU/(2*c_squ));