Matlab实现格子玻尔兹曼方法

Matlab实现格子玻尔兹曼方法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%cylinder.m:Flowaroundacyliner,usingLBM%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Thisprogramisfreesoftware;youcanredistributeitand/or%modifyitunderthetermsoftheGNUGeneralPublicLicense%aspublishedbytheFreeSoftwareFoundation;eitherversion2%oftheLicense,or(atyouroption)anylaterversion.%Thisprogramisdistributedinthehopethatitwillbeuseful,%butWITHOUTANYWARRANTY;withouteventheimpliedwarrantyof%MERCHANTABILITYorFITNESSFORAPARTICULARPURPOSE.Seethe%GNUGeneralPublicLicenseformoredetails.%YoushouldhavereceivedacopyoftheGNUGeneralPublic%Licensealongwiththisprogram;ifnot,writetotheFree%SoftwareFoundation,Inc.,51FranklinStreet,FifthFloor,%Boston,MA02110-1301,USA.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear%GENERALFLOWCONSTANTSlx=250;ly=51;obst_x=lx/5+1;%positionofthecylinder;(exactobst_y=ly/2+1;%y-symmetryisavoided)obst_r=ly/10+1;%radiusofthecylinderuMax=0.02;%maximumvelocityofPoiseuilleinflowRe=100;%Reynoldsnumbernu=uMax*2.*obst_r/Re;%kinematicviscosityomega=1./(3*nu+1./2.);%relaxationparametermaxT=400000;%totalnumberofiterationstPlot=5;%cycles%D2Q9LATTICECONSTANTSt=[4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36];cx=[0,1,0,-1,0,1,-1,-1,1];cy=[0,0,1,0,-1,1,1,-1,-1];opp=[1,4,5,2,3,8,9,6,7];col=[2:(ly-1)];[y,x]=meshgrid(1:ly,1:lx);obst=(x-obst_x).^2+(y-obst_y).^2=obst_r.^2;obst(:,[1,ly])=1;bbRegion=find(obst);%INITIALCONDITION:(rho=0,u=0)==fIn(i)=t(i)fIn=reshape(t'*ones(1,lx*ly),9,lx,ly);%MAINLOOP(TIMECYCLES)forcycle=1:maxT%MACROSCOPICVARIABLESrho=sum(fIn);ux=reshape(...(cx*reshape(fIn,9,lx*ly)),1,lx,ly)./rho;uy=reshape(...(cy*reshape(fIn,9,lx*ly)),1,lx,ly)./rho;%MACROSCOPIC(DIRICHLET)BOUNDARYCONDITIONS%Inlet:PoiseuilleprofileL=ly-2;y=col-1.5;ux(:,1,col)=4*uMax/(L*L)*(y.*L-y.*y);uy(:,1,col)=0;rho(:,1,col)=1./(1-ux(:,1,col)).*(...sum(fIn([1,3,5],1,col))+...2*sum(fIn([4,7,8],1,col)));%Outlet:Zerogradientonrho/uxrho(:,lx,col)=rho(:,lx-1,col);uy(:,lx,col)=0;ux(:,lx,col)=ux(:,lx-1,col);%COLLISIONSTEPfori=1:9cu=3*(cx(i)*ux+cy(i)*uy);fEq(i,:,:)=rho.*t(i).*...(1+cu+1/2*(cu.*cu)...-3/2*(ux.^2+uy.^2));fOut(i,:,:)=fIn(i,:,:)-...omega.*(fIn(i,:,:)-fEq(i,:,:));end%MICROSCOPICBOUNDARYCONDITIONSfori=1:9%LeftboundaryfOut(i,1,col)=fEq(i,1,col)+...18*t(i)*cx(i)*cy(i)*(fIn(8,1,col)-...fIn(7,1,col)-fEq(8,1,col)+fEq(7,1,col));%RightboundaryfOut(i,lx,col)=fEq(i,lx,col)+...18*t(i)*cx(i)*cy(i)*(fIn(6,lx,col)-...fIn(9,lx,col)-fEq(6,lx,col)+fEq(9,lx,col));%BouncebackregionfOut(i,bbRegion)=fIn(opp(i),bbRegion);end%STREAMINGSTEPfori=1:9fIn(i,:,:)=...circshift(fOut(i,:,:),[0,cx(i),cy(i)]);end%VISUALIZATIONif(mod(cycle,tPlot)==0)u=reshape(sqrt(ux.^2+uy.^2),lx,ly);u(bbRegion)=nan;imagesc(u');axisequaloff;drawnowendend