格子玻尔兹曼MATLAB运用(LBGK-D2Q9-poiseuille-channel2D)

%GianniSchenaJuly2005,[e-mail]%LatticeBoltzmannLBE,geometry:D2Q9,model:BGK%ApplicationtopermeabilityinporousmediaRestart=false%torestartfromanearlierconvergencelogical(Restart);ifRestart==false;closeall,clearall%startfromscratchandclean...Restart=false;%typeofchannelgeometry;%oneoftheflollowingflags==truePois_test=true,%noobstaclesinthe2Dchannel%poroussystemsobs_regolare=false%obs_irregolare=false%tic%IN%|vvvv|+y%|vvvv|^%|vvvv||-+x%OUT%Poresin2D:WetandDrylocations(Wet==1,Dry==0)wXh_Dry=[3,1];wXh_Wet=[3,4];ifobs_regolare,%withinternalobstaclesA=repmat([zeros(wXh_Dry),ones(wXh_Wet)],[1,3]);A=[A,zeros(wXh_Dry)];B=ones(size(A));C=[A;B];D=repmat(C,4,1);D=[B;D]endifobs_irregolare,%withintobstaclesA1=repmat([zeros(wXh_Dry),ones(wXh_Wet)],[1,3]);A1=[A1,zeros(wXh_Dry)];B=ones(size(A1));C1=repmat([ones(wXh_Wet),zeros(wXh_Dry)],[1,3]);C1=[C1,ones(wXh_Dry)];E=[A1;B;C1;B];D=repmat(E,2,1);D=[B;D]endif~Pois_testfigure,imshow(D,[])Channel2D=D;Len_Channel_2D=size(Channel2D,1);%LengthWidth=size(Channel2D,2);%shouldnotbehodChannel_2D_half_Width=Width/2,end%testwithoutobstacles(i.e.2Dchannel&noobstacles)ifPois_test%over-writesthedefinitionoftheporespaceclearChannel2DLen_Channel_2D=36,%lunghezzacanale2dChannel_2D_half_Width=8;Width=Channel_2D_half_Width*2;Channel2D=ones(Len_Channel_2D,Width);%definewetarea%Channel2D(6:12,6:8)=0;%putfluidobstacleimshow(Channel2D,[]);end[NrMc]=size(Channel2D);%NumberrowsandMunbercolumns%porosityporosity=nnz(Channel2D==1)/(Nr*Mc)%FLUIDPROPERTIES%physicalpropertiescs2=1/3;%cP_visco=0.5;%[cP]1CPDinamicwaterviscosity20Cdensity=1.;%fluiddensityLky_visco=cP_visco/density;%latticekinematicviscosityomega=(Lky_visco/cs2+0.5).^-1;%omega:relaxationfrequency%Lky_visco=cs2*(1/omega-0.5),%latticekinematicviscosity%dPdL=Pressure/dL;%Externalpressuregradient[atm/cm]uy_fin_max=-0.2;%dPdL=abs(2*Lky_visco*uy_fin_max/(Channel_2D_half_Width.^2));dPdL=-0.0125;uy_fin_max=dPdL*(Channel_2D_half_Width.^2)/(2*Lky_visco);%PoiseuilleGradient;%maxpoiseuillefinalvelocityontheflowprofileuy0=-0.001;ux0=0.0001;%linearvel..inizialization%%uy_fin_max=-0.2;%maxpoiseuillefinalvelocityontheflowprofile%omega=0.5,cs2=1/3;%omega:relaxationfrequency%Lky_visco=cs2*(1/omega-0.5),%latticekinematicviscosity%dPdL=abs(2*Lky_visco*uy_fin_max/(Channel_2D_half_Width.^2));%PoiseuilleGradient;%uyf_av=uy_fin_max*(2/3);;%averagefluidvelocityontheprofilex_profile=([-Channel_2D_half_Width:+Channel_2D_half_Width-1]+0.5);uy_analy_profile=uy_fin_max.*(1-(x_profile/Channel_2D_half_Width).^2);%analyticalvelocityprofileav_vel_t=1.e+10;%inizialization(t=0)%PixelSize=5;%[Microns]%dL=(Nr*PixelSize*1.0E-4);%samplehight[cm]%%EXPERIMENTALSET-UP%inletandoutletbuffersinb=2,oub=2;%inletandoutletbuffersthickness%addfluidattheinlet(top)andoutlet(down)inlet=ones(inb,Mc);outlet=ones(oub,Mc);Channel2D=[[inlet];Channel2D;[outlet]];%addfluxinanddown(EtoW)[NrMc]=size(Channel2D);%updatesize%boundariesrelatedtotheexperimentalsetupwb=2;%wallthicknessChannel2D=[zeros(Nr,wb),Channel2D,zeros(Nr,wb)];%addwalls(nofluidleak)[NrMc]=size(Channel2D);%updatesizeuy_analy_profile=[zeros(1,wb),uy_analy_profile,zeros(1,wb)];%takeintoaccountwallsx_pro_fig=[[x_profile(1)-[wb:-1:1]],[x_profile,[1:wb]+x_profile(end)]];%Figureplotsanalyticalparabolicprofilefigure(20),plot(x_pro_fig,uy_analy_profile,'-'),gridon,title('Analyticalparab.profileforPoiseuilleplanarflowinachannel')%VISUALIZEPORESPACE&FLUIDOSTACLES&MEDIALAXISfigure,imshow(Channel2D);title('Vasselgeometry');Channel2D=logical(Channel2D);%obstaclesforBounceBack(infrontofthegrain)Obstacles=bwperim(Channel2D,8);%perimeterofthegrainsforbouncebackBound.Cond.border=logical(ones(Nr,Mc));border([1:inb,Nr-oub:Nr],[wb+2:Mc-wb-1])=0;Obstacles=Obstacles.*(border);figure,imshow(Obstacles);title('Fluidobstacles(inthefluid)');%Medial_axis=bwmorph(Channel2D,'thin',Inf);%figure,imshow(Medial_axis);title('Medialaxis');figure(10)%usedtovisualizeevolutionofrhofigure(11)%usedtovisualizeuxfigure(12)%usedtovisualizeuy(i.e.top-down)%INDICES%Wetlocationsetc.[iabw1jabw1]=find(Channel2D==1);%indicesi,j,ofactivelatticelocationsi.e.porelena=length(iabw1);%numberofactivelocationi.e.ofporespacelatticecellsija=(jabw1-1)*Nr+iabw1;%equivalentsingleindex(i,j)-ijaforactivelocations%absolute(singleindex)positionoftheobstaclesinforbouncebackinChannel2D%Obstacles[iobsjobs]=find(Obstacles);lenobs=length(iobs);ijobs=(jobs-1)*Nr+iobs;%asabove%Medialaxisoftheporespace[imajma]=find(Medial_axis);lenma=length(ima);ijma=(jma-1)*Nr+ima;%asabove%Internalwetlocations:wet&~obstables%(i.e.internalwetlatticelocationnonincontactwithdraylocations)[iawintjawint]=find((Channel2D==1&~Obstacles));%indicesi,j,ofactivelatticelocationslenwint=length(iawint);%numberofinternal(i.e.notborder)wetlocationsijaint=(jawint-1)*Nr+iawint;%equivalentsinglNxM=Nr*Mc;%DIRECTIONS:ENWSNENWSWSEZERO(ZERO:RestParticle)%y^%625^NWNNE%391...+x-+yWRPE%748SWSSE%-y%x&ycomponentsofvelocities,+xistoest,+yistonordEast=1;North=2;West=3;South=4;NE=5;NW=6;SW=7;SE=8;RP=9;N_c=9;%numberofdirections%versorsD2Q9C_x=[10-101-1