%% main program for 1-D foam transport in porous media %%

clear all
clc

    NX=50; %total number of grid blocks in x direction       
    deltax=1/NX; %step length in x direction
    deltat=0.005*deltax; %step length in t direction
    NT=round(5/deltat); %total number of grid blocks in t direction
    TIME_TOTAL=5; %total dimensionless time
    
    T_max=0; %surfactant adsorption plateau (mg/m^2)
    Tr_max=0; %tracer (NaCl) adsorption plateau (mg/m^2)
    alpha=12.2; %parameter in Langmuir isotherm (L/g)
    poro=0.36; %porosity of the formation
    rou_s=2.71; %density of the formation (g/cm^3)
    a_s=0.1; %specific surface area of the formation (m^2/g)
    
    Cg_inj=0;%injected surfactant concentration in CO2 phase (g/L)
    Psgw=0; %partition coefficient between gas phase and aqueous phase 
    Cw_inj=10; %injected surfactant concentration in aqueous phase (g/L)
    NaCl_w_0=1.0; %injected  tracer cocentration (dimensionless)
  
    f_g_BC=1; %injected gas fraction (foam quality)
    
    S_wc=0.07; %connate water saturation  
    S_gr=0; %residual gas saturation
    krw0=0.79; %end-point relative permeability for aqueous phase
    krg0=1.0; %end-point relative permeability for gas phase
    visw=0.001*1e3; %in cp
    visg=0.00002*1e3; %in cp
    nw=1.96; %exponent in water relative permeability function
    ng=2.29; %exponent in gas relative permeability function
    fmmob=47094; %reference mobility reduction factor
    fmdry=0.1006; %critical water saturation for foam collapse
    epdry=500; %dry-out parameter
    
    S_g=zeros(1,NX); % initialize gas phase saturation (invading phase) 
    S_g(1,:)=S_gr; %IC for gas saturation
    S_w=1-S_g; %IC for water saturation 

    Csw=Cw_inj*ones(1,NX); %IC for aqueous urfactant concentration(g/L)
    Css=rou_s*a_s*T_max*alpha*Csw./(1+alpha*Csw); % IC for surfactant adsorption(g/L)

    NaCl_w=zeros(1,NX); %IC for aqueous tracer concentration (dimensionless)   
    NaCl_s=rou_s*a_s*Tr_max*alpha*NaCl_w./(1+alpha*Cw_inj*NaCl_w);       
      
    Fwater=0.5+(atan(epdry*(S_w-fmdry)))/pi; %dry-out function
    fmsurf=10; %minimal surfactant concentration (g/L) to reach full strength foam
    epsurf=1; %exponent in surfactant-cocentration-dependent function
    Fsurf=(Csw/fmsurf).^epsurf; %contribution of surfactant to foam stability    
    FM=1./(1+fmmob*Fwater.*Fsurf); %mobility reduction factor in the foam model

    S_nor=(S_g-S_gr)/(1-S_gr-S_wc);%normalized gas phase saturation (invading phase)
    kri=FM.*(S_nor.^ng); %IC for invading (gas) phase relative permeability
    krd=(1-S_nor).^nw; %IC for displaced (aqueous) phase relative permeability
    M=krg0/visg/(krw0/visw); %end-point mobility ratio 
    lamda=kri+krd/M; %IC for lamda(t,x) = kri + krd/M
       
    pot=zeros(1,NX);  %initialize flow potential
    pc=zeros(1,NX);  %initialize dimensionless capillary pressure       
    Npc=0; %capillary number in capillary pressure model
    C1=0.17; %parameter C1 in capillary pressure model
    C2=0.29; %parameter C2 in capillary pressure model

    u_w=(1-f_g_BC)*ones(1,NX+1); %dimensionless aqueous phase velocity at the WALL of each grid block
    %u_w(I-1/2) is in u_w(I) and u_w(I+1/2) is in u_w(I+1)
    u_g=f_g_BC*ones(1,NX+1); %dimensionless gas phase velocity at the WALL of each grid block

    COUNTER=1;
    TIME=0;
    TIME_HIS(COUNTER)=TIME;
    f_g_HIS(COUNTER)=0;
    f_surf_HIS(COUNTER)=0;
    f_NaCl_HIS(COUNTER)=0;     
    f_total_HIS(COUNTER)=0;
    p_drop_HIS(COUNTER)=0;
    VIS_FOAM(COUNTER)=0; %foam apparent viscosity
 
    Pe_w=1e100; %Define Peclect number for aqueous dispersion Pe= L*uBC/phi/Ksw
    Pe_g=1e100; %Define Peclect number for gaseous dispersion Pe= L*uBC/phi/Ksg       

    SAG=0; %SAG=1 if Surfactant-Alternating-Gas or SAG=0 if co-injection 

%% enter the loop for caculation in the next time step
while TIME<=TIME_TOTAL

    if SAG==1
        CYCLE=0.3; %SAG cycle throughput
        fg=2/3; %gas fraction per cycle
    for JJJJ=1:(round(TIME_TOTAL/CYCLE)+1) %loop for SAG process
        AAAA=JJJJ*CYCLE;    
        if TIME>=(AAAA-CYCLE) & TIME<(AAAA-0.0)
            if TIME>=(AAAA-CYCLE) & TIME<(AAAA-CYCLE*fg)
                f_g_BC=0;    
            else
                f_g_BC=1;
            end
        end
    end  
     u_w(1)=1-f_g_BC;
    end
      
    S_w_old=1-S_g; %water saturation in the current time step for calculation of surfactant concentration
    S_g_old=S_g; %gas saturation in the current time step for calculation of surfactant concentration
    u_w_old=u_w; %water fractional flow in the current time step for calculation of surfactant concentration
    u_g_old=u_g; %gas fractional flow in the current time step for calculation of surfactant concentration
    
    %% solve Tridiagonal System by LU decomposition and backward substitution  START%%
    UPPER=ones(1,NX-1);   
    for uu=1:(NX-1)
        UPPER(uu)=lamda(uu);
    end

    MIDDLE=ones(1,NX);
    for dd=2:(NX-1)
        MIDDLE(dd)= -(lamda(dd-1)+lamda(dd));
    end    
    MIDDLE(1)=-lamda(1);
    MIDDLE(NX)=-(lamda(NX-1)+2*lamda(NX));
        
    if pot(NX)<0 %prevent reverse flow of invading phase       
        kri(NX)=0;         
        krd(NX)=1;         
        lamda(NX)=kri(NX)+krd(NX)/M;         
        d(NX)=-(lamda(NX-1)+2*lamda(NX));     
    end
        
    LOWER=ones(1,NX-1);      
    for ll=1:(NX-1)      
        LOWER(ll)=lamda(ll);     
    end   
    [u,d,l]=decomt(UPPER,MIDDLE,LOWER); %find LU decomposition Tridiagonal Matrix by "decomt.m"
     
    b=zeros(1,NX);  
    for k=2:(NX-1)  
        b(k)=1/M*(krd(k)*(pc(k+1)-pc(k))-krd(k-1)*(pc(k)-pc(k-1)));  
    end    
    b(1)=-1/NX+1/M*(krd(1)*(pc(2)-pc(1))); 
    b(NX)=1/M*(krd(NX-1)*pc(NX-1)-(krd(NX-1)+2*krd(NX))*pc(NX));
 
    pot=solvet(u,d,l,b); %find flow potential at TIME=TIME(old)+deltat with backward substitution by "solvet.m"
  
    %% solve Tridiagonal System by LU decomposition and backward substitution  END %%
        
   S_g(1)=S_g(1)+deltat/(deltax)^2*(kri(1)*(pot(2)-pot(1))+deltax*f_g_BC);  %BC at x=0; find saturation at TIME=TIME(old)+deltat
   for i=2:(NX-1)
   S_g(i)=S_g(i)+deltat/(deltax)^2*(kri(i)*(pot(i+1)-pot(i))-kri(i-1)*(pot(i)-pot(i-1)));
   end
   S_g(NX)=S_g(NX)+deltat/(deltax)^2*(kri(NX)*2*(0-pot(NX))-kri(NX-1)*(pot(NX)-pot(NX-1)));  %BC at x=1
   
   S_nor=(S_g-S_gr)/(1-S_gr-S_wc); % normalized gas saturation for capillary pressure calculation

       for kk=1:NX  %%prevent the generation of complex number in krg and/or krw by round-off errors
       if S_nor(kk)<0
           S_nor(kk)=0;
       elseif S_nor(kk)>1
           S_nor(kk)=1;
       end
       end
       
    S_g=S_nor*(1-S_gr-S_wc)+S_gr;
    S_w=1-S_g;
    Fwater=0.5+(atan(epdry*(S_w-fmdry)))/pi;
    krd=(1-S_nor).^nw; %find aqueous relative permeability at TIME=TIME(old)+deltat
          
    for m=1:NX %capillary pressure model using Leveret-J function for geometrically similar systems
            if S_nor(m)<=1e-3    
                pc(m)=Npc*C1./exp(C2./log(1e-3))*(S_nor(m)/1e-3);   
            else if S_nor(m)>0.999               
                    pc(m)=Npc*C1./exp(C2./log(0.999));      
                else                    
                    pc(m)=Npc*C1./exp(C2./log(S_nor(m)));     
                end                
            end          
    end
    
     u_w(1)=1-f_g_BC; %In-flow BC
     u_w(NX+1)=-krd(NX)/M*2.*(0-pot(NX)+pc(NX))/deltax; 
            
    for iii=1:(NX-1) %find water fractional flow at grid block wall at TIME=TIME(old)+deltat
        u_w(iii+1)=-krd(iii)/M.*((pot(iii+1)-pc(iii+1))-(pot(iii)-pc(iii)))/deltax;
    end
    
    u_g=1-u_w;
%if TIME>=1 
%   Cw_inj=0; %%FINITE SLUG from TIME=0 to TIME=1%%
%end  

    Csg=Psgw*Csw; %surfactant concentration in gas phase
    for ii=1:NX           
        Css(ii)=rou_s*a_s*T_max*alpha*Csw(ii)/(1+alpha*Csw(ii)); 
        para1= alpha*(S_w(ii)+Psgw*S_g(ii));
       
         if ii==1                        
            N_temp=(Csw(ii)*S_w_old(ii)+Csg(ii)*S_g_old(ii)+(1-poro)/poro*Css(ii))+deltat*(1/Pe_w...
                *(S_w_old(ii)*(Csw(ii+1)-Csw(ii))-(1-f_g_BC)*2*(Csw(ii)-Cw_inj))/(deltax^2)+...
               1/Pe_g*(S_g_old(ii)*(Csg(ii+1)-Csg(ii))-f_g_BC*2*(Csg(ii)-Cg_inj))/(deltax^2)...
               -(Csw(ii)*u_w_old(ii+1)+Csg(ii)*u_g_old(ii+1)-Cw_inj*u_w_old(ii)-Cg_inj*u_g_old(ii))/deltax);        
        elseif ii==NX
             N_temp=(Csw(ii)*S_w_old(ii)+Csg(ii)*S_g_old(ii)+(1-poro)/poro*Css(ii))+deltat*(1/Pe_w...
                *(S_w_old(ii)*(Csw(ii)-Csw(ii))-S_w_old(ii-1)*(Csw(ii)-Csw(ii-1)))/(deltax^2)+...
               1/Pe_g*(S_g_old(ii)*(Csg(ii)-Csg(ii))-S_g_old(ii-1)*(Csg(ii)-Csg(ii-1)))/(deltax^2)...
               -(Csw(ii)*u_w_old(ii+1)+Csg(ii)*u_g_old(ii+1)-Csw(ii-1)*u_w_old(ii)-Csg(ii-1)*u_g_old(ii))/deltax);     
        else
            N_temp=(Csw(ii)*S_w_old(ii)+Csg(ii)*S_g_old(ii)+(1-poro)/poro*Css(ii))+deltat*(1/Pe_w...
                *(S_w_old(ii)*(Csw(ii+1)-Csw(ii))-S_w_old(ii-1)*(Csw(ii)-Csw(ii-1)))/(deltax^2)+...
               1/Pe_g*(S_g_old(ii)*(Csg(ii+1)-Csg(ii))-S_g_old(ii-1)*(Csg(ii)-Csg(ii-1)))/(deltax^2)...
               -(Csw(ii)*u_w_old(ii+1)+Csg(ii)*u_g_old(ii+1)-Csw(ii-1)*u_w_old(ii)-Csg(ii-1)*u_g_old(ii))/deltax);
        end
        para2=S_w(ii)+Psgw*S_g(ii)+(1-poro)/poro*rou_s*a_s*T_max*alpha-alpha*N_temp;
        para3=-N_temp;
        
        if para1==0          
             Csw_NEW(ii)=-para3/para2;        
            else
             Csw_NEW(ii)=(-para2+sqrt(para2^2-4*para1*para3))/(2*para1); 
             %find surfactant concentration at TIME=TIME(old)+deltat
             %a new temporary variable Csw_NEW is used to avoid it to replace the value C_w(ii) in the old time step of the iteration
        end
              
        if Csw_NEW(ii)>Cw_inj %round-off error can cause serious oscillation        
           Csw_NEW(ii)=Cw_inj;   
        elseif Csw_NEW(ii)<0             
                Csw_NEW(ii)=0;            
        end
    end         
        Csw=Csw_NEW; %surfactant concentration at TIME=TIME(old)+deltat
 
%if TIME>=1
%   NaCl_w_0=0.0; %FINITE SLUG from TIME=0 to TIME=1.0
%end  

    for jj=1:NX
        NaCl_s(jj)=rou_s*a_s*Tr_max*alpha*NaCl_w(jj)/(1+alpha*Cw_inj*NaCl_w(jj)); 
        para1= alpha*Cw_inj*S_w(jj)/deltat;
         if jj==1                        
            para2= (S_w(jj)+(1-poro)/poro*rou_s*a_s*Tr_max*alpha-(NaCl_w(jj)*S_w_old(jj)...
                +(1-poro)/poro*NaCl_s(jj))*alpha*Cw_inj)/deltat...
                + alpha*Cw_inj*((NaCl_w(jj)*u_w_old(jj+1)-NaCl_w_0*u_w_old(jj))/deltax...
                -1/Pe_w*(S_w_old(jj)*(NaCl_w(jj+1)-NaCl_w(jj))-(1-f_g_BC)*2*(NaCl_w(jj)-NaCl_w_0))/(deltax^2));
            para3= -(NaCl_w(jj)*S_w_old(jj)+(1-poro)/poro*NaCl_s(jj))/deltat...
                + (NaCl_w(jj)*u_w_old(jj+1)-NaCl_w_0*u_w_old(jj))/deltax...
                -1/Pe_w*(S_w_old(jj)*(NaCl_w(jj+1)-NaCl_w(jj))-(1-f_g_BC)*2*(NaCl_w(jj)-NaCl_w_0))/(deltax^2);
            elseif jj==NX
            para2= (S_w(jj)+(1-poro)/poro*rou_s*a_s*Tr_max*alpha-(NaCl_w(jj)*S_w_old(jj)...
                +(1-poro)/poro*NaCl_s(jj))*alpha*Cw_inj)/deltat...
                + alpha*Cw_inj*((NaCl_w(jj)*u_w_old(jj+1)-NaCl_w(jj-1)*u_w_old(jj))/deltax...
                -1/Pe_w*(S_w_old(jj)*(NaCl_w(jj)-NaCl_w(jj))-S_w_old(jj-1)*(NaCl_w(jj)-NaCl_w(jj-1)))/(deltax^2));
            para3= -(NaCl_w(jj)*S_w_old(jj)+(1-poro)/poro*NaCl_s(jj))/deltat...
                + (NaCl_w(jj)*u_w_old(jj+1)-NaCl_w(jj-1)*u_w_old(jj))/deltax...
                -1/Pe_w*(S_w_old(jj)*(NaCl_w(jj)-NaCl_w(jj))-S_w_old(jj-1)*(NaCl_w(jj)-NaCl_w(jj-1)))/(deltax^2);     
         else            
            para2= (S_w(jj)+(1-poro)/poro*rou_s*a_s*Tr_max*alpha-(NaCl_w(jj)*S_w_old(jj)...
                +(1-poro)/poro*NaCl_s(jj))*alpha*Cw_inj)/deltat...
                + alpha*Cw_inj*((NaCl_w(jj)*u_w_old(jj+1)-NaCl_w(jj-1)*u_w_old(jj))/deltax...
                 -1/Pe_w*(S_w_old(jj)*(NaCl_w(jj+1)-NaCl_w(jj))-S_w_old(jj-1)*(NaCl_w(jj)-NaCl_w(jj-1)))/(deltax^2));
            para3= -(NaCl_w(jj)*S_w_old(jj)+(1-poro)/poro*NaCl_s(jj))/deltat...
                + (NaCl_w(jj)*u_w_old(jj+1)-NaCl_w(jj-1)*u_w_old(jj))/deltax...
                -1/Pe_w*(S_w_old(jj)*(NaCl_w(jj+1)-NaCl_w(jj))-S_w_old(jj-1)*(NaCl_w(jj)-NaCl_w(jj-1)))/(deltax^2);
        end
        
        if para1==0          
            NaCl_w_NEW(jj)=-para3/para2;        
            else
            NaCl_w_NEW(jj)=(-para2+sqrt(para2^2-4*para1*para3))/(2*para1); %find NaCl concentration at TIME=TIME(old)+deltat
        end
              
        if NaCl_w_NEW(jj)>1 %round-off error can cause serious oscillation
            NaCl_w_NEW(jj)=1;     
        elseif NaCl_w_NEW(jj)<0             
                NaCl_w_NEW(jj)=0;                     
        end
    end   
         NaCl_w=NaCl_w_NEW;
    
         for mm=1:NX 
         if Csw(mm)<fmsurf     
             Fsurf(mm)=(Csw(mm)/fmsurf).^epsurf;
         else
             Fsurf(mm)=1;
         end
         end
         
    FM=1./(1+fmmob*Fwater.*Fsurf);
    kri=FM.*(S_nor.^ng); %find gasous relative permeability at TIME=TIME(old)+deltat
    lamda=kri+krd/M; %find lamda at TIME=TIME(old)+deltat
    
    COUNTER=COUNTER+1;
    TIME=TIME+deltat;
    TIME_HIS(COUNTER)=TIME;
    
    if (TIME-0.5)<=1e-6 & (TIME-0.5)>=-1e-6
        TIME1=TIME
        S_time1=S_g;
        C_w_time1=Csw/Cw_inj;
        NaCl_w_time1=NaCl_w;
        pot_time1=pot;
        pc_time1=pc;
        pwd_time1=pot_time1-pc_time1;
    end
    
    if (TIME-1.0)<=1e-6 & (TIME-1.0)>=-1e-6
        TIME2=TIME
        S_time2=S_g;
        C_w_time2=Csw/Cw_inj;
        NaCl_w_time2=NaCl_w;
        pot_time2=pot;
        pc_time2=pc;
        pwd_time2=pot_time2-pc_time2;
    end
    
    if (TIME-5.0)<=1e-6 & (TIME-5.0)>=-1e-6  
        TIME3=TIME
        S_time3=S_g;                  
        C_w_time3=Csw/Cw_inj;
        NaCl_w_time3=NaCl_w;  
        pot_time3=pot;                                   
        pc_time3=pc;
        pwd_time3=pot_time3-pc_time3;
    end       
   
    f_g_HIS(COUNTER)=kri(NX).*pot(NX)./(kri(NX).*pot(NX)+krd(NX).*(pot(NX)-pc(NX))/M); 
    f_total_HIS(COUNTER)=2*(kri(NX).*pot(NX)+krd(NX).*(pot(NX)-pc(NX))/M)/deltax;   
    f_surf_HIS(COUNTER)=Csw(NX)/Cw_inj; %surfactant effluent history in aqueous phase (g/L)
    f_NaCl_HIS(COUNTER)=NaCl_w(NX);
    p_drop_HIS(COUNTER)=pot(1)*NX/(NX-0.5);
    VIS_FOAM(COUNTER)=p_drop_HIS(COUNTER)*visg/krg0; %foam apparent viscosity
    
end

x1=deltax/2:deltax:(1-deltax/2);

prof1=[S_time1' C_w_time1' NaCl_w_time1' x1'];
prof2=[S_time2' C_w_time2' NaCl_w_time2'];
prof3=[S_time3' C_w_time3' NaCl_w_time3'];
prof4=[TIME_HIS' f_g_HIS' f_surf_HIS' f_NaCl_HIS' p_drop_HIS' VIS_FOAM'];
prof5=[pot_time1' pot_time2' pot_time3'];
prof6=[pwd_time1' pwd_time2' pwd_time3'];


save prof1.dat prof1 /ascii
save prof2.dat prof2 /ascii
save prof3.dat prof3 /ascii
save prof4.dat prof4 /ascii
save prof5.dat prof5 /ascii
save prof6.dat prof6 /ascii

comp_time = cputime