Matlab实现格子玻尔兹曼方法(LatticeBoltzmannMethod,LBM)模拟

Matlab实现格子玻尔兹曼方法(LatticeBoltzmannMethod,LBM)模拟%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)42+(y-obst_y)42fIn(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.A2+uy42);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(ux42+uy.A2),lx,ly);u(bbRegion)=nan;imagesc(u);axisequaloff;drawnowendend