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_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

Contact us