Gianni SchenaJuly.docx

Gianni SchenaJuly.docx

《Gianni SchenaJuly.docx》由会员分享，可在线阅读，更多相关《Gianni SchenaJuly.docx（15页珍藏版）》请在冰豆网上搜索。

GianniSchenaJuly

%GianniSchenaJuly2005,schena@units.it

%LatticeBoltzmannLBE,geometry:

D2Q9,model:

BGK

%Applicationtopermeabilityinporousmedia

Restart=false%torestartfromanearlierconvergence

logical（Restart）;

ifRestart==false;

closeall,clearall%startfromscratchandclean...

Restart=false;

%typeofchannelgeometry;

%oneoftheflollowingflags==true

Pois_test=true,%noobstaclesinthe2Dchannel

%poroussystems

obs_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,%withinternalobstacles

A=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]

end

ifobs_irregolare,%withintobstacles

A1=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]

end

if~Pois_test

figure,imshow（D,[]）

Channel2D=D;

Len_Channel_2D=size（Channel2D,1）;%Length

Width=size（Channel2D,2）;%shouldnotbehod

Channel_2D_half_Width=Width/2,

end

%testwithoutobstacles（i.e.2Dchannel&noobstacles）

ifPois_test

%over-writesthedefinitionoftheporespace

clearChannel2D

Len_Channel_2D=36,%lunghezzacanale2d

Channel_2D_half_Width=8;Width=Channel_2D_half_Width*2;

Channel2D=ones（Len_Channel_2D,Width）;%definewetarea

%Channel2D（6:

12,6:

8）=0;%putfluidobstacle

imshow（Channel2D,[]）;

end

[NrMc]=size（Channel2D）;%NumberrowsandMunbercolumns

%porosity

porosity=nnz（Channel2D==1）/（Nr*Mc）

%FLUIDPROPERTIES

%physicalproperties

cs2=1/3;%

cP_visco=0.5;%[cP]1CPDinamicwaterviscosity20C

density=1.;%fluiddensity

Lky_visco=cP_visco/density;%latticekinematicviscosity

omega=（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;

%maxpoiseuillefinalvelocityontheflowprofile

uy0=-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）;;%averagefluidvelocityontheprofile

x_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）;%analyticalvelocityprofile

av_vel_t=1.e+10;%inizialization（t=0）

%PixelSize=5;%[Microns]

%dL=（Nr*PixelSize*1.0E-4）;%samplehight[cm]

%

%EXPERIMENTALSET-UP

%inletandoutletbuffers

inb=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

%boundariesrelatedtotheexperimentalsetup

wb=2;%wallthickness

Channel2D=[zeros（Nr,wb）,Channel2D,zeros（Nr,wb）];%addwalls（nofluidleak）

[NrMc]=size（Channel2D）;%updatesize

uy_analy_profile=[zeros（1,wb）,uy_analy_profile,zeros（1,wb）];%takeintoaccountwalls

x_pro_fig=[[x_profile

（1）-[wb:

-1:

1]],[x_profile,[1:

wb]+x_profile（end）]];

%Figureplotsanalyticalparabolicprofile

figure（20）,plot（x_pro_fig,uy_analy_profile,'-'）,gridon,

title（'Analyticalparab.profileforPoiseuilleplanarflowinachannel'）

%VISUALIZEPORESPACE&FLUIDOSTACLES&MEDIALAXIS

figure,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）%usedtovisualizeevolutionofrho

figure（11）%usedtovisualizeux

figure（12）%usedtovisualizeuy（i.e.top->down）

%INDICES

%Wetlocationsetc.

[iabw1jabw1]=find（Channel2D==1）;%indicesi,j,ofactivelatticelocationsi.e.pore

lena=length（iabw1）;%numberofactivelocationi.e.ofporespacelatticecells

ija=（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,ofactivelatticelocations

lenwint=length（iawint）;%numberofinternal（i.e.notborder）wetlocations

ijaint=（jawint-1）*Nr+iawint;%equivalentsingl

NxM=Nr*Mc;

%DIRECTIONS:

ENWSNENWSWSEZERO（ZERO:

RestParticle）

%y^

%625^NWNNE

%391...+x->+yWRPE

%748SWSSE

%-y

%x&ycomponentsofvelocities,+xistoest,+yistonord

East=1;North=2;West=3;South=4;NE=5;NW=6;SW=7;SE=8;RP=9;

N_c=9;%numberofdirections

%versorsD2Q9

C_x=[10-101-1-110];

C_y=[010-111-1-10];C=[C_x;C_y]

%BOUNCEBACKSCHEME

%aftercollisionthefluidelementsdensitiesfaresentbacktothe

%latticenodetheycomefromwithoppositedirection

%indicesoppositeto1:

8forfastinversionafterbounce

ic_op=[34127856];%i.e.4isoppositeto2etc.

%PERIODICBOUNDARYCONDITIONS-reinjectionrules

yi2=[Nr,1:

Nr,1];%thisdefinitionallowsimplemeningPeriodBoundCond

%yi2=[1,Nr,2:

Nr-1,1,Nr];%re-injthesecondlasttoasfirst

%directionalweights（densityweights）

w0=16/36.;w1=4/36.;w2=1/36.;

W=[w1w1w1w1w2w2w2w2w0];

%cconstants（soundspeedrelated）

cs2=1/3;cs2x2=2*cs2;cs4x2=2*cs2.^2;

f1=1/cs2;f2=1/cs2x2;f3=1/cs4x2;

f1=3.,f2=4.5;f3=1.5;%coef.ofthefequil.

%declarativestatemets

f=zeros（Nr,Mc,N_c）;%arrayoffluiddensitydistribution

feq=zeros（Nr,Mc,N_c）;%fatequilibrium

rho=ones（Nr,Mc）;%macro-scopicdensity

temp1=zeros（Nr,Mc）;

ux=zeros（Nr,Mc）;uy=zeros（Nr,Mc）;uyout=zeros（Nr,Mc）;%dimensionlessvelocities

uxsq=zeros（Nr,Mc）;uysq=zeros（Nr,Mc）;usq=zeros（Nr,Mc）;%higherdegreevelocities

%initializationarrays:

startvaluesinthewetarea

foria=1:

lena%statvaluesintheactivecellsonly;0outside

i=iabw1（ia）;j=jabw1（ia）;

f（i,j,:

）=1/9;%uniformdensitydistributionforastart

end

uy（ija）=uy0;ux（ija）=ux0;%initializefluidvelocities

rho（ija）=density;

%EXTERNAL（Body）FORCESe.g.inletpressureorinlet-outletgradient

%directions:

ENWSNENWSWSEZERO

force=-dPdL*（1/6）*1*[0-101-1-1110]';%;

%...ENESNENWSWSERP...

%thepressurepushesthefluiddowni.e.NtoS

%While..MAINTIMEEVOLUTIONLOOP

StopFlag=false;%i.e.logical（0）

Max_Iter=3000;%maxallowednumberofiteration

Check_Iter=1;Output_Every=20;%frequencyofcheck&output

Cur_Iter=0;%currentiterationcounterinizialization

toler=1.0e-8;%tollerancetodeclareconvegence

Cond_path=[];%recordingvaluesoftheconvergencecriterium

density_path=[];%recordingaver.densityvaluesforconvergence

end%endsifrestart

if（Restart==true）

StopFlag=false;Max_Iter=Max_Iter+3000;toler=1.0e-12;

end

while（~StopFlag）

Cur_Iter=Cur_Iter+1%iterationcounterupdate

%densityandmoments

rho=sum（f,3）;%density

ifCur_Iter>1%useinizializationuxuytostart

%Moments...Note:

C_x（9）=C_y（9）=0

ux=zeros（Nr,Mc）;uy=zeros（Nr,Mc）;

foric=1:

N_c-1;

ux=ux+C_x（ic）.*f（:

:

ic）;uy=uy+C_y（ic）.*f（:

:

ic）;

end

%uy=f（:

:

2）+f（:

:

5）+f（:

:

6）-f（:

:

4）-f（:

:

7）-f（:

:

8）;%inshort!

%ux=f（:

:

1）+f（:

:

5）+f（:

:

8）-f（:

:

3）-f（:

:

6）-f（:

:

7）;%inshort!

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ux（ija）=ux（ija）./rho（ija）;uy（ija）=uy（ija）./rho（ija）;

uxsq（ija）=ux（ija）.^2;uysq（ija）=uy（ija）.^2;

usq（ija）=uxsq（ija）+uysq（ija）;%

%weighteddensities:

restparticle,principalaxis,diagonals

rt0=w0.*rho;rt1=w1.*rho;rt2=w2.*rho;

%Equilibriumdistribution

%maindirections（+cross）

feq（ija）=rt1（ija）.*（1+f1*ux（ija）+f2*uxsq（ija）-f3*usq（ija））;

feq（ija+NxM*（2-1））=rt1（ija）.*（1+f1*uy（ija）+f2*uysq（ija）-f3*usq（ija））;

feq（ija+NxM*（3-1））=rt1（ija）.*（1-f1*ux（ija）+f2*uxsq（ija）-f3*usq（ija））;

%feq（ija+NxM*（3）=f（ija）-2*rt1（ija）*f1.*ux（ija）;%muchfaster...!

!

feq（ija+NxM*（4-1））=rt1（ija）.*（1-f1*uy（ija）+f2*uysq（ija）-f3*usq（ija））;

%diagonals（Xdiagonals）（ic-1）

feq（ija+NxM*（5-1））=rt2（ija）.*（1+f1*（+ux（ija）+uy（ija））+f2*（+ux（ija）+uy（ija））.^2-f3.*usq（ija））;

feq（ija+NxM*（6-1））=rt2（ija）.*（1+f1*（-ux（ija）+uy（ija））+f2*（-ux（ija）+uy（ija））.^2-f3.*usq（ija））;

feq（ija+NxM*（7-1））=rt2（ija）.*（1+f1*（-ux（ija）-uy（ija））+f2*（-ux（ija）-uy（ija））.^2-f3.*usq（ija））;

feq（ija+NxM*（8-1））=rt2（ija）.*（1+f1*（+ux（ija）-uy（ija））+f2*（+ux（ija）-uy（ija））.^2-f3.*usq（ija））;

%restparticle（.）ic=9

feq（ija+NxM*（9-1））=rt0（ija）.*（1-f3*usq（ija））;

%Collision（betweenfluidelements）omega=relaxationfrequency

f=（1.-omega）.*f+omega.*feq;

%%%%%%%%%%%%%%%%%%%%