Gianni SchenaJuly.docx

上传人:b****7 文档编号:11190119 上传时间:2023-02-25 格式:DOCX 页数:15 大小:19.75KB
下载 相关 举报
Gianni SchenaJuly.docx_第1页
第1页 / 共15页
Gianni SchenaJuly.docx_第2页
第2页 / 共15页
Gianni SchenaJuly.docx_第3页
第3页 / 共15页
Gianni SchenaJuly.docx_第4页
第4页 / 共15页
Gianni SchenaJuly.docx_第5页
第5页 / 共15页
点击查看更多>>
下载资源
资源描述

Gianni SchenaJuly.docx

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

Gianni SchenaJuly.docx

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;

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

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 初中教育 > 初中作文

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1