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;
%%%%%%%%%%%%%%%%%%%%
升级会员 