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

Avi_Recording2.m
%grap_out_level=4;

CM=zeros(Nr,Mc,3);
rhoRpB=rhoR(ija)+rhoB(ija);
CM(ija)=rhoR(ija)./rhoRpB; % first layer R
%CM(ija+NxM)=0; % 2nd B
CM(ija+NxM*2)=rhoB(ija)./rhoRpB; % 3rd B

%%
waitstring=...
    ['iter.: ', num2str(Cur_Iter,'%10.0f'),' of ',num2str(Max_Iter,'%10.0f') ,' or ',...
    num2str(Cur_Iter/Max_Iter*100,'%8.2f'), ' pct.',' output each ' num2str(Output_Every,'%10.0f')];
%xlabel(waitstring);
%%

% saturation
if grap_out_level;
    Sw=zeros(Nr,1);
    Sw=sum(CM(:,:,3),2); % blue locations level 3
    %nnz_i=sum(Channel2D,2); % Void (defined in vessel)
    Sw=Sw./nnz_i; % Blue/ Void
    Swm=cumsum(Sw); Swm=Swm./(1.:1:single(Nr))'; % saturazione media
    
    % posizione minimo Sw
    [minSw,posmSw]=min(Sw(Sw>0.00001));
    
    
    
    if(1)
        
        if Bkt_flag==0;
            
            H_finger= max( posmSw, H_finger_old)
            H_finger_old=H_finger;
            if H_finger==Nr;
                Bkt_flag==1; % second check
                %(see also fractional flow at outlet)
            end
        else
            H_finger==Nr;
        end
        
        % calcolo profilo fw = QB/ Qtot
        fw=repmat(nan,Nr,1);
        %fw=zeros(Nr,1);
        for ir=max(1,H_finger-20):H_finger-1;
            iauy=find(iawet == ir);
            uyr=uy(iauy) ; %r= row
            
            %[~,ijir_r]=intersect(ija,ija(iauy)); uyr=uy(ijir_r) ; %r= row
            
            uyr(uyr<=0)=0;
            jv=find(Channel2D(ir,:)>0); % void
            rhoBr= rhoB(ir,jv)'; rhoRr=rhoR(ir,jv)';
            QwB=sum(uyr.*rhoBr); QwR=sum(uyr.*rhoRr); % scalari
            fw(ir)=QwB./(QwB+QwR);
        end
        
    end
    
    
end

if grap_out_level>1 %% Pressioni normalizzate , Relative injectivity .^1
    
    % pressiore ratio (i.e. rho(i,:) )
    Pres_rho =sum(rho,2)./max(sum(rho~=0,2),1);
    Pres_rhoB=sum(rhoB,2)./max(sum(rhoB~=0,2),1);
    Pres_rhoR=sum(rhoR,2)./max(sum(rhoR~=0,2),1);
    
    [max_Pres,posmxPre]=max(Pres_rho);  Pres_rho=Pres_rho./max_Pres;
    [max_PresB,posmxPreB]=max(Pres_rhoB); Pres_rhoB=Pres_rhoB./max_PresB;
    [max_PresR,posmxPreR]=max(Pres_rhoR); Pres_rhoR=Pres_rhoR./max_PresR;
    
    % stesso normalizzatore
    mxmxPres=max([max_Pres;max_PresB;max_PresR]);
    max_Pres=mxmxPres;max_PresB=mxmxPres;max_PresR=mxmxPres;
    % yellow [1 1 0]
    % [1 0 1] magenta
    
end,%if(1)% Pressioni normalizzate


switch grap_out_level
    case 1 % simple fluid displacement visualization
        %iptsetpref('imshowInitialMagnification','fit')
        H_avi= imshow(CM,[]);
        %set(gca,'units','normalized','outerposition',[0 0 1 1])
        colormap(flag),     %axis tight
        xlabel(waitstring); 
        axis equal 
        title([' Fluid displacement   ',  num2str(Nr) ,' x ', num2str(Mc)] ),
        
        % iframe = getframe(gca);   mov=addframe(mov,iframe); %
        
        if(0) % add border at interfaces
          hold on 
           % BW = im2bw(CM, graythresh(CM));
            BW = im2bw(CM, 0.15); % hight value for an 'internal' border
           % figure, imshow(BW,[]);
           [B,L] = bwboundaries(BW,'noholes');
          % imshow(label2rgb(L, @jet, [.5 .5 .5]))
           hold on
           for k = 1:length(B)
           boundary = B{k};
           plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 2)
          end
        end 
        
        
    case 2 % displacement & plots
        cla
        l1=0.01; b1=0.1;wd1=0.5; h1=wd1.*Nr/Mc;
        if (h1 > 0.95 | h1 < 0.5 ) , h1=0.6; wd1=h1*single(Mc)/single(Nr); end
        subplot('Position',[l1,b1,wd1,h1]);% [left bottom width height].
        imshow(CM,[]); xlabel(waitstring);
        sdx=0.05;
        l2=l1+wd1+sdx; b2=b1; wd2=0.15; h2=h1;
        
        subplot('Position',[l2,b2,wd2,h2]);% Sw [left bottom width height].
        set(gca,'nextplot','replacechildren');
        plot(Sw(end:-1:1),1:Nr,'b'); Ylim([1,Nr]); hold on
        plot(Swm(end:-1:1),1:Nr,'k'); hold on
        plot(fw(end:-1:1),(1:Nr),'c-.'); hold on
        plot([0. ; max(Swm(H_finger),0.1) ],[Nr-posmSw ; Nr-posmSw],'r')
        title('Sw, fw -.-');
        drawnow;
        l3=l2+wd2+sdx; b3=b1; wd3=wd2; h3=h1;
        
        subplot('Position',[l3,b3,wd3,h3]);% Pressure [left bottom width height].
        set(gca,'nextplot','replacechildren');
        plot(Pres_rho(end:-1:1),1:Nr,'c'); ylim([1,Nr]); xlim([0,1]);
        hold on,
        plot(Pres_rhoR(end:-1:1),1:Nr,'r'); hold on
        plot(Pres_rhoB(end:-1:1),1:Nr,'b');
        plot([0. ; 1.],[Nr-posmSw ; Nr-posmSw],'r'); % hor red sign
        set(gca, 'YTickLabel','')
        title(' P / Pmax');
        drawnow;
        
    case 3 % rotated quiver only to exploit a wide screen
        expf = 1;  % SET EXPANSION FACTOR FOR QUIVER
        %ncs  = 100; % NUMBER OF CONTOURS
        U=sqrt(ux.^2 + uy.^2);  vmax=max(U>0).^-1;  U=U.*vmax;
        %minuy=min(uy(:));maxuy=max(uy(:)); step=(maxuy-minuy)/ncs;
        %l1=0.005; b1=0.005;wd1=0.9; h1=wd1.*Nr/Mc; sdx=0.001;
        % rotated
        l1=0.005; b1=0.01; wd1=0.95; h1=wd1.*Mc/Nr;
        
        % if (h1 > 0.95 | h1 < 0.5 ) , h1=0.8; wd1=h1*Mc/Nr; end
        
        subplot('Position',[l1,b1,wd1,h1]);% [left bottom width height].
        set(gca,'nextplot','replacechildren');
        %[cs,h] = contourf(uyy,minuy:step:maxuy);
        %set(h,'edgecolor','none'); colormap(jet); colorbar
        %hold on
        %jump  densita frecce
        if(1) % (1) ROTATES
            Rotated_Red_and_Blue_Quiver;
            axis([-expf Nr+expf -expf Mc+expf])
            
            title('Outlet on the Left , Inlet on the Right side')
            
        else
            Red_and_Blue_Quiver;
            axis([-expf Mc+expf -expf Nr+expf])
            
            title(' Inlet ')
            xlabel('  Outlet'),
            
        end
        %axis tight
        %axis equal
        axis image
        
        %axis off ,
        
        
    case 4 % quiver & fluids & 
        expf = 1;  % SET EXPANSION FACTOR FOR QUIVER
        ncs  = 100; % NUMBER OF CONTOURS
        U=sqrt(ux.^2 + uy.^2);  vmax=max(U>0).^-1;  U=U.*vmax;
        %minuy=min(uy(:));maxuy=max(uy(:)); step=(maxuy-minuy)/ncs;
        l1=0.005; b1=0.1;wd1=0.5; h1=wd1.*double(Nr)./double(Mc); sdx=0.001;
        if (h1 > 0.95) , h1=0.8; wd1=h1*double(Mc)./double(Nr); end
        if (h1 < 0.4 ) , h1=0.5; wd1=h1*double(Mc)./double(Nr); end

        subplot('Position',[l1,b1,wd1,h1]);% [left bottom width height].
        set(gca,'nextplot','replacechildren');
        %[cs,h] = contourf(uyy,minuy:step:maxuy);
        %set(h,'edgecolor','none'); colormap(jet); colorbar
        %hold on
        %jump  per freccie piu rade
        Red_and_Blue_Quiver;
        axis equal
        axis([-expf Mc+expf -expf Nr+expf])
        %axis off ,
        title('Velocity arrows')
        %
        l2=l1+wd1+sdx; b2=b1; wd2=wd1; h2=h1;
        subplot('Position',[l2,b2,wd2,h2]);% [left bottom width height].
        imshow(CM,[]);  colormap(flag),  axis tight ,
       % title('Fluids'), 
        title([' Fluid displacement   ',  num2str(Nr) ,' x ', num2str(Mc)] ),
        xlabel(waitstring);
        
    case 41
        
         expf = 1;  % SET EXPANSION FACTOR FOR QUIVER
         ncs  = 100; % NUMBER OF CONTOURS
      %  U=sqrt(ux.^2 + uy.^2);  vmax=max(U>0).^-1;  U=U.*vmax;
        %minuy=min(uy(:));maxuy=max(uy(:)); step=(maxuy-minuy)/ncs;
        l1=0.005; b1=0.1;wd1=0.3; h1=wd1.*double(Nr)./double(Mc); sdx=0.001;
        if (h1 > 0.95) , h1=0.8; wd1=h1*double(Mc)./double(Nr); end
        if (h1 < 0.4 ) , h1=0.4; wd1=h1*double(Mc)./double(Nr); end

        subplot('Position',[l1,b1,wd1,h1]);% [left bottom width height].
        set(gca,'nextplot','replacechildren');
        
        U=sqrt(ux.^2 + uy.^2);  title('sqrt (ux^2 + uy^2 )');
        %U=uy;  title('uy');
        % U(ija)=uy; title(' uy ');
        U2d=nan(Nr,Mc); U2d(ija)=U(:);
       % U2d=M_Grad; 
        imshow(U2d,[]); colormap(jet), axis tight
        colorbar;
      
        
        %
        l2=l1+wd1+sdx; b2=b1; wd2=wd1; h2=h1;
        subplot('Position',[l2,b2,wd2,h2]);% [left bottom width height].
        imshow(CM,[]);  colormap(jet),  axis tight ,
       % title('Fluids'), 
        title([' Fluid displacement   ',  num2str(Nr) ,' x ', num2str(Mc)] ),
        xlabel(waitstring);
         
        
        
        
    case 5 % subplot (2,x,x) ...  sqrt (ux.^2 + uy.^2 ) etc
        
        %subplot(2,4,1) % fluid
        wimg=0.08;left=0.05; bot=0.05;half=0.5;
        himg=wimg*single(Nr)/single(Mc); dxr=0.05;
        subplot('Position',[left,half,wimg,himg]);
        imshow(CM,[]); %colormap(flag),
        axis tight,% title(' Fluid displacement '),
        title([' Fluid displacement   ',  num2str(Nr) ,' x ', num2str(Mc)] ),
        xlabel(waitstring);
        %
        %subplot(2,4,2) % quiver
        left2=left+wimg+dxr;
        
        subplot('Position',[left2,half,wimg,himg]);
        set(gca,'nextplot','replacechildren');
        expf = 1;
        U=sqrt(ux.^2 + uy.^2); U2d=zeros(Nr,Mc); U2d(ija)=U;
        vmax=max(U>0).^-1;
        % jump per frecce piu rade 
        Red_and_Blue_Quiver; axis tight
        axis image ,
        %axis equal
        axis([-expf Mc+expf -expf Nr+expf])
        %set(gca, 'XTickLabel',[],'YTickLabel',[])
        
        title('Velocity arrows')
        
        if(1)
        %subplot(2,4,3);% Sw [left bottom width height].
        left3=left2+wimg+dxr;
        wplo=wimg;
        subplot('Position',[left3,half,wplo,himg]);
        set(gca,'nextplot','replacechildren');
        
        Sw_and_fw_vertical_plot; grid on;
        
        %set(gca, 'YTickLabel','')
        drawnow;
        
        
        
        end
        
        %subplot(2,4,4);% Sw [left bottom width height].
        left4=left3+wplo+dxr;
        
        subplot('Position',[left4,half,wplo,himg]);
        set(gca,'nextplot','replacechildren');
        
        Pressures_vertical_Plot;
        drawnow;
       
        set(gca, 'YTickMode', 'Auto'),
        title('Pres ratio : P / Pmax');
        %set(gca, 'YTickLabel','')
        drawnow;
        
        % SECOND (LOWER ) LINE OF FIGURE
        
        %subplot(2,4,5) % module
        bot=0.05;
        subplot('Position',[left,bot,wimg,himg]);
        U2d=zeros(Nr,Mc); U2d(ija)=U(:);
        imshow(U2d,[]); colormap(jet), axis tight
        colorbar;
        title('sqrt (ux^2 + uy^2 )');
        
        %subplot(2,4,6)% uy
        subplot('Position',[left2,bot,wimg,himg]);
        
        U2d(ija)=uy;
        imshow(U2d,[]); colormap(jet), axis tight
        colorbar; title('uy : vertical component')
        
        subplot('Position',[left3,bot,wimg,himg]);
        %subplot(2,4,7)
        U2d(ija)=ux(:);
        imshow(U2d,[]); colormap(jet), axis tight
        colorbar; title('ux : horiz component')
        
        if(1)
        % aggiungi arrow grandi a destra 
        himg=.9; wimg=Mc/Nr;
        subplot('Position',[0.6,bot,wimg,himg]);
         
        set(gca,'nextplot','replacechildren');
        expf = 1;
        U=sqrt(ux.^2 + uy.^2); U2d=zeros(Nr,Mc); U2d(ija)=U;
        vmax=max(U>0).^-1;
        % jump per frecce piu rade 
        Red_and_Blue_Quiver; axis tight
        axis image ,
        %axis equal
        axis([-expf Mc+expf -expf Nr+expf])
        %set(gca, 'XTickLabel',[],'YTickLabel',[])
        end
        
        
    case 6 % subplot (2,x,x) ...  permeability
        
        %UPPER ... time evolution graphics
        left=0.02;bot=0.55;wimg=0.3;himg=wimg.*single(Nr)/single(Mc); dsx=0.005;
        
        subplot('Position',[left,bot,wimg,himg ]) % displacement of fluids
        imshow(CM,[]); xlabel(waitstring); %axis image 
        axis equal 
        title([' Fluid displacement   ',  num2str(Nr) ,' x ', num2str(Mc)] ),
       % title(' Fluid displacement ')
        
                
        subplot('Position',[left+wimg+dsx,bot,0.01,himg]) % plot sw and fw
        set(gca,'nextplot','replacechildren');
        plot([0,1],[in_manom, in_manom],'r-'); ylim([1,Nr]);
        hold on 
        plot([0,1],[i_fwL,i_fwL],'c-'); hold on, 
        plot([0,1],[out_manom,out_manom],'r-'); 
        plot([0.; 1],[H_finger ; H_finger],'r<');
        set(gca, 'XTickLabel',''),ylim([1,Nr]);
        set(gca,'YDir','reverse'), 
        title('pres.gauges'), ylabel(' H, finger position ' ),
        
        subplot('Position',[0.40,bot,0.05,himg]) % plot sw and fw
        set(gca,'nextplot','replacechildren');
        Sw_and_fw_vertical_plot;    %set(gca, 'YTickLabel','')
        drawnow;
        
        subplot('Position',[0.40+0.10,bot,0.05,himg]) % pressures
        set(gca,'nextplot','replacechildren');
        Pressures_vertical_Plot;
        drawnow;
        
        % % densit 
        %subplot('Position',[1-wimg-0.15,bot,wimg.*0.5,himg.*0.5 ]) % displacement of fluids
        subplot('Position',[0.58,bot-0.014,wimg.*0.5,himg.*0.7 ]) % displacement of fluids
        set(gca,'nextplot','replacechildren');
        %imshow(rho,[0 , max(rho(:))]); %xlabel(waitstring); %axis image
        
        U2d=nan(Nr,Mc); U2d(ija)=rho(ija);
        imshow(U2d,[min(rho(ija)) , max(rho(ija))]);  colormap(jet), axis tight
        %xlabel(waitstring); %axis image 
        title(' pressure ')
        axis equal , 
        cbar=colorbar;
       % cbar=colorbar('Location','EastOutside'); 
       % cpos=get(cbar,'Position'); cpos(3)=cpos(3).*0.8;  set(cbar,'Position',cpos); 
        
        
         % velocita'  
        %subplot('Position',[1-wimg-0.01,bot,wimg.*0.5,himg.*0.5 ]) % displacement of fluids
        subplot('Position',[0.74,bot-0.014,wimg.*0.5,himg.*0.7 ]) % displacement of fluids
        set(gca,'nextplot','replacechildren');
        
        %[max(uy(:)),min(uy(:))]
        
        U2d=nan(Nr,Mc);  U2d(ija)=uy;
        
        imshow(U2d,[]); colormap(jet), 
        axis tight, % colorbar;
        title('uy : vertical comp.')
        % positive if top to down 
        axis equal , colorbar;
        
        
        
        
        %LOWER ..
        left=0.05;bot=0.05;wimg=0.4;himg=0.10; dist_pl=0.05;
        
        % Flow Rates vs time 
        subplot('Position',[left,bot+2*(himg+dist_pl),wimg,himg ]) % fw versus iteration
        %set(gca,'nextplot','add');
        ylt=0.50;
        w=logspace(0,6,7);semilogx(w,0.5.*ones(1,7),'--k','LineWidth',1);
        grid on , hold on
        plot(Cur_Iter,.5.*(qv_B(1)+qv_B(2)),'b+'); grid on, hold on ,       
        plot(Cur_Iter,(qv_B(i_fwL)),'c.'); grid on, hold on ,  
        plot(Cur_Iter,(qv_R(i_fwL)),'r.'); grid on, hold on ,  
        %ylim([0 ylt]);
        ylabel('q, flow rates')
        title('volumetric flow rates at inlet (Blue) and outlet')
        
        % sotto : H finger
        subplot('Position',[left,bot,wimg,himg ]) % fw versus iteration
        
        if(0)
           
            plot(Cur_Iter,frac_flow_wat,'b.'); hold on
            plot(Cur_Iter,frac_flow_oil,'r.'); hold on
            plot(Cur_Iter,H_finger/Nr,'k.'), xlim([1 Max_Iter]);
        else
            %set(gca,'nextplot','add');
            w=logspace(0,6,7); semilogx(w,0.5.*ones(1,7),'--k','LineWidth',1);
            grid on , hold on
            semilogx(Cur_Iter,H_finger/Nr,'k.'), grid on, hold on ,
            semilogx(Cur_Iter,frac_flow_wat,'b.'); grid on; hold on ,
            semilogx(Cur_Iter,frac_flow_oil,'r.'); grid on, hold on,
            semilogx(Cur_Iter,1-cos(pi*0.5*min(Cur_Iter,tau0_uy)./tau0_uy),'g.'); grid on,
        end
        
        hold on
        
        xlabel('Iteration, i.e. time')
        ylabel('H,fwL,foL,1-cos t')
        title('fractional flows');
        
        %%
        
        % appena piu sopra : pressure drop etc
        subplot('Position',[left,bot+himg+dist_pl,wimg,himg ]) % injectivity versus iteration
        if(0)
            set(gca,'nextplot','add');
            plot(Cur_Iter, cur_DP_val,'r.'),
        else
            w=logspace(0,6,7); semilogx(w,0.5.*ones(1,7),'--k','LineWidth',1); grid on,
            hold on ,
            semilogx(Cur_Iter, cur_DP_val,'r.'); grid on, hold on,
            semilogx(Cur_Iter, P_man_in,'b.'); grid on, hold on,
            semilogx(Cur_Iter, P_man_out,'g.'); grid on, hold on,
            semilogx(Cur_Iter,Rel_Inj,'k.'); grid on, hold on,
            
        end
        
        %xlabel('Iteration, i.e. time'),
        ylabel('pressures'), 
        title('pressure-in (Blue), -out , -drop (Red) & Rel. Inject. (Black)' );
        
        if (1)
            if Bkt_flag;
                text(1.5,0.8,'Breakthrough occurred! ')
            else
                %  text(1.5,0.2,'awaiting for Breakthrough ')
            end
        end
        
        %%
        
        subplot('Position',[0.5,bot,0.4,0.4]) % permeability versus saturation after breakthrougt
        %set(gca,'nextplot','add');
        set(gca,'nextplot','replacechildren');
        if Bkt_flag;
            %set(gca,'nextplot','add');
            if (0)
                % usa saturazione media campione
                plot(Satu_sample,O_rel_per,'r.'); grid on, hold on
                plot(Satu_sample,W_rel_per,'b.'); grid on, hold on
                plot(Satu_sample,PermW_2_PermO,'kx'); grid on, hold on
                xlabel('Aver. sample Phase Saturation, Sw(Av)'),
                title('Relative Permeabilities vs Average Saturation' )
                ylabel('Kr displaced (Red) , displacing (Blue)')
                
            else
                if(1)
                    
                    plot(v_Sw,v_KrR,'r.'); grid on, hold on
                    plot(v_Sw,v_KrB,'b.'); grid on, hold on
                    %plot(Sw(i_fwL),PermW_2_PermO,'kx'); grid on, hold on
                    xlabel('Terminal Phase Saturation, Sw x=L, i.e. at outlet'),
                    title('Relative Permeabilities vs Terminal Saturation' )
                    ylabel('Kr displaced (Red) , displacing (Blue)')

                    if( Sw(i_fwL)>0.20 ), 
                        fitting_Corey; end
                    
                end
                
            end
            
            xlim([0 1]); ylim([0 1]);
            grid on,
            text(0.1,0.95,'recording after Breakthrough ')
            %text(0.1,0.9,'awaiting-for-Breakthrough ')
        else
            %set(gca,'nextplot','replacechildren');
            text(0.1,0.95,'awaiting for the Breakthrough of the invading fluid')
        end % end if Bkt_flag;
        xlabel('Terminal Phase Saturation, SwL i.e. at outlet'),
        title('Relative Permeabilities,  Krel X=K X / Kabs vs. Saturation X=Blue or X=Red' )
        
% MAIN MOVIE TITLE 
mtit('Lattice Boltzmann simulation for Relative Permeabilities (JBN method)',...
    'fontsize',14,'color',[1 0 0],...
	     'xoff',-.2,'yoff',.1); hold on,
%
    case 9 % grap_out_level =9 
        
        grap_out_level_9; 
        % perimeter of bubble for changing surface tension 
        

    otherwise
        disp('Unknown or invalid graphical output parameter.')
end


%

iframe = getframe(gcf);
mov=addframe(mov,iframe); % add to recording movie

pause (0.1)

Contact us