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

Red_and_Blue_Stream_D2Q9b.m
tic
for i=1:10000
% Forward Propagation step & % Bounce Back (collision fluid with obstacles)
%f(:,:,9) = f(:,:,9); % Rest element do not move

%temp1=f;  

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); 
    % VECTORIAL
    % vectorialize from wet location that are NOT on the border to other wet locations
  
                       
                       %temp1xB(ija) = tempB(ija+NMic1); % all wet locations 
                        %temp1xR(ija) = tempR(ija+NxM*(ic-1)); % all wet locations 
                        
                        % 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+1); %j2 = j2(j2+1); % new i,j locat corrected for PBC
                        i2j2a_ic = (j2-1).*Nr + (i2 + NMic1); %ijaint=uint32(ijaint); 
                        %i2j2a_ic = ((j2-1).*Nr + (i2 + NxM.*(ic-1)));
                        
                       % Blue(i2j2a_ic) = temp1xB(ijaint);
                       %% Red(i2j2a_ic) = temp1xR(ijaint);
                       
                        Blue(i2j2a_ic)= tempB(ijaint+NMic1); 
                        Red(i2j2a_ic) = tempR(ijaint+NMic1); 
                        
                        % VECTORIAL 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
                        %Blue(ija_ic2(tru_Cha)) = temp1xB(ijobs(tru_Cha)); %new loc = old loc
                        %Red(ija_ic2(tru_Cha)) = temp1xR(ijobs(tru_Cha)); %new loc = old loc
                       
                         Blue(ija_ic2(tru_Cha)) = tempB(ijobs(tru_Cha)+NMic1); 
                         Red(ija_ic2(tru_Cha))  = tempR(ijobs(tru_Cha)+NMic1); 
                       
                        % from old location going into new wet , then same
                        % state but move
                        notruC=~tru_Cha;
                       % Blue(i2j2a_ic(notruC)) = temp1xB(ijobs(notruC)); 
                       % Red(i2j2a_ic(notruC)) = temp1xR(ijobs(notruC)); 
                       
                        Blue(i2j2a_ic(notruC)) = tempB(ijobs(notruC)+NMic1); 
                        Red(i2j2a_ic(notruC))  = tempR(ijobs(notruC)+NMic1); 
end ; %  for ic direction
end
toc

Contact us