资源描述:
《Matlab实现玻尔兹曼晶格模拟》由会员上传分享,免费在线阅读,更多相关内容在应用文档-天天文库。
1、Matlab实现格子玻尔兹曼方法(LatticeBoltzmannMethod,LBM)模拟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;%Reynoldsnumb
2、ernu=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
3、,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%MACROSCO
4、PICVARIABLESrho=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
5、)=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
6、).*...(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
7、));%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)])
8、;end%VISUALIZATIONif(mod(cycle,tPlot)==0)u=reshape(sqrt(ux.^2+uy.^2),lx,ly);u(bbRegion)=nan;imagesc(u');axisequaloff;drawnowendend