clc
clear all
close all
x=100
y=200
tau=0.5*ones(x,y);
step=100
u0=0.05;
%% D2Q9 setup
ex(1,1,1:9)=[0 1 0 -1 0 1 -1 -1 1];
ey(1,1,1:9)=[0 0 1 0 -1 1 1 -1 -1];
c0=[0,0];c1=[1,0];c2=[0,1];c3=[-1,0];c4=[0,-1];c5=[1,1];c6=[-1,1];c7=[-1,-1];c8=[1,-1];
w(1,1,1:9)=[4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36];
%% setup solid boundary
[i,j]=meshgrid(1:y,1:x);
Solid=false(x,y);
Solid=Solid | (j-x/3).^2+(i-y/5).^2<(x/20)^2;
Solid=Solid | (j-x/2).^2+(i-y/2).^2<(x/20)^2;
Solid=Solid | (j-x/1.5).^2+(i-y/1.25).^2<(x/20)^2;
Solid(1,:)=true;
Solid(x,:)=true;
Fluid=~Solid;
Inlet=Fluid & i==1;
Outlet=i==y;
FluidAll=repmat(Fluid,[1,1,9]);
%% used for streaming indexing
Ind=reshape(1:x*y*9,x,y,9);
Ind(:,:,1+1)=circshift(Ind(:,:,1+1),c1);
Ind(:,:,2+1)=circshift(Ind(:,:,2+1),c2);
Ind(:,:,3+1)=circshift(Ind(:,:,3+1),c3);
Ind(:,:,4+1)=circshift(Ind(:,:,4+1),c4);
Ind(:,:,5+1)=circshift(Ind(:,:,5+1),c5);
Ind(:,:,6+1)=circshift(Ind(:,:,6+1),c6);
Ind(:,:,7+1)=circshift(Ind(:,:,7+1),c7);
Ind(:,:,8+1)=circshift(Ind(:,:,8+1),c8);
%% used for bounceback indexing
Blank=false(x,y);
SolidBB1=cat(3,Blank,Solid,Solid,Blank,Blank,Solid,Solid,Blank,Blank);
SolidBB2=cat(3,Blank,Blank,Blank,Solid,Solid,Blank,Blank,Solid,Solid);
temp=Ind(SolidBB1);Ind(SolidBB1)=Ind(SolidBB2);Ind(SolidBB2)=temp;
%% viscosity
Re=Inf;
Cs=5e-3
viscosity=u0*y/Re
tau0=3*viscosity+.5;
%% LES setup
exx(1,1,1:9)=[0 1 0 1 0 1 1 1 1];
eyy(1,1,1:9)=[0 0 1 0 1 1 1 1 1];
exy(1,1,1:9)=[0 0 0 0 0 1 -1 1 -1];
%% MRT setup
M=[ 1 -4 4 0 0 0 0 0 0
1 -1 -2 1 -2 0 0 1 0
1 -1 -2 0 0 1 -2 -1 0
1 -1 -2 -1 2 0 0 1 0
1 -1 -2 0 0 -1 2 -1 0
1 2 1 1 1 1 1 0 1
1 2 1 -1 -1 1 1 0 -1
1 2 1 -1 -1 -1 -1 0 1
1 2 1 1 1 -1 -1 0 -1];
Mt=M';
M0=diag(1./sum(M.*M))*Mt;
S0=[0,1.63,1.14,0,1.92,0,1.92];
% S0=[0,1.64,1.54,0,1.7,0,1.7];
S0=ones(x*y,1)*S0;
%% UI control
restart = uicontrol('position',[20 20 80 20],'style','toggle','string','restart');
quit = uicontrol('position',[120 20 80 20],'style','toggle','string','close');
pausebutton = uicontrol('position',[220 20 80 20],'style','toggle','string','pause');
bgkmrt = uicontrol('position',[320 20 80 20],'style','toggle','string','switch to bgk');
%%
while get(quit,'value') == 0
set(restart,'value',0);
set(pausebutton,'value',0);
%% initial condition
nTimestep=0;
ux=zeros(x,y);
uy=u0*ones(x,y);
forcex=zeros(x,y);
forcey(1:x,1:y)=2e-5;
rho=ones(x,y);
usq=ux.*ux+uy.*uy;
c=(bsxfun(@times,ex,ux)+bsxfun(@times,ey,uy));
density=bsxfun(@times,rho,w) .* bsxfun(@plus,1+3*c+4.5*c.*c,-1.5*usq);
%% loop
while get(restart,'value')==0 && get(quit,'value')==0
while get(pausebutton,'value')==1 && get(restart,'value')==0 && get(quit,'value')==0
set(pausebutton,'string','continue');pause(.01);end;set(pausebutton,'string','pause');%% restart and quit are volatile, so extra test is required
%% calc macro
rho=sum(density,3);
ux=sum(bsxfun(@times,density,ex),3)./rho;
uy=sum(bsxfun(@times,density,ey),3)./rho;
%% forcing
ux=ux+forcex./rho;
uy=uy+forcey./rho;
%% collision
usq=ux.*ux+uy.*uy;
c=(bsxfun(@times,ex,ux)+bsxfun(@times,ey,uy));
density_eq=bsxfun(@times,rho,w) .* bsxfun(@plus,1+3*c+4.5*c.*c,-1.5*usq);
density_delta=density_eq-density;
moment=reshape(density,[x*y,9])*M; ux_=ux(:);uy_=uy(:);usq_=usq(:);
moment_eq=bsxfun(@times,rho(:),[ones(x*y,1),-2+3*usq_,1-3*usq_,ux_,-ux_,uy_,-uy_,ux_.*ux_-uy_.*uy_,ux_.*uy_]);
%% eddy viscosity with Smagorinsky model
qxx=sum(bsxfun(@times,density_delta,exx),3);
qyy=sum(bsxfun(@times,density_delta,eyy),3);
qxy=sum(bsxfun(@times,density_delta,exy),3);
q=sqrt(qxx.*qxx+qyy.*qyy+2*qxy.*qxy);
tau=tau0+.5*(sqrt(tau0.*tau0+36*Cs*q)-tau0);
%% collision cont'd
if(get(bgkmrt,'value')==0)
set(bgkmrt,'string','switch to bgk')
S=[S0,1./tau(:),1./tau(:)];
else
set(bgkmrt,'string','switch to mrt')
S=1./tau(:)*ones(1,9);
end
moment=moment+bsxfun(@times,moment_eq-moment,S);
density_temp=reshape(moment*M0,[x,y,9]);
% density_temp=density+bsxfun(@rdivide,density_delta,tau);
density(FluidAll)=density_temp(FluidAll);
%% stream & bounce back
density=density(Ind);
%%
nTimestep=nTimestep+1;
if(~mod(nTimestep,step))
uxplot=uy.*Fluid;
uyplot=ux.*Fluid;
umag=sqrt(uxplot.*uxplot+uyplot.*uyplot);
vmax=max(umag(:));
vel=(256/vmax)*(umag);
vor=curl(uxplot,uyplot)*20000+128;
subplot(2,1,1)
pcolor(vel)
xlabel('Velocity')
shading interp
axis equal
axis([1 y 1 x])
title(sprintf('%d, tau=%f, u_m_a_x=%f, Vorticity_m_a_x=%f, Mass=%f',nTimestep,tau(x/2,y/5),vmax,max(vor(:)),sum(rho(:))));
subplot(2,1,2)
pcolor(vor)
xlabel('Vorticity')
shading interp
axis equal
axis([1 y 1 x])
drawnow
end
end
end
close(gcf)
