Joule-Thomson inversion curve with the Peng-Robinson EOS

Computes the Joule-Thomson coefficient inversion curve by using the Peng-Robinson equation of state.
54 Downloads
Updated 29 Sep 2021

View License

This is a bare-bones method to plot the Joule-Thomson (coefficient) inversion curve based on Eq. (2) in:
Colina, C.M., & Olivera-Fuentes, C. (1998). Prediction of the Joule–Thomson inversion curve of air from cubic equations of state. Cryogenics, 38(7), 721–728. doi:10.1016/s0011-2275(98)00036-8.
The Joule-Thomson inversion curve is plotted with the (pr, Tr) points for which µJT = 0. And, according to Eq. (2), µJT is null if and only if (∂z/∂T) = 0 at constant pressure. Therefore, to the determine the inversion curve one must find the reduced temperatures Tr for fixed, reduced pressures pr that satisfy (∂z/∂T) = 0. The Peng-Robinson equation of state is used to compute the compressibility factor, z (Ref: https://doi.org/10.1021/i160057a011).
Run the script as JTCInvert_PR(number), where number is the component ID: 1 for CO2, 2 for Nitrogen; 3 for H2O, 4 for CH4, 5 for C2H6, 6 for C3H8 and 7 for SF6. Example: JTCInvert_PR(4) plots the inversion curve of the methane.
The script can easily be adapted to include more components. Critical pressures, temperatures and acentric factors for several substances can be found in:
https://webbook.nist.gov/chemistry/fluid/
function JTCInvert_PR(component)
%
% A bare-bones method to plot the Joule-Thomson (coefficient) inversion
% curve based on Eq. (2) in:
% Colina, C. M., & Olivera-Fuentes, C. (1998). Prediction of the Joule–
% Thomson inversion curve of air from cubic equations of state.
% Cryogenics, 38(7), 721–728. doi:10.1016/s0011-2275(98)00036-8.
% That is, finding the reduced temperatures Tr for fixed, reduced pressures
% pr, that satisfy (∂z/∂T) = 0. The Peng-Robinson equation of state is
% used (Ref: https://doi.org/10.1021/i160057a011).
% Run the script as JTCInvert_PR(number), where number is the component
% ID: 1 for CO2, 2 for Nitrogen; 3 for H2O, 4 for CH4, 5 for C2H6, 6 for
% C3H8 and 7 for SF6. Example: JTCInvert_PR(4) plots the inversion curve
% of the methane.
% The script can easily be adapted to include more components. Critical
% pressures, temperatures and acentric factors for several substances
% can be found in:
% https://webbook.nist.gov/chemistry/fluid/
% August, 2021.
format long;
% Data is a 7x4 matrix whose columns are:
% Component ID, critical temperature, critical pressure and acentric
% factor.
Data = [1, 304.13, 73.773, 0.22394; ...
2, 126.19, 33.958, 0.03700; ...
3, 647.14, 220.640, 0.3443; ...
4, 190.56, 45.992, 0.01100; ...
5, 305.33, 48.718, 0.09930; ...
6, 369.83, 42.477, 0.15240; ...
7, 318.73, 37.546, 0.21000];
% Those are the components this script is currently capable of
% processing:
Comp = ["CO2"; "N2"; "H2O"; "CH4"; "C2H6"; "C3H8"; "SF6"];
nome = Comp(component, 1);
Tc = Data(component, 2);
pc = Data(component, 3);
omega = Data(component, 4);
np = 45;
nT = 40;
% Upper limits for Tr (5.5) and pr (15.0) should be enough.
Tr = linspace(0.5, 5.5, nT);
pr = linspace(0.1, 15.0, np);
W = zeros(np*nT+1,3);
% Tinv and Pinv must not be (0, 0, 0). This is to prevent warnings.
Tinv = zeros(1,3);
pinv = zeros(1,3);
reps = 1.490116e-08;
k = 1;
for ip=1:np
for iT=1:nT
deltaT = Tr(iT)*reps;
% Compute z(pr, Tr+ΔTr), z(pr, Tr-ΔTr) and the correspondent
% central difference (numerical derivative, first order) for
% a constant pr.
z1 = zpengrobinson(pr(ip), Tr(iT)-deltaT, pc, Tc, omega);
z2 = zpengrobinson(pr(ip), Tr(iT)+deltaT, pc, Tc, omega);
dzdT = (z2 - z1)/(2.0*Tc*deltaT);
% Store pr, Tr and dzdT values in W.
W(k,1) = pr(1,ip);
W(k,2) = Tr(1,iT);
W(k,3) = dzdT;
k = k + 1;
end
end
j = 1;
for i=2:np*nT-3
%Change of sign between two consecutive dzdT values indicate the
%presence of a root, that is, ∂z/∂Tr = 0 at constant pr.
if W(i,3)*W(i+1,3) <= 0
%Interpolate to find the zero using MATLAB's interp1.
val = interp1([W(i,3); W(i+1,3)], [W(i,2);W(i+1,2)], 0.0);
pinv(j) = W(i,1);
Tinv(j) = val;
j = j + 1;
end
end
str1 = strcat("Joule-Thomson inversion curve (Peng-Robinson) for ",nome);
figure('Name','JTC Inversion curve','Position',[400,150,750,550])
plot(pinv, Tinv, "s", linewidth=1, markersize=3, color="#dd0000");
grid on;
xlabel(texlabel("p_r"), "FontName","Times", "FontAngle","Oblique", fontsize=16);
ylabel(texlabel("T_r"), "FontName","Times", "FontAngle","Oblique", fontsize=16);
title(str1, "FontName","Times", "FontAngle","Oblique", fontsize=12);
end
function result = zpengrobinson(pr_, Tr_, pc, Tc, omega)
%
% Peng, D.-Y., & Robinson, D. B. (1976). A New Two-Constant Equation of
% State. Industrial & Engineering Chemistry Fundamentals, 15(1), 59–64.
% doi:10.1021/i160057a011
%
% Parameters
% ----------
% pr_ : Reduced pressure.
% Tr_ : Reduced temperature.
% pc : Critical pressure (bar).
% Tc : Critical temperature (K).
% omega : Acentric factor
%
% Returns
% -------
% z : Compressibility factor.
R = 8.3144721345e-02; % bar*m3/(kmol*K)
TT = Tr_*Tc;
pp = pr_*pc;
m = 3.74640e-01 + 1.54226e0*omega - 2.69920e-01*omega^2;
ac = 4.57236e-01*(R*Tc)^2/pc;
a = ac*(1.0 + m*(1.0 - sqrt(TT/Tc)))^2;
b = 7.77961e-02*R*Tc/pc;
AA = a*pp/(R*TT)^2;
BB = b*pp/(R*TT);
z = 9.5e-01;
erro = 1.0e+03;
c0 = -(AA*BB - BB^2 - BB^3);
c1 = (AA - 2.0*BB - 3.0*BB^2);
c2 = -(1.0 - BB);
while erro >= 1.0e-08
f = z^3 + c2*z^2 + c1*z + c0;
df = 3.0e0*z^2 + 2.0e0*c2*z + c1;
zf = z - f/df;
erro = abs((zf - z)/z);
z = zf;
end
result = z;
end

Cite As

Fausto Arinos de Almeida Barbuto (2024). Joule-Thomson inversion curve with the Peng-Robinson EOS (https://www.mathworks.com/matlabcentral/fileexchange/98234-joule-thomson-inversion-curve-with-the-peng-robinson-eos), MATLAB Central File Exchange. Retrieved .

MATLAB Release Compatibility
Created with R2021a
Compatible with any release
Platform Compatibility
Windows macOS Linux

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
Version Published Release Notes
1.0.2

Constants multiplying "ac" and "b" were slightly modified.

1.0.1

Minor correction in one of the script's comment lines.

1.0.0