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_Demo_ILBGK_costriction.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 
dir_avi='D:\AVI_ILB\'; 
avi_fname='ImLaBo_t'; % Immiscible Lattice Boltzmann Test 
avi_fname=[dir_avi, avi_fname]; 

grap_out_level=1;

Len_Channel_2D=128; % 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=-1.e-2 ; % 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.39; % volumetric fraction of Red fluid
Sigma=0.000002; % surface tension parameter
Ades_R_B=-1; 



Red_and_Blue_Fluids; % properties

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

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';

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=55000; % max allowed number of iteration
        Check_Iter=200; Output_Every=400; % 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
        
        
        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
       

Contact us