Mixture critical point calculator based on Heidemann & Khalil/ Michelsen & Heidemann.

critpt_eos.matlab

function [T,P,v] = critpt_eos(z, Tci, Pci, w, kij, lij, varargin)
% CRITPT_EOS Calculates the critical temperature, pressure and volume of
% a multicomponent mixture using the method developed by Heidemann and
% Khalil (1980) and Michelsen and Heidemann (1981). Can use both
% Peng-Robinson and Soave-Redlich-Kwong EOS. Two parameter van der Waals 
% mixing rules are used. Initial critical property estimates are generated 
% by the program unless specified by the user (see below).
%
% Inputs:
% n   - vector of molar compositions of components 
% Tci - Critical temperatures of pure components (K)
% Pci - Critical pressures of pure components (bar)
% w   - Acentric factors of pure components 
% 
% Optional:
% kij - Binary interaction parameters of n species in vector form
%       ex. kij = [k12 k13 k14 k15 k23 k24 k25 k34 k35 k45]
% lij - Binary interaction parameters of n species in vector form
%       ex. lij = [l12 l13 l14 l15 l23 l24 l25 l34 l35 l45]
% 
% Optional Parameters:
% eos - Flag indicating which EOS to use (1 = PR, 2 = SRK; default = PR)
% Ti  - Starting temperature for Tc iteration (K)
% vi  - Starting molar volume for v iteration (cm^3/mol)
% ver - Verbose output (default = 0, 1 = on)
% 
% Outputs:
% T  - Critical temperature of mixture (K)
% P  - Critical pressure of mixture (bar) 
% v  - Critical volume of mixture (cm^3/mol)
%
% Example Usage:
% -- CO2 and ethanol -- 
% Tci = [304.2, 516.4]; % (K)
% Pci = [73.82, 63.84]; % (bar)
% w   = [0.288, 0.637]; 
% kij = 0.06; % Secuianu et al. (2008)
% n = ...
% [T, P, v] = CRITPT_EOS(n, Tci, Pci, w, kij, 'ver', 1);

% Copyright Hythem Sidky (c) 2012
% 
% This work is licensed under a Creative Commons 
% Attribution 3.0 Unported (CC BY 3.0) Unported License.
% http://creativecommons.org/licenses/by/3.0/
%
% Changes:
% 04/24/14 - Modified interaction parameter input to packed format.
% 07/07/12 - Corrected error in summation.
% 07/04/12 - Initial file.

% References: 
% 1 - Heidemann, Robert A., and Ahmed M. Khalil. "The Calculation of 
%     Critical Points." AIChE Journal 26.5 (1980): 769-79. 
% 2 - Michelsen, Michael L., and Robert A. Heidemann. "Calculation of 
%     Critical Points from Cubic Two-constant Equations of State." AIChE 
%     Journal 27.3 (1981): 521-23.

    % Get component count 
    n = length(z);
    
    % Parse user inputs 
    p = inputParser;
    % - Required 
    p.addRequired('n',  @(x)~isempty(x));
    p.addRequired('Tci',@(x)length(x) == n);
    p.addRequired('Pci',@(x)length(x) == n);
    p.addRequired('w',  @(x)length(x) == n);
    
    % - Optional
    p.addOptional('kij', 0, @(x)length(x) == factorial(n)/(2*factorial(n-2)));
    p.addOptional('lij', 0, @(x)length(x) == factorial(n)/(2*factorial(n-2)));
    % - Optional Parameters
    p.addParamValue('Ti', 0,@(x)x>0);
    p.addParamValue('vi', 0,@(x)x>0);
    p.addParamValue('eos',1,@(x)x==1 || x==2);
    p.addParamValue('ver',0,@(x)x==0 || x==1);
    
    % - Parse
    p.parse(z,Tci,Pci,w,kij,lij,varargin{:});
    z   = p.Results.n;
    Tci = p.Results.Tci;
    Pci = p.Results.Pci;
    w   = p.Results.w;
    kij = p.Results.kij;
    lij = p.Results.lij;
    eos = p.Results.eos;
    Ti  = p.Results.Ti;
    vi  = p.Results.vi;
    ver = p.Results.ver;
    
    % Expand packed interaction parameter vector 
    km = zeros(n);
    lm = zeros(n);
    for i=1:n
        for j=i:2
            if i==j
                km(i,j) = 0;
                lm(i,j) = 0;
            else
                km(i,j) = kij(j+(2*n-i)*(i-1)/2-i);
                lm(i,j) = lij(j+(2*n-i)*(i-1)/2-i);
                km(j,i) = km(i,j);
                lm(j,i) = lm(i,j);
            end
        end
    end
    
    % Override previous entries 
    kij = km;
    lij = lm;
    
    % If user supplied a mol frac of 0 for any component, just remove it.
    % This makes things easy when generating a critical curve in an
    % external call
    j = find(z == 0);
    n = length(j); % Update count
    for i=1:n
        z(j(i))   = [];
        Tci(j(i)) = [];
        Pci(j(i)) = [];
        w(j(i))   = [];
        kij(j(i),:) = [];
        kij(:,j(i)) = [];
        lij(j(i),:) = [];
        lij(:,j(i)) = [];
    end
    
    % Now let's check if the result of the above process leaves us with a
    % pure compound. If that's the case then just return supplied
    % properties. User doesn't supply v so return 0.
    if length(z) == 1
        if ver
            fprintf('Only a single compound supplied.');
            fprintf(' Returning pure properties\n');
        end
        T  = Tci; 
        P  = Pci;
        v = 0;
        return;
    end
    
    % Define R and normalize n
    R = 83.14472; % (bar cm3/mol-K)
    z = z./sum(z);
    
    % Initial guess of Tc, if user didn't supply one. We will evaluate 
    % initial v guess later because we need to calculate some thermo 
    % properties first.
    T = Ti;
    if ~T
        T = 1.5*sum(z.*Tci);
    end
    
    % Declare thermodynamics properties anonymous functions 
    props = @(n,T,v)calc_eos_props(n, T, v, Tci, Pci, w, kij, lij, eos);
    
    % Generate initial v guess by evaluating b. Below is not ideal but we
    % really only need b. It doesn't matter what v is (obviously b is
    % independent of it).
    v = vi;
    [~,~,~,~,~,b,~] = props(z, T, 0);
    if ~v
        v = 4*b;
    end

    % Display initial conditions 
    if ver
        fprintf('\nInitial estimates: Tc = %.3f K, v = %.3f (cm^3/mol)\n', ...
            T, v);
        fprintf('\nSolving for initial Tc where det(Q)=0\n');
    end
   
    % Get an initial evaluation on Tcrit then begin finding zero of cubic
    % form Cf which contains a nested routine for calc_T_crit
    f = @(T)calc_T_crit(z,T,v,props);
    T = fzero(f,T);
    
    % Display initial calculated Tc
    if ver
        fprintf('Initial Tc = %.3f K\n',T);
        fprintf('\nBegin soving both criticality conditions\n');
    end
   
    % Solve for critical volume
    f = @(v)calc_v_crit(z,T,v,props,ver);   
    v = fzero(f,v);
    
    % Re-evaluate critical volume for appropriate temperature
    [~,T,~] = f(v);
    
    % Return critical temperature, pressure and volume
    [~,D,~,~,a,b,~] = props(z, T, v);
    P = R*T/(v-b)-a/((v+D(1)*b)*(v+D(2)*b));
    
    % Display final results
    if ver
        fprintf('\nFinal Results:\n');
        fprintf('Tc = %.3f K, Pc = %.3f bar, v = %.3f (cm^3/mol)\n', ...
            T, P, v);
    end

% calc_v_crit - Calculates the cubic form Cf for a given T. Returns Cf,
% iterated temperature T from inside loop, and perturbation vector N.
function [Cf, T, N] = calc_v_crit(n, T, v, props, ver)
    
    % Total number of moles
    C  = length(n);
    nt = sum(n);    
    
    % Define R
    R = 83.14472; % (bar cm3/mol-K)
    
    % Evaluate Tcrit for given v 
    tcrit = @(T)calc_T_crit(n,T,v,props);
    T = fzero(tcrit,T);
    
    % Re-evaluate tcrit for det(Q), then get thermo properties
    [~, Q] = tcrit(T);
    [F, ~, A, B, a, b, aij] = props(n, T, v);

    % Calculate a non-zero perturbation vector in component mols (N) based
    % on Q*N=0 using SVD then normalize the vector.
    [U,~,~] = svd(Q);
    N = U(:,end);
    N = N'/norm(N);
    
    % Calculate additional critical point parameters (underscore indicates
    % a bar on variable in paper.
    N_ = sum(N);
    A_ = sum(N.*A);
    B_ = sum(N.*B);
    a_ = 0;
    for i=1:C
        for j=1:C
            a_ = a_ + N(i)*N(j)*aij(i,j);
        end
    end
    a_ = a_/a;
    
    % Calculate cubic form Cf. Split into terms for readability. Note in
    % [2], "RT" is multipled by first term. However in original paper [1]
    %  it is divided through. 
    t1 = -sum(N.^3./n.^2)+3*N_*(B_*F(1))^2+2*(B_*F(1))^3;
    t2 = 3*B_^2*(2*A_-B_)*(F(3)+F(6))-2*B_^3*F(4)-3*B_*a_*F(6);
    Cf = 1/nt^2*(t1)+a/(nt^2*b*R*T)*t2;
    Cf = Cf*((v-b)/(2*b))^2; % Dimensionless scaling term [1]
    
    % Display iteration
    if ver
        fprintf('Iteration T = %.3f K, v = %.3f (cm^3/mol)\n', T, v);
    end

% calc_T_crit - Calculates hessian matrix Q and its determinant detQ for a
% given temperature. The name of the function is misleading because it
% doesn't actually calculate Tc. It's simply used in it's calculation.
function [detQ, Q] = calc_T_crit(n, T, v, props)
    
    % Total number of moles
    C  = length(n);
    nt = sum(n);    
    
    % Define R
    R = 83.14472; % (bar cm3/mol-K)
    
    % Kronecker delta
    delta = eye(C);
    
    % Get thermodynamics properties 
    [F, ~, A, B, a, b, aij] = props(n, T, v);
    
    % Calculate matrix Q and determinant. 
    Q = zeros(C);
    for i=1:C
        for j=1:C
            % Split terms up for readability. Note in [2], "RT" is 
            % multipled by first term, however in original paper [1] it is 
            % divided through. This scaling helps with root finding.
            t1 = (1/nt)*(delta(i,j)/n(i)+(B(i)+B(j))*F(1)+B(i)*B(j)*F(1)^2);
            t2 = (B(i)*B(j)-A(i)*B(j)-A(j)*B(i))*F(6);
            t3 = a/(b*R*T*nt)*(B(i)*B(j)*F(3)-aij(i,j)/a*F(5)+t2);
            Q(i,j) = (t1+t3)*T/100; % Dimensionless scaling factor [1]
        end
    end
    detQ = det(Q);
    
% calc_eos_props - Calculates the thermodynamic properties including EOS
% parameters for the input components. Returns mixture properties and
% auxillary functions used in critical calculation.
function [F, D, A, B, a, b, aij] = calc_eos_props(n, T, v, Tci, Pci, w, kij, lij, eos) 
    
    % Total number of moles
    C  = length(n);
    nt = sum(n); 
    
    % Define reduced temperature and R
    Tri = T./Tci;
    R   = 83.14472; % (bar cm3/mol-K)
    
    % - EOS assignments. Switch between PR and SRK EOS depending on user
    %   input. These are condition independent parameters.
    if eos==1 % Peng-Robinson
        u0  = 2;
        w0  = -1;
        psi = 0.45724; 
        omg = 0.07780;
        m  = 0.37464+1.54226.*w-0.26992*w.^2;
    else % Soave-Redlich-Kwong
        u0  = 1;
        w0  = 0;
        psi = 0.42748;
        omg = 0.08664;
        m  = 0.480+1.574.*w-0.176*w.^2;
    end
    D1 = (u0+sqrt(u0^2-4*w0))/2;
    D2 = (u0-sqrt(u0^2-4*w0))/2;
    
    % - Calculate EOS pure component parameters
    ai = (R.*Tci).^2*psi./Pci.*(1+m.*(1-sqrt(Tri))).^2;
    bi = omg*R*Tci./Pci;
    
    % - Calculate EOS mixture parameters
    a   = 0;
    b   = 0;
    
    aij = zeros(C);
    bij = zeros(C);
    for i=1:C
        for j=1:C
            aij(i,j)=sqrt(ai(i)*ai(j))*(1-kij(i,j));
            bij(i,j)=sqrt(bi(i)*bi(j))*(1-lij(i,j));
            a = a + n(i)*n(j)*aij(i,j)/nt^2;
            b = b + n(i)*n(j)*bij(i,j)/nt^2;
        end
    end
    
    % Override b. comment this line to have geometric mixing for 
    % b as well.
    b  = sum(n.*bi/nt);
    
    B = bi./b;      % beta
    A = zeros(1,C); % alpha
    for i=1:C
        for j=1:C
            A(i) = A(i)+n(j)*aij(i,j);
        end
        A(i) = A(i)/a;
    end
        
    % Dimensionless volume and auxillary functions
    K  = v/b;
    F1 = 1/(K-1);
    F2 = 2/(D1-D2)*(D1/(K+D1)-D2/(K+D2));
    F3 = 1/(D1-D2)*((D1/(K+D1))^2-(D2/(K+D2))^2);
    F4 = 1/(D1-D2)*((D1/(K+D1))^3-(D2/(K+D2))^3);
    F5 = 2/(D1-D2)*log((K+D1)/(K+D2));
    F6 = F2 - F5;
    
    % Return matrices 
    F = [F1, F2, F3, F4, F5, F6];
    D = [D1, D2];             
        
        

Hi, Mr. Sidky, it's me again.
I found that the calculation result seems to be the same whether I use the only kij or both kij and lij in the vdW mixing rule, why is that?

%% CO2+propane
Tci = [304.1282, 369.89];     % (K)
Pci = [73.773, 42.512];         % (bar)
w   = [0.22394, 0.1521];
kij = 0.131; 
lij = 0.2;
n = [0.5, 0.5];

[Tcm, Pcm] = critpt_eos(n, Tci, Pci, w, kij, 'ver', 0)
[Tcm2, Pcm2] = critpt_eos(n, Tci, Pci, w, kij, lij, 'ver', 0)

Tcm =

  333.8199


Pcm =

   63.3009


Tcm2 =

  333.8199


Pcm2 =

   63.3009

@RayleighSun If you would like to consider the lij, simply comment on line No. 339.

@RayleighSun If you would like to consider the lij, simply comment on line No. 339.

Thanks, it works!