胳膊格子玻尔兹曼MATLAB运用(LBGK_D2Q9_poiseuille_channel2D)

2021年1月21日发(作者：abracadabra)
% 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