turbulence.zip_MRT-LBM_bgk_mrt_turbulence matlab_zip

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)