LBM模拟多孔介质流

clear
%---------------------
%随机四参数生成多孔介质
%----------------------

% GENERAL FLOW CONSTANTS
set（0‘defaultfigurecolor‘‘w‘）;
lx     = 400;      % number of cells in x-direction
ly     = 100;      % number of cells in y-direction
obst_x = lx/5+1;   % position of the cylinder; （exact
obst_y = ly/2+3;   % y-symmetry is avoided）
obst_r = ly/10+1;  % radius of the cylinder
uMax   = 0.1;      % maximum velocity of Poiseuille inflow
Re     = 100;      % Reynolds number
nu     = uMax * 2.*obst_r / Re;  % kinematic viscosity
omega  = 1. / （3*nu+1./2.）;      % relaxation parameter
maxT   = 400;  % total number of iterations
tPlot  = 10;      % cycles
d14 = 0.1; 
d58 = d14/4;  
n = 0.99;         %（1-n）为孔隙度                                                          
cdd = 0.01;       %生长概率                                                         
numtotal_need= （1-n） * lx * ly;     %所需生长相总数                                          
numsoild = 0;  
arrgrid = zeros;                                                            
soild = zeros;                                                               
 
 
for i=1: lx  
   for j = 1:ly 
     if rand（ ） < cdd    %随机函数来生成生长相
       numsoild= numsoild+ 1; 
       arrgrid（ij） = 1; 
       soild（numsoild 1） = i; %储存生长相的i坐标值
       soild（numsoild 2） = j; %储存生长相的j坐标值
        
     end % end if 
   end % end j 
end % end i 
Tnumsoild= numsoild; 
 
 
while Tnumsoild< numtotal_need                                                                 
   for index_soild= 1: Tnumsoild 
    index_i= soild（index_soild 1）; 
    index_j= soild（index_soild 2）; 
  
     if index_j< ly 
       i= index_i;  %生长相向右生长
       j = index_j + 1;  
       if （arrgrid（ i j ）== 0 && rand（ ） < d14） 
          numsoild= numsoild+ 1; arrgrid（ i j ） = 1; 
          soild（ numsoild 1） = i; soild（numsoild 2） = j; 
       end  
     end 
 
     if index_i> 1 
       i= index_i-1; j = index_j; % 生长相向上生长
       if （arrgrid（ i j ）== 0 && rand（ ） < d14） 
          numsoild= numsoild+ 1; arrgrid（ i j ） = 1; 
          soild（ numsoild 1） = i; soild（numsoild 2） = j; 
       end  
     end  
 
     if index_j> 1 
       i= index_i; j = index_j-1; % 生长相向左生长
       if （arrgrid（ i j ）== 0 && rand（ ） < d14） 
          numsoild= numsoild+ 1; arrgrid（ i j ） = 1; 
          soild（ numsoild 1） = i; soild（numsoild 2） = j; 
       end  
     end  
 
     if index_i< lx 
       i= index_i+1; j = index_j; % 生长相向下生长
        if （arrgrid（ i j ）== 0 && rand（ ） < d14） 
          numsoild= numsoild+ 1; arrgrid（ i j ） = 1; 
          soild（ numsoild 1） = i; soild（numsoild 2） = j; 
       end  
     end  
 
 
     if index_j< ly && index_i>1 
       i= index_i-1; j = index_j+1; % 生长相向左下生长
        if （arrgrid（ i j ）== 0 && rand（ ） < d14） 
          numsoild= numsoild+ 1; arrgrid（ i j ） = 1; 
          soild（ numsoild 1） = i; soild（numsoild 2） = j; 
       end  
     end