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

Main_Time_Evolution_Loop_Lishc.m
         Cur_Iter=1;
        while(~StopFlag)
            Cur_Iter=Cur_Iter+1; % iteration counter update
            f=Red+Blue;
            
            Velocities_ux_uy;
            
            rhoR=sum(Red,3); rhoB=sum(Blue,3); % colour density
            
            rhoR_B(ija)=rhoR(ija)-rhoB(ija); %
            
            rhoR_B(ijPdry)=Ades_R_B; % necessary after filtering 
            
            Local_Colour_Gradient_and_Divergence; % grad calculated on rhoR_B=rhoR-rhoB; 
            
            ux_uy_Lishchuk_correct;% correction 
            
            F_BGK_EQUI_D2Q9_mat2d;
            
            %%
            %Collision (between fluid elements) omega=1 , relaxation frequency
            Fast_Red_and_Blue_Collision_step; 
            %%
           Body_Forces,
            %%

            rhoR=sum(Red,3); rhoB=sum(Blue,3); % colour density
            rho(ija)=rhoR(ija)+rhoB(ija); % (viene ric + tardi fa f)
            f=Red+Blue;
            
            %f(iawet,jawet,:)= Red(iawet,jawet,:)+Blue(iawet,jawet,:);
            %f_eq_Red_plus_Blue;

            % REDISTRUBUTION INTO Red & Blue
            %%%%%%%%%%%%%%%%
            % 2 COLOURS : MULTIPHASE SECTION
            % Two Colours ! (see Gunstensen, Keller , Rothman, Zaleski, Zanetti)
            %rhoR_B=rhoR-rhoB;
            %rhoR_B(ija)=(rhoR(ija)-rhoB(ija)); %

            % colouring the dry locations to establish adhesion or repulsion for given colour
            % R-B=1 adesion of Red , R-B=-1 adesion of Blue ,  R-B=0 Neutral
            % relevant locations are on the perimeter of the dry areas
           % rhoR_B(ijPdry)=Ades_R_B; 

           
            %Local_Colour_Gradient; % grad calculated on rhoR_B=rhoR-rhoB;
            
            %%Surface_Tension_opt; % add perturbation
            
            Select_Colour_Redistribution_Method_2010d;
            %
            % STREAM (Red & Blue separately)
%            tempR=Red; %Red_Stream_D2Q9_optim;
%            tempB=Blue; %Blue_Stream_D2Q9_optim;

           Red_and_Blue_Stream_D2Q9_2011;

           %%%
           rhoR=sum(Red,3); rhoB=sum(Blue,3); % colour density
           rho=rhoR+rhoB; %%%%
           rho(ija) = rhoR(ija)+rhoB(ija);
           f=Red+Blue;
           %%%%
          
% ends Forward Propagation &  Bounce Back Sections for Red & Blue

            % re-calculate  uy as uyout for convergence
            % monitor the average density evolution
            %vect=rho(:); vect(vect>0);cur_density=mean(vect(vect>0));
            uyout= zeros(Nr,Mc);
            for ic=1:N_c-1;
                uyout= uyout + double(C_y(ic)).*f(:,:,ic) ; % flow dim.less velocity out
            end

            % Convergence check on velocity values
            if (mod(Cur_Iter,Check_Iter)==0) ; % check for convergence every 'Check_Iter' iterations

                Cur_Iter
                % variables monitored
                % mean density and
                %vect=rho(ija); cur_density=mean(vect');
                % mean velocity
                vect=uyout(ija); av_vel_tp1= mean(vect')  ;

                Condition=abs( abs(av_vel_t/av_vel_tp1 )-1), % should --> 0

                Cond_path=[Cond_path, Condition]; % records the convergence path (value)
                %density_path=[density_path, cur_density];
                %
                av_vel_t=av_vel_tp1; %

                if (Condition < toler) || (Cur_Iter > Max_Iter)
                    StopFlag=true;
                    display( 'Stop iteration: Convergence met or iteration exeeding the max allowed' )
                    display( ['Current iteration: ',num2str(Cur_Iter),...
                        ' Max Number of iter: ',num2str(Max_Iter)] )
                    toc
                    break % Terminate execution of WHILE .. exit the time evolution loop.

                end    % if(Condition < toler | ...

            end

            if ( (mod(Cur_Iter,Output_Every)==0) || Cur_Iter<10  );  % Output from loop every ...
                % and also at each iteration between 1 and 10
                %if (Cur_Iter>60) ;  % Output from loop every ...

              %  figure(10); imshow(rho,[0.1 0.9]); title(' rho'); % visualize density evolution
               % figure(11); imshow(-ux,[ ]); title(' ux' ); % visualize fluid velocity horizontal
               % figure(12); imshow(-uy,[ ]); title(' uy' ); % visualize fluid velocity down
                %figure(14), imshow(-uyout,[]), title('uyout'); % vis vel flow out
                up=5; % linear section to visualize up from the lower row
                %figure(15), hold off, feather(ux(Nr-up,:),uy(Nr-up,:)),
                %figure(15), hold on , plot(uy_analy_profile,'r-');
               
                % show the color field composition with R/B  (french colormap)
                %rho(ija)=rhoB(ija)+rhoR(ija);
                display('recording avi frame')
                Avi_Recording2;
                
            end % every


                  end %  End main time Evolution Loop

Contact us