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

Run_Maze_solution_ILB.m
% Gianni Schena  March 2007, schena@units.it
% Lattice Boltzmann LBE, geometry: D2Q9, model: BGK
% 2D channel with Poiseuille flow through ( North to South )
% 2 immiscible fluids R-K model (Rothman - Keller)

close all, clear all % start from scratch and clean ...
tic

%   IN
% |vvvv|    - y
% |vvvv|     |
% |vvvv|     | -> + x
%  OUT       v +y
% FLOW
% Laminar parabolic flow in a 2d channel

% AVI File Mame to be generated 
avi_fname='maze_avi'; % Immiscible Lattice Boltzmann Test 
dir_avi='C:\'; 
avi_fname=[dir_avi, avi_fname]; 
grap_out_level=41; % level of graphics output
%

Len_Channel_2D=64; % Length
Channel_2D_half_Width=16; % Half_Width
Width=Channel_2D_half_Width*2+1; % add 1 for symmetry 
Channel2D=ones(Len_Channel_2D,Width); % define wet area
[Nr Mc]=size(Channel2D); % Number rows and Munber columns
uy0=0; ux0=0; %  linear vel .. inizialization

uy_fin_max=-6.e-1 ; % max poiseuille final  velocity on the flow profile

uyf_av=uy_fin_max*(2/3); % average fluid velocity on the profile
omega=1; cs2=1/3; % omega: relaxation frequency
Lky_visco=cs2*(1/omega - 0.5) ; % lattice kinematic viscosity

dPdL = ( 2*Lky_visco*uy_fin_max/(Channel_2D_half_Width.^2) ); % Poiseuille Gradient;

%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

Poiseuille_Profile; % uy_analy_profile= ... 

PixelSize= 5; % [Microns]
dL=(Nr*PixelSize*1.0E-4); % sample hight [cm]

% COLOUR REFILL METHOD ALTERNAVIVES 

%Refill_method='Latvo_Kakko'  ; % see M.Latvo kakko and Rothmann papers 
Refill_method='Keller'  ; % see Keller and Rothmann paper 


fct_Red=0.01; % volumetric fraction of Red fluid
Sigma=0.0001; % surface tension parameter
Ades_R_B=-1; 



Red_and_Blue_Fluids; % properties

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

ILB_vessels; % include file

%figure(10) % used to visualize evolution of rho
%figure(11) % used to visualize ux
%figure(12) % used to visualize uy (i.e. top -> down)

Compute_Absolute_Indices_for_Wet_Locations; 

D2Q9_Lattice; % Lattice & related matter 

Red_and_Blue_Fluids; 

Initialize_Velocities; % vectors linear vel .. inizialization


% INITIAL COLOUR DISTRIBUTION INTO WET LOCATIONS
% initialization arrays : start Kolour values in the wet area
%Initial_Colour_Distribution = 'Rand_Unif'
%Initial_Colour_Distribution ='Vert_Interf';
%Initial_Colour_Distribution ='Horiz_Interf';
%Initial_Colour_Distribution ='Red_Lane';
%Initial_Colour_Distribution ='Red_Buble2D';
%Initial_Colour_Distribution = 'red_stripe';
Initial_Colour_Distribution = 'red_spot';

ILB_initial_RB_colours; % include file 

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % distributions
        % Blue or Red !
        f=Red+Blue; % f is 0 in the dry areas
        rhoR=sum(Red,3); rhoB=sum(Blue,3); % colour density
        rho=rhoR+rhoB; %

        Avi_Initialization2; 
        
        uy(ija)=uy0; ux(ija)=ux0; % initialize fluid velocities

        % While .. MAIN TIME EVOLUTION LOOP
        av_vel_t=1.e+10; % inizialization (t=0)
        StopFlag=false; % i.e. logical(0)
        Max_Iter=1000; % max allowed number of iteration
        Check_Iter=100; Output_Every=100; % frequency of check & output
        Cur_Iter=0; % current iteration counter inizialization
        toler=1.0e-7; % tollerance to declare convegence
        Cond_path=[]; % recording values of the convergence criterium
        density_path=[]; % recording aver. density values for convergence
        
        %
        Define_Body_Force;
        
        Main_Time_Evolution_Loop,
        
        % Output & Draw after the end of the time evolution

        figure, plot(Cond_path(2:end)); title('convergence path')
        %figure, plot(density_path(2:end)); title('density convergence path')
        %figure, plot( [uy(Nr-up,:)-uy_analy_profile] ); title('LB - Analytical solution')

        toc
     Close_Avi
     
%     Red_Bubbles_size_distribution
     
       

Contact us