Gianni SchenaJulyWord文件下载.docx
《Gianni SchenaJulyWord文件下载.docx》由会员分享，可在线阅读，更多相关《Gianni SchenaJulyWord文件下载.docx（15页珍藏版）》请在冰点文库上搜索。
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）;
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,
%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,[]）;
[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;
%Lky_visco=cs2*（1/omega-0.5）,%latticekinematicviscosity
%dPdL=abs（2*Lky_visco*uy_fin_max/（Channel_2D_half_Width.^2））;
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）
%updatesize
%boundariesrelatedtotheexperimentalsetup
wb=2;
%wallthickness
Channel2D=[zeros（Nr,wb）,Channel2D,zeros（Nr,wb）];
%addwalls（nofluidleak）
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）;
Fluidobstacles（inthefluid）'
）;
Medial_axis=bwmorph（Channel2D,'
thin'
Inf）;
%
figure,imshow（Medial_axis）;
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;
%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）
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
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;
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）;
foric=1:
N_c-1;
ux=ux+C_x（ic）.*f（:
:
ic）;
uy=uy+C_y（ic）.*f（:
end
%uy=f（:
2）+f（:
5）+f（:
6）-f（:
4）-f（:
7）-f（:
8）;
%inshort!
%ux=f（:
1）+f（:
8）-f（:
3）-f（:
7）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;
%%%%%%%%%%%%%%%%%%%%
升级会员 