Code covered by the BSD License  

Highlights from
immiscible LB

immiscible LB

by

 

23 Jul 2009 (Updated )

Implements Immiscible Lattice Boltzmann (ILB, D2Q9) method for two phase flows

ILB_vessels.m
%
% EXPERIMENTAL SET-UP
%Vessel = 'Channel2D'; % N to S laminar flow
%Vessel = 'Box2D'; % for bubble test
%Vessel = 'Cylinder'; % for cylindrical vessel / tank 


switch lower(Vessel)
    case {'channel2d'}
        display('Channel 2D' );
        % 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
        Cur_iter=1;
        Poiseuille_Profile;  uy_analy_profile= vet_profile;
        
        % Figure plots analytical parabolic profile
        figure(20), plot((1:Mc)-0.5,uy_analy_profile,'-'), grid on,
        title('Analytical parab. profile for Poiseuille planar flow in a channel')

        figure, imshow(Channel2D); title('Vassel geometry');
        Channel2D=logical(Channel2D);
        % obstacles for Bounce Back ( fluid 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');

        % Perimeter of dry areas
        Perim_dry=bwperim(~Channel2D,8);
        % exclude borders
        border=logical(zeros(Nr,Mc)); border(1:Nr,2:Mc-1)=1;
        border=logical(border);
        Perim_dry= Perim_dry.*(border);
        figure, imshow(Perim_dry); title('Perimeter dry');
        
        numeri=(1:1:Mc);
        jcol_free_outlet_range=(Channel2D(Nr,:)==1);
        jcol_free_outlet_range=numeri(jcol_free_outlet_range) ;
        
        
        % EXTERNAL (Body) FORCES e.g. inlet pressure or inlet-outlet gradient
        % directions: E N W S NE NW SW SE ZERO
        %Bforce = 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
        
        Define_Body_Force; 
 
    case {'box2d'}
        % vessel for bubble tests ... not implemented yet
        inb=3, oub=3; % inlet and outlet buffers thickness
        wb=3; % wall thickness

      
       
        
        if(1) % simple box 
            
         % close i.e. seal inlet and outlet
         [Nr Mc]=size(Channel2D); % update size
         inlet=zeros(inb,Mc); outlet=zeros(oub,Mc);
         Channel2D=[ [inlet]; Channel2D ;[outlet] ] ; %
         
        else % load image of porous media 
            
        % Read_Vessel_from_image_file; % reads tif file 
         inlet=zeros(inb,Mc); outlet=zeros(oub,Mc);% seal inlet & outlet 
         Channel2D=[ [inlet]; A ;[outlet] ] ;  
         
        end
        
        [Nr Mc]=size(Channel2D); % update size

        % boundaries related to the experimental set up
        Channel2D=[zeros(Nr,wb), Channel2D , zeros(Nr,wb)]; % add walls (no fluid leak)
        [Nr Mc]=size(Channel2D); % update size

        figure, imshow(Channel2D); title('Vassel geometry');
        Channel2D=logical(Channel2D);

        % obstacles for Bounce Back ( fluid in front of the grain)
        Obstacles=bwperim(Channel2D,8); % perimeter of the grains for bounce back Bound.Cond.
        figure, imshow(Obstacles); title(' Fluid obstacles (in the fluid)' );

        Medial_axis=bwmorph(Channel2D,'thin',Inf); %
        figure, imshow(Medial_axis); title('Medial axis');

        % Perimeter of dry areas for bounce back
        Perim_dry=bwperim(~Channel2D,8);
        Perim_dry=imclearborder(Perim_dry,4);
        figure, imshow(Perim_dry); title('Perimeter dry');

        % no ext Bforce for bubble test
        Bforce=0;
        
     case {'cylinder'}    
        display('Cylinder 2D' );
        % vessel for bubble tests ... not implemented yet
        inb=5; oub=5; % inlet and outlet buffers thickness
        wb=5; % wall thickness

        % close i.e. seal inlet and outlet
        inlet=zeros(inb,Mc); outlet=zeros(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
        Channel2D=[zeros(Nr,wb), Channel2D , zeros(Nr,wb)]; % add walls (no fluid leak)
        [Nr Mc]=size(Channel2D); % update size
        
        Xcentre=round(Mc/2); Ycentre=round(Nr/2); Dcyl=Nr-2*wb;
        x = [1:Mc]; x = repmat(x, Nr, 1);  
        y = [1:Nr]; y= repmat(y', 1, Mc);  
        Cylind= ( (x-Xcentre).^2 + (y-Ycentre).^2 ) <= (Dcyl/2).^2 ;
        
        Channel2D=Cylind;
        
        figure, imshow(Channel2D); title('Vassel geometry: Cylinder');
        
        Channel2D=logical(Channel2D);

        % Obstacles for Bounce Back ( fluid in front of the grain)
        Obstacles=bwperim(Channel2D,8); % perimeter of the grains for bounce back Bound.Cond.
        figure, imshow(Obstacles); title(' Fluid obstacles (in the fluid)' );

        Medial_axis=bwmorph(Channel2D,'thin',Inf); %
        figure, imshow(Medial_axis); title('Medial axis');

        % Perimeter of dry areas for bounce back
        Perim_dry=bwperim(~Channel2D,8);
        Perim_dry=imclearborder(Perim_dry,4);
        figure, imshow(Perim_dry); title('Perimeter dry');

        % no ext Bforce for bubble test
        Bforce=zeros(9,1);
        
        case {'doublet_of_pores'}    
     
%%  % doublet condotto idraulico con curve secche 
if(0)   
A=[]; A=zeros(80,60);
% condotto verticale tutta altezza
% posizione poro 
Dbv=10; Dbh=10; %  Dist v & h from border
% dimensione poro 
hvp=40; whp=40;
A(Dbv:Dbv+hvp,Dbh:Dbh+whp)=1;
Npore=10; N=Npore; ap=4, bp=6; 
% core of pore 
A(Dbv+5:Dbv+hvp-4,Dbh+3:Dbh+whp-8)=0;
% 
A(1:Dbv,20:29)=1;
A(Dbv+hvp:end,20:29)=1; 
end 
%% doublet smussato ..
if (1)
N=128; % 

%Npore=40;
xax = (1:N)-ceil(N/2); x = repmat(xax, N, 1); 
yax = (1:N)-ceil(N/2); y= repmat(yax', 1, N);  

ap=  26; bp=38;
exp=7; % esponente di prova per variare la forma della superellisse
A= (abs(x/ap).^exp+abs(y/bp).^exp) <= 1 ; % mask the inner circle of diam 

Ncwv=58; % posizione nera 
wv=14; % larghezza condotto 

A(:,Ncwv+1:N-Ncwv)=1; % condotto 

ap= 12 ; bp=29;
cx=3; cy=0.5;
exp_int=5; % esponente per variare la forma della superellisse interna
B= ((abs((x-cx)/ap).^exp_int+abs((y-cy)/bp).^exp_int) > 1 ); % mask the inner circle of diam 
A = A &  B;


ap= 16 ; bp=30; % ingrossamento del poro di diametro maggiore
cx=-17; cy=0;
exp_int=1.5; 
C= ((abs((x-cx)/ap).^exp_int+abs((y-cy)/bp).^exp_int) <= 1 ); % mask the inner circle of diam 
A = A | C;
figure, imshow(A,[]);

%A=A(:,40:end-40); % ritaglia 
%A(1:15,10:end-9)=1; % condizione "lago" sopra 
%imshow(A,[])

%A(N-40:N,20:end-20)=1; % % lago sotto 
% imshow(A,[])

A=A(1:end-10,:); % per rendere il canale pi corto
%A=A(:,70:end-70); % rifila i bordi dx e sn

figure, imshow(A,[])
end

%%

%imshow(A,[])
[Nr Mc]=size(A);
new_dims=[0.9 0.5].*[Nr, Mc] ;      
Channel2D=imresize(A,new_dims,'nearest' );        
[Nr,Mc] = size(Channel2D); % update size  
imshow(Channel2D,[])

Obstacles=bwperim(Channel2D,8); % perimeter of the grains for bounce back Bound.Cond.
Obstacles=Obstacles(3:Nr-1,:);Channel2D=Channel2D(3:Nr-1,:);[Nr Mc]=size(Channel2D);

figure, imshow(Obstacles); title(' Fluid obstacles (in the fluid)' );

Perim_dry=bwperim(~Channel2D,8);
figure, imshow(Perim_dry); title('Perimeter dry');
Medial_axis=bwmorph(Channel2D,'thin',Inf); %
figure, imshow(Medial_axis); title('Medial axis');

 case {'duct_constriction'} 
     
N=128; % 

%Npore=40;
xax = (1:N)-ceil(N/2); x = repmat(xax, N, 1); 
yax = (1:N)-ceil(N/2); y= repmat(yax', 1, N);  

ap=  16; bp=32; % assi 
exp_c=0.9;
A= (abs(x/ap).^exp_c + abs(y/bp).^exp_c) <= 1 ; % mask the inner circle of diam 

Nh=N/2; 
A=~[A(:,Nh:end),A(:,1:Nh-1) ]; 

A(:,Nh-46:Nh+46)=[]; % taglia stiscia centrale 
%imshow(A,[]) 

wA=size(A,2); A=[ones(1,wA); A; ones(1,wA)];

LA=size(A,1); 

A=[false(LA,2),A,false(LA,2)];

A=logical(A);

%imshow(A); 
[Nr Mc]=size(A);
new_dims=[1.3 1].*[Nr, Mc] ;  % nuove dimensioni       
Channel2D=imresize(A,new_dims,'nearest' );        
[Nr,Mc] = size(Channel2D); % update size  
imshow(Channel2D,[])     
     
Obstacles=bwperim(Channel2D,8); % perimeter of the grains for bounce back Bound.Cond.
Obstacles=Obstacles(3:Nr-1,:);
figure, imshow(Obstacles); title(' Fluid obstacles (in the fluid)' );

Channel2D=logical(Channel2D(3:Nr-1,:)); 
[Nr Mc]=size(Channel2D); % taglia ricalcola 


Perim_dry=bwperim(~Channel2D,8);
figure, imshow(Perim_dry); title('Perimeter dry');
Medial_axis=bwmorph(Channel2D,'thin',Inf); %
figure, imshow(Medial_axis); title('Medial axis');

Define_Body_Force; 

case {'flow_focus'} 
    
    Nr=128; Mc=65;
    A=ones(Nr,Mc); A=[zeros(Nr,2), A, zeros(Nr,2)];
    L_can=54; w_can=10; % canale Lunghezza e wide 
    sps=2;
    P_can=ceil((Mc-w_can)/2); % posizione da mrgine SX
    
    A(3:L_can,P_can:P_can+sps)=0; % SX cha wall 
    A(3:L_can,P_can+sps+w_can:P_can+w_can+2*sps)=0; %DX cha wall 
    % uscita
    %H_can=3; A(Nr-H_can:Nr,[P_can:P_can+2,P_can+w_can:P_can+w_can+2])=0;
    
    % focus costrizione
    Pos_cost=ceil(Nr/2);%posizione del nozzle in verticale da alto 
    %Pos_cost=75; 
    L_cost=0;  % allargamento  del canale di lancio 
    [Nr,Mc]=size(A);
    A((Pos_cost:Pos_cost+2), [ 1:P_can+sps-L_cost])=0; % left obstac 
    
    A((Pos_cost:Pos_cost+2), [P_can+w_can+1*sps+L_cost:Mc])=0;% right 
    
    imshow(A)
    [Nr Mc]=size(A);
    new_dims=[1.3 1].*[Nr, Mc] ;  % nuove dimensioni
    Channel2D=imresize(A,new_dims,'nearest' );
    [Nr,Mc] = size(Channel2D); % update size
    imshow(Channel2D,[])
    
    imshow(Channel2D,[])     
     
Obstacles=bwperim(Channel2D,8); % perimeter of the grains for bounce back Bound.Cond.
Obstacles=Obstacles(3:Nr-1,:);
figure, imshow(Obstacles); title(' Fluid obstacles (in the fluid)' );

Channel2D=logical(Channel2D(3:Nr-1,:)); 
[Nr Mc]=size(Channel2D); % taglia ricalcola 


Perim_dry=bwperim(~Channel2D,8);
figure, imshow(Perim_dry); title('Perimeter dry');
Medial_axis=bwmorph(Channel2D,'thin',Inf); %
figure, imshow(Medial_axis); title('Medial axis');
  
case {'flow_in_pores'} 
    
dove='C:\'
%fil_name_por='ellipse_60vf_1ab_theta20.tif'; 

Read_Vessel_from_image_file;
if(0)
fil_name_por='Nik_Bona.gif'; 
%fil_name_por='pore2d.tif'; 
%fil_name_por='ellipse_60vf_2ab_theta10.tif';
fil_name_por=[dove,fil_name_por], 
A=imread(fil_name_por); % peloso 
A=A(:,1:end); [Nr, MC]=size(A);
A=imfill(A,'holes'); 
se = strel('disk',2) ; %A=imerode(A,se);  
A=im2bw(A,0.35); A=~A; %A=bwfill(A);
A = imclose(A,se); % lisciato 
figure, imshow(A,[]);
end % Nick bona 

if(0) % nicola bona da spinelli 
fil_name_por='BW_nic.tif'; 
fil_name_por=[dove,fil_name_por];
A=imread(fil_name_por); % peloso 
A=im2bw(A);
A=A(:,1:281);  
%A=A(120:140,282:(282+70)); % piccolo 

porosity=nnz(A)./prod(size(A)); 
[Nr, Mc]=size(A);
%A=imfill(A,'holes'); 
se = strel('disk',2) ; A=imdilate(A,se);  
%A=im2bw(A,0.35); A=~A; %A=bwfill(A);
%A = imclose(A,se); % lisciato 
%A=~imfill(~A,'holes'); 
figure, imshow(A,[]);   
end

%
if(0)
se = strel('disk',8) ;       % disk, radius 15
%A = imclose(A,se); % lisciato 
A=imdilate(A,se); 
se = strel('disk',8) ; 
A=imerode(A,se); 
% close == dilation followed by an erosion,
end
%
figure, imshow(A,[]);
[Nrin,Mcin]=size(A); 

%Nr= 256; Mc=128; % overide size taken from main 
% Buffers 
Top_buf=6; Bot_buf=3;
A=A(1:Nr,1:Mc); 
A=[ones(Top_buf,Mc);A;ones(Bot_buf,Mc)]; 
[Nr,Mc]=size(A);
A=[zeros(Nr,2),A,zeros(Nr,2)];
imshow(A,[ ]); 
 
[Nr Mc]=size(A);
Magnif=[1 , 1]; 
new_dims=Magnif.*[Nr, Mc] ;  % nuove dimensioni
Channel2D=imresize(A,new_dims,'nearest' );
[Nr,Mc] = size(Channel2D); % update size
Top_buf=Top_buf.*Magnif(1); Bot_buf=Bot_buf.*Magnif(1);
imshow(Channel2D,[])

% attenzione taglio 
Obstacles=bwperim(Channel2D,8); % perimeter of the grains for bounce back Bound.Cond.
Obstacles=Obstacles(3:Nr-1,:);
figure, imshow(Obstacles); title(' Fluid obstacles (in the fluid)' );

Channel2D=logical(Channel2D(3:Nr-1,:)); 
[Nr Mc]=size(Channel2D); % taglia ricalcola 

Perim_dry=bwperim(~Channel2D,8);
figure, imshow(Perim_dry); title('Perimeter dry');
Medial_axis=bwmorph(Channel2D,'thin',Inf); %
figure, imshow(Medial_axis); title('Medial axis');

end % swithch IBL_vessels

nnz_i=sum(Channel2D,2); % Void (defined in vessel)

Contact us