LBM顶盖流(matlab)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%cavity2d.m:2Dcavityflow,simulatedbyaLBmethod%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%LatticeBoltzmannsample,Matlabscript%Copyright(C)2006-2008JonasLatt%Address:RueGeneralDufour24,1211Geneva4,Switzerland%E-mail:[e-mail]%%Implementationof2dcavitygeometryandZou/Heboundary%conditionbyAdrianoSciacovelli%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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=128;ly=128;uLid=0.05;%horizontallidvelocityvLid=0;%verticallidvelocityRe=100;%Reynoldsnumbernu=uLid*lx/Re;%kinematicviscosityomega=1./(3*nu+1./2.);%relaxationparametermaxT=40000;%totalnumberofiterationstPlot=10;%cyclesforgraphicaloutput%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];lid=[2:(lx-1)];[y,x]=meshgrid(1:ly,1:lx);obst=ones(lx,ly);obst(lid,2:ly)=0;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)BOUNDARYCONDITIONSux(:,lid,ly)=uLid;%lidx-velocityuy(:,lid,ly)=vLid;%lidy-velocityrho(:,lid,ly)=1./(1+uy(:,lid,ly)).*(...sum(fIn([1,2,4],lid,ly))+2*sum(fIn([3,6,7],lid,ly)));%MICROSCOPICBOUNDARYCONDITIONS:LID(Zou/HeBC)fIn(5,lid,ly)=fIn(3,lid,ly)-2/3*rho(:,lid,ly).*uy(:,lid,ly);fIn(9,lid,ly)=fIn(7,lid,ly)+1/2*(fIn(4,lid,ly)-fIn(2,lid,ly))+...1/2*rho(:,lid,ly).*ux(:,lid,ly)-1/6*rho(:,lid,ly).*uy(:,lid,ly);fIn(8,lid,ly)=fIn(6,lid,ly)+1/2*(fIn(2,lid,ly)-fIn(4,lid,ly))-...1/2*rho(:,lid,ly).*ux(:,lid,ly)-1/6*rho(:,lid,ly).*uy(:,lid,ly);%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%MICROSCOPICBOUNDARYCONDITIONS:NO-SLIPWALLS(bounce-back)fori=1:9fOut(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(:,ly:-1:1)'./uLid);colorbaraxisequaloff;drawnowendend