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

Local_Colour_Gradient_and_Divergence.m
 % Local Colour Gradient CG 
 
           min_M_Grad =1.e-3;
 
            % max gradient is perpend to R/B interface
            % local i.e. microscopic colour field 
            % Gradient by linear filtering rhoRed - rhoBlue 
            %
            %rhoR_B(ijPdry)=Ades_RNB; % enforce wettability / adesivity
            % just done 
            
            if(0)
           %rhoR_B(ijPdry)=Ades_R_B;
           rhoR_B_1=[rhoR_B(end:-1:end-1,:);rhoR_B;rhoR_B(1:2,:)];
           CGx=filter2(Cx_2d, rhoR_B_1) ; % Grad x
           CGy=filter2(-Cy_2d,rhoR_B_1) ; % Grad y 
           
           CGx=CGx(3:end-2,:);
           CGy=CGy(3:end-2,:);
            end 
            
            %
            if(0) % 
            CGx=imfilter(rhoR_B, Cx_2d); % Grad x
            CGy=imfilter(rhoR_B,-Cy_2d); % Grad y 
            rhoR_B(ijPdry)=Ades_R_B;
            end
            
            
            if(1) % 
             rhoR_B=single(rhoR_B); 
            % figure, imshow(rhoR_B,[]);
         %   CGx=imfilter([rhoR_B(end:-1:end-1,:);rhoR_B;rhoR_B(1:2,:)], Cx_2d); % Grad x
         %   CGy=imfilter([rhoR_B(end:-1:end-1,:);rhoR_B;rhoR_B(1:2,:)],-Cy_2d); % Grad y 
             %CGx=CGx(3:end-2,:);
            %CGy=CGy(3:end-2,:);
            
            CGx=imfilter(rhoR_B, Cx_2d,'circular','same'); % Grad x
            CGy=imfilter(rhoR_B,-Cy_2d,'circular','same'); % Grad y 
            
           % rotate 180 for conv2 
            %CGx=conv2(rhoR_B,-Cx_2d,'same'); % Grad x
           % CGy=conv2(rhoR_B,+Cy_2d,'same'); % Grad y 
         
           
            end
            
            
            
            if(0)
            %use ROIFILT2(H,I,BW)
            Channel2D_1=[Channel2D(end:-1:end-1,:);Channel2D;Channel2D(1:2,:)];
            CGx=roifilt2(Cx_2d, [rhoR_B(end:-1:end-1,:);rhoR_B;rhoR_B(1:2,:)], Channel2D_1); % Grad x
            CGy=roifilt2(-Cy_2d,[rhoR_B(end:-1:end-1,:);rhoR_B;rhoR_B(1:2,:)], Channel2D_1); % Grad y 
            
            CGx=CGx(3:end-2,:);
            CGy=CGy(3:end-2,:);
            
            
            %CGx = roifilt2(Cx_2d,rhoR_B,Channel2D);
            %CGy = roifilt2(-Cy_2d,rhoR_B,Channel2D);
            end
            
            % Note the change of sign: first row of rhoR_B is on top of 2nd
            % row and alike for 3rd threfore the -1 row of Cy_2d is over the +1
           % CGx=single(CGx) ; CGy=single(CGy);
            
            M_Grad=sqrt(CGx.^2+CGy.^2); % Mod CG (Magnitude )
            %M_Grad=abs(CGx) + abs(CGy); % Mod CG (faster Magnitude ext. as abs value)
           
         %   figure, imshow(M_Grad,[]); title('M grad');
            
          %  figure, imshow(single(CGy),[]);
            
            % Normalization 
           %  CGy=CGy./(M_Grad+single(M_Grad==0));   CGx=CGx./(M_Grad+single(M_Grad==0)); 
          
           % CGy=CGy./M_Grad;   CGx=CGx./M_Grad; 
            
            phi_Grad= atan2(CGy,CGx); % angle with +x (East) direction
            
            %Conc_fct(ija)=1.-abs( rhoR_B(ija)./rho(ija) ); % conc gradient 
            
         % CGy=CGy./M_Grad;   CGx=CGx./M_Grad; 
            
        % CGy=CGy./(M_Grad+single(M_Grad==0));    CGx=CGx./(M_Grad+single(M_Grad==0));  
       
 
    % divergence 
    CGxx=imfilter(CGx, Cx_2d,'circular','same'); % Grad
    CGyy=imfilter(CGy,-Cy_2d,'circular','same'); % Grad y ^2
        

    CGxy=imfilter(CGx, -Cy_2d,'circular','same'); % Grad y Grad x
    CGyx=imfilter(CGy,  Cx_2d,'circular','same'); % Grad x Grad y 

    K_Lish = CGx.*CGy.*(CGyx + CGxy) -(CGy.^2).*CGxx - (CGx.^2).*CGyy ;
   
 %   K_Lish= -(CGxx + CGyy) ;  %figure, imshow(K_Lish,[]); title('lishchuk');
    
%figure, imshow(K_Lish,[]); title('lishchuk');

% Angles of the Grad measured in anti-clockwise % m


ijagm0=find(phi_Grad(:)<0); 
phi_Grad(ijagm0) = 2*pi+phi_Grad(ijagm0);

% where in the active wet location the gradient is ... 

%Flag_Grad=false(lena,1);

%Flag_Grad =( M_Grad(ija)-min_M_Grad > 0 & rhoR(ija) >0 & rhoB(ija) >0); 
% OR 
Flag_Grad = M_Grad(ija)-min_M_Grad > 0; 

ijagr=ija(Flag_Grad); 
lenag=int32(length(ijagr));

tC2=double(C2(:,1:8));
Grad_List=[CGx(ijagr) , CGy(ijagr)]*tC2;
CGx(~ijagr)=0; CGy(~ijagr)=0; 

K_Lish(~ijagr)=0;
%figure, imshow(K_Lish,[]); title('lishchuk');

Contact us