Code covered by the BSD License

# immiscible LB

by

### G. Sken (view profile)

23 Jul 2009 (Updated )

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

D2Q9_Lattice.m
```% DIRECTIONS: E N W S NE NW SW SE RP (RP:Rest Particle)
%    y^
%  6 2 5           ^         NW  N  NE
%  3 9 1 ... +x-> +y         W   RP  E
%  7 4 8                     SW  S  SE
%   -y
% x & y components of velocities , +x is to est , +y is to south
% Notice. Matlab origin 0,0 is top left
East=1; North=2; West=3; South=4; NE=5; NW=6; SW=7; SE=8; RP=9;
N_c=9 ; % number of directions
% versors D2Q9
C_x=int32([1, 0, -1, 0, 1, -1, -1,  1, 0]);
C_y=int32([0, 1, 0, -1, 1,  1, -1, -1, 0]);
C2=[C_x;C_y];

if(0)
Cx_2d=zeros(3,3); Cx_2d=[-ones(3,1),zeros(3,1),ones(3,1)]; % for Colour Gradient x calculation
Cy_2d=zeros(3,3); Cy_2d=[ones(1,3);zeros(1,3);-ones(1,3)]; % For CGy calculation
else
Cx_2d=zeros(3,3); Cx_2d(2,1)=-1; Cx_2d(2,3)=1;
Cy_2d=zeros(3,3); Cy_2d(1,2)=1; Cy_2d(3,2)=-1;
end

% BOUNCE BACK SCHEME
% after collision the fluid elements densities f are sent back to the
% lattice node they come from with opposite direction
% indices opposite to 1:8 for fast inversion after bounce

ic_op = int32([3 4 1 2 7 8 5 6 9]); %   i.e. 4 is opposite to 2 etc.

ic_90= [2 3 4 1 6 7 8 5 ]; %ic_90(8)=5 i.e.5 is 90 degrees anti-clockwise from 8
ic_rad=single([0, pi/2, pi, 3/2*pi, pi/4 , 3/4*pi 5/4*pi 7/4*pi]); % 8 ic angles (rad)

% cos(2*(ic_rad-CG_phi)) is mass and moment conserving

% PERIODIC BOUNDARY CONDITIONS PBC - reinjection rules
yi2=int32([Nr , 1:Nr , 1]'); % this definition allows implemening Period Bound Cond
xi2=int32([Mc , 1:Mc , 1]');
%yi2=[1, Nr , 2:Nr-1 , 1,Nr]; % re-inj the second last to as first

% directional weights (density weights)
w0=16/36. ; w1=4/36. ; w2=1/36.;
W=[ w1 w1 w1 w1 w2 w2 w2 w2 w0]; gi=[-2. -2. -2. -2.  4. 4. 4. 4. 1];
wigi=W.*gi;
% tuning for micro-curents due 2 surface tension pert
st_lambda=[sqrt(3)/2*[1 1 1 1] 1 1 1 1 ]; % see Dupin, Holliday, Care
%c constants (sound speed related)
cs2=1/3; cs2x2=2*cs2; cs4x2=2*cs2.^2;
f1=1/cs2; f2=1/cs2x2; f3=1/cs4x2;
f1=3.; f2=4.5; f3=1.5; % coef. of the f equil.