Code covered by the BSD License  

Highlights from
DAK Equation of State

DAK Equation of State

by

 

05 Jul 2012 (Updated )

Tabular and graphical implementation of the Dranchuk-Abbou Kassem equation of state

DAKtable(sg,minP,maxP,Pstep,T)
% AUTHOR: CHIGOZIE ANIEMENA

function DAKtable(sg,minP,maxP,Pstep,T)


% Calculates the compressibility factor of natural gases for a range of pressures (Psia) 'minP' to
% 'maxP' in steps 'Pstep' (for a given temperature T(degree Farenheit) and specific gravity sg) using the Dranchuk-Abbou Kassem
% Equation of state. The range of validity is 0.2<Ppr<30 and 1.0<Tpr<3.0 (Craft
% &Hawkins, 1991)

%To calculate for a single pressure point P, set minP = maxP = Pstep
format short g
i = 1;
for P = minP:Pstep:maxP % 
    % Estimation of Psuedocritical pressure and temperature based on
    % suttons correlations (Craft & Hawkins, 1991)
    
    Ppc = 756.8 - (131*sg) - (3.6*(sg^2));
    Tpc = (169.2 + (349.5*sg) - (74*(sg^2)));
    
    %Calculation of Pseudo reduced pressure and temperature
    Ppr = P / Ppc;
    Tpr = (T+459.67)/ Tpc;
    
    % Correlation Constants, c4 has to be updated within the root finding
    % algorithm
    A1 = 0.3265;
    A2 = -1.0700;
    A3 = -0.5339;
    A4 = 0.01569;
    A5 = -0.05165;
    A6 = 0.5475;
    A7 = -0.7361;
    A8 = 0.1844;
    A9 = 0.1056;
    A10 = 0.6134;
    A11 = 0.7210;
    c1 = A1 + (A2/Tpr) + (A3/(Tpr^3))+ (A4/(Tpr^4))+ (A5/(Tpr^5));
    c2 = A6 + (A7/Tpr) + (A8/(Tpr^2));
    c3 = A9*((A7/Tpr) + (A8/(Tpr^2)));

   % Implementation of bisection root finding method (The Newton-Raphson method becomes unstable and slower 
   % at low temperature and higher specific gravities, the bisection method
   % appears more robust for all cases in this specific application)
    brack1 = 1e-2;
    brack2 = 4;
    itr= 0;
    while 2<3 
        zguess = (brack1 + brack2) / 2;
        Rrguess = (0.27*Ppr) /(zguess*Tpr);
        c4 = (A10)*(1+(A11*(Rrguess^2)))*((Rrguess^2)/(Tpr^3))*(exp(-A11*(Rrguess^2)));
        syms z
        Rr = (0.27*Ppr) /(z*Tpr);
        f = z+ (c3*(Rr^5)) - (c2*(Rr^2)) - (c1*(Rr^1)) - c4 - 1;

        fzguess = subs(f,z,zguess);
        
        convergence  = abs(fzguess-0);
        conv(itr+1) = convergence;
        disp_convergence = conv'; % remove ";" to display convergence progression for each pressure point
        
        if convergence <=10^-4
            itr = itr + 1;
            break
        end

        if fzguess < 0
            brack1 = zguess;
        else
            brack2 = zguess;
        end
        
      itr = itr + 1;
    
    end
    Z(i) = zguess;
    Pressure(i) = P;
    I__________P__________Z_____________I = [Pressure' Z']
    i=i+1;
    
%Uncomment following comments to view convergence progression for each
%pressure point



%     plot(disp_convergence)
%     xlabel('Convergence','fontsize',12)
%     ylabel('Iteration','fontsize',12)
%     hold on
%     pause

    clear conv
end

% hold off

figure;plot(Pressure,Z, 'r*','MarkerSize',6)
xlabel('Pressure','fontsize',12);
ylabel('Z','fontsize',12)


end

Contact us