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

Red_and_Blue_Stream_D2Q9_2010.m
```
tempR=Red;
tempB=Blue;

i32_1=int32(1);

% Forward Propagation step &  Bounce Back (collision fluid with obstacles)

%sizBl=size(Blue);
sizBl=[Nr , Mc, 9];
%Blue(sizBl)=Blue(sizBl); % preallocation trick
%Red(sizBl)=Red(sizBl); % preallocation trick

for ic=1:1:N_c-1, % select velocity layer

%ic=uint32(ic);

ic2=ic_op(ic); % selects the layer of the velocity opposite to ic for BB

NMic1= NxM*(ic-1); NMmNr= NMic1-Nr;
% VECTORIAL
% section 1: cares for wet locations that are NOT on the border

% internal location i.e. away from borders
% i=iawint; j=awint; %

i2 = iawint+C_y(ic);
j2 = jawint+C_x(ic); % new i,j location
i2 = yi2(i2+i32_1); %j2 = j2(j2+1); % new i,j locations corrected for PBC

%i2j2a_ic = (j2-int32(1)).*Nr + i2 + NMic1;
%i2j2a_ic = bsxfun(@times,(j2-1),Nr) + bsxfun(@plus,i2,NMic1);

i2j2a_ic0 = j2.*Nr+i2+NMmNr;
idx=ijaint+NMic1; % ijaint defined in pre computed absolute indeces

if(1)

Blue(sizBl)=Blue(sizBl); % preallocation trick
Blue(i2j2a_ic0)=tempB(idx);

Red(sizBl)=Red(sizBl); % preallocation trick
Red(i2j2a_ic0) =tempR(idx);

else
%                             %   (VERY SLOW in 2010a)
for i=1:length(idx)
%ix=; jx=i2j2a_ic0(i);
Blue(i2j2a_ic0(i))=tempB(idx(i));
Red(i2j2a_ic0(i)) =tempR(idx(i));

end
end
% section 2: cares for wet location touching the border

i2 = iobs+C_y(ic); j2 = jobs+C_x(ic);
i2=yi2(i2+1); %j2 = j2(j2+1); % Per BC
% i2j2a where to go i.e. new location

i2j2a = (j2-1).*Nr+ i2; % new can be wet or dry
tru_Cha = (Channel2D(i2j2a)==0) ; % dry subset of new loc , logical

i2j2a_ic = i2j2a+NMic1 ; % new loc for state dep (same state)
ija_ic2  = (ijobs + NxM*(ic2-1)); % old loc but new state

% from old location going into new dry , then
% change state only

Red(ija_ic2(tru_Cha))  = tempR(ijobs(tru_Cha)+NMic1);
Blue(ija_ic2(tru_Cha)) = tempB(ijobs(tru_Cha)+NMic1);

% from old location going into new wet , then same
% state but move
notruC=~tru_Cha;

Blue(i2j2a_ic(notruC)) = tempB(ijobs(notruC)+NMic1);
Red(i2j2a_ic(notruC))  = tempR(ijobs(notruC)+NMic1);

end ; %  for ic direction```