Code covered by the BSD License

# DAK Equation of State

### Chigozie Aniemena (view profile)

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

```