bibo-胳膊
% Gianni Schena July 2005, schena@
% 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(wX h_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),o nes(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
.
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
bibo-胳膊
bibo-胳膊
bibo-胳膊
bibo-胳膊
bibo-胳膊
bibo-胳膊
bibo-胳膊
bibo-胳膊
本文更新与2021-01-21 05:10,由作者提供,不代表本网站立场,转载请注明出处:https://www.bjmy2z.cn/gaokao/542585.html