格子玻尔兹曼MATLAB运用（LBGK_D2Q9_poiseuille_channel2D）

格子玻尔兹曼MATLAB运用(LBGK_D2Q9_poiseuille_channel2D),二维的，适合初学者使用。

% 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=[31];wXh_Wet=[34];
 
if obs_regolare % with internal obstacles 
    
A=repmat（[zeros（wXh_Dry）ones（wXh_Wet）][13]）;A=[Azeros（wXh_Dry）];
B=ones（size（A））; 
C=[A;B]  ; D=repmat（C41）;
D=[B;D]
end
 
if obs_irregolare % with int obstacles 
A1=repmat（[zeros（wXh_Dry）ones（wXh_Wet）][13]）; 
A1=[A1zeros（wXh_Dry）]  ;
B=ones（size（A1））; 
C1=repmat（[ones（wXh_Wet）zeros（wXh_Dry）][13]）; C1=[C1ones（wXh_Dry）];
E=[A1;B;C1;B]; 
D=repmat（E21）;
D=[B;D]
end
 
if ~Pois_test
figureimshow（D[]） 
Channel2D=D;
Len_Channel_2D=size（Channel2D1）; % Length
Width=size（Channel2D2）; % 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_2DWidth）; % define wet area
%Channel2D（6:126: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 ）; % a

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件      16020  2018-04-24 10:10  D2Q92.m

----------- ---------  ---------- -----  ----

                16020                    1