Joule-Thomson inversion curve with the Peng-Robinson EOS

Computes the Joule-Thomson coefficient inversion curve by using the Peng-Robinson equation of state.

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