基于MATLAB实现的LBE模型在二维多孔介质流体渗流中的应用，基于geometry+使用说明文档.zip

% Gianni Schena July 2005, schena@units.it % Lattice Boltzmann LBE, geometry: D2Q9, model: BGK % Application to permeability in porous media Restart=false % to restart from an earlier convergence logical(Restart); if Restart==false; close all, clear all % start from scratch and clean ... Restart=false; % type of channel geometry ; % one of the flollowing flags == true Pois_test=true, % no obstacles in the 2D channel % porous systems obs_regolare=false % obs_irregolare=false % tic % IN % |vvvv| + y % |vvvv| ^ % |vvvv| | -> + x % OUT % Pores in 2D : Wet and Dry locations (Wet ==1 , Dry ==0 ) wXh_Dry=[3,1];wXh_Wet=[3,4]; if obs_regolare, % with internal obstacles 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 if obs_irregolare, % with int obstacles 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); % should not be hod Channel_2D_half_Width=Width/2, end % test without obstacles (i.e. 2D channel & no obstacles) if Pois_test %over-writes the definition of the pore space clear Channel2D Len_Channel_2D=36, % lunghezza canale 2d Channel_2D_half_Width=8; Width=Channel_2D_half_Width*2; Channel2D=ones(Len_Channel_2D,Width); % define wet area %Channel2D(6:12,6:8)=0; % put fluid obstacle imshow(Channel2D,[]); end [Nr Mc]=size(Channel2D); % Number rows and Munber columns % porosity porosity=nnz(Channel2D==1)/(Nr*Mc) % FLUID PROPERTIES % physical properties cs2=1/3; % cP_visco=0.5; % [cP] 1 CP Dinamic water viscosity 20 C density=1.; % fluid density Lky_visco=cP_visco/density; % lattice kinematic viscosity omega=(Lky_visco/cs2+0.5).^-1; % omega: relaxation frequency %Lky_visco=cs2*(1/omega - 0.5) , % lattice kinematic viscosity %dPdL= Pressure / dL;% External pressure gradient [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); % Poiseuille Gradient; % max poiseuille final velocity on the flow profile uy0=-0.001; ux0=0.0001; % linear vel .. inizialization % % uy_fin_max=-0.2; % max poiseuille final velocity on the flow profile % omega=0.5, cs2=1/3; % omega: relaxation frequency % Lky_visco=cs2*(1/omega - 0.5) , % lattice kinematic viscosity % dPdL = abs( 2*Lky_visco*uy_fin_max/(Channel_2D_half_Width.^2) ); % Poiseuille Gradient; % uyf_av=uy_fin_max*(2/3);; % average fluid velocity on the profile 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 ); % analytical velocity profile av_vel_t=1.e+10; % inizialization (t=0) %PixelSize= 5; % [Microns] %dL=(Nr*PixelSize*1.0E-4); % sample hight [cm] % % EXPERIMENTAL SET-UP % inlet and outlet buffers inb=2, oub=2; % inlet and outlet buffers thickness % add fluid at the inlet (top) and outlet (down) inlet=ones(inb,Mc); outlet=ones(oub,Mc); Channel2D=[ [inlet]; Channel2D ;[outlet] ] ; % add flux in and down (E to W) [Nr Mc]=size(Channel2D); % update size % boundaries related to the experimental set up wb=2; % wall thickness Channel2D=[zeros(Nr,wb), Channel2D , zeros(Nr,wb)]; % add walls (no fluid leak) [Nr Mc]=size(Channel2D); % update size uy_analy_profile=[zeros(1,wb), uy_analy_profile, zeros(1,wb) ] ; % take into account walls x_pro_fig=[[x_profile(1)-[wb:-1:1]], [x_profile, [1:wb]+x_profile(end)] ]; % Figure plots analytical parabolic profile figure(20), plot(x_pro_fig,uy_analy_profile,'-'), grid on, title('Analytical parab. profile for Poiseuille planar flow in a channel') % VISUALIZE PORE SPACE & FLUID OSTACLES & MEDIAL AXIS figure, imshow(Channel2D); title('Vassel geometry'); Channel2D=logical(Channel2D); % obstacles for Bounce Back ( in front of the grain) Obstacles=bwperim(Channel2D,8); % perimeter of the grains for bounce back Bound.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(' Fluid obstacles (in the fluid)' ); % Medial_axis=bwmorph(Channel2D,'thin',Inf); % figure, imshow(Medial_axis); title('Medial axis'); figure(10) % used to visualize evolution of rho figure(11) % used to visualize ux figure(12) % used to visualize uy (i.e. top -> down) % INDICES % Wet locations etc. [iabw1 jabw1]=find(Channel2D==1); % indices i,j, of active lattice locations i.e. pore lena=length(iabw1); % number of active location i.e. of pore space lattice cells ija= (jabw1-1)*Nr+iabw1; % equivalent single index (i,j)->> ija for active locations % absolute (single index) position of the obstacles in for bounce back in Channel2D % Obstacles [iobs jobs]=find(Obstacles);lenobs=length(iobs); ijobs= (jobs-1)*Nr+iobs; % as above % Medial axis of the pore space [ima jma]=find(Medial_axis); lenma=length(ima); ijma= (jma-1)*Nr+ima; % as above % Internal wet locations : wet & ~obstables % (i.e. internal wet lattice location non in contact with dray locations) [iawint jawint]=find(( Channel2D==1 & ~Obstacles)); % indices i,j, of active lattice locations lenwint=length(iawint); % number of internal (i.e. not border) wet locations ijaint= (jawint-1)*Nr+iawint; % equivalent singl NxM=Nr*Mc; % DIRECTIONS: E N W S NE NW SW SE ZERO (ZERO:Rest Particle) % y^ % 6 2 5 ^ NW N NE % 3 9 1 ... +x-> +y W RP E % 7 4 8 SW S SE % -y % x & y components of velocities , +x is to est , +y is to nord East=1; North=2; West=3; South=4; NE=5; NW=6; SW=7; SE=8; RP=9; N_c=9 ; % number of directions % versors D2Q9 C_x=[1 0 -1 0 1 -1 -1 1 0]; C_y=[0 1 0 -1 1 1 -1 -1 0]; C=[C_x;C_y] % BOUNCE BACK SCHEME % after collision the fluid elements densities f are sent back to the % lattice node they come from with opposite direction % indices opposite to 1:8 for fast inversion after bounce ic_op = [3 4 1 2 7 8 5 6]; % i.e. 4 is opposite to 2 etc. % PERIODIC BOUNDARY CONDITIONS - reinjection rules yi2=[Nr , 1:Nr , 1]; % this definition allows implemening Period Bound Cond %yi2=[1, Nr , 2:Nr-1 , 1,Nr]; % re-inj the second last to as first % directional weights (density weights) w0=16/36. ; w1=4/36. ; w2=1/36.; W=[ w1 w1 w1 w1 w2 w2 w2 w2 w0]; %c constants (sound speed related) 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. of the f equil. % declarative statemets f=zeros(Nr,Mc,N_c); % array of fluid density distribution feq=zeros(Nr,Mc,N_c); % f at equilibrium rho=ones(Nr,Mc); % macro-scopic density temp1=zeros(Nr,Mc); ux=zeros(Nr,Mc); uy=zeros(Nr,Mc); uyout=zeros(Nr,Mc); % dimensionless velocities uxsq=zeros(Nr,Mc); uysq=zeros(Nr,Mc); usq=zeros(Nr,Mc); % higher degree velocities % initialization arrays : start values in the wet area for ia=1:lena % stat values in the active cells only ; 0 outside i=iabw1(ia); j=jabw1(ia); f(i,j,:)=1/9; % uniform density distribution for a start end uy(ija)=uy0; ux(ija)=ux0; % initialize fluid velocities rho(ija)=density; % EXTERNAL (Body) FORCES e.g. inlet pressure or inlet-outlet gradient % directions: E N W S NE NW SW SE ZERO force = -dPdL*(1/6)*1*[0 -1 0 1 -1 -1 1 1 0]'; %; %... E N E S NE NW SW SE RP ... % the pressure pushes the fluid down i.e. N to S % While .. MAIN TIME EVOLUTION LOOP StopFlag=false; % i.e. logical(0) Max_Iter=3000; % max allowed number of iteration Check_Iter=1; Output_Every=20; % frequency of check & output Cur_Iter=0; % current iteration counter inizializa