4 replies to this topic

[suman93]

Hello Everyone,

I am trying to calculate critical Temperature and Pressure of binary and ternary mixtures using equation of state. In particular, I am using Peng-Robinson Equation of State. I combined this with Van-der Waals critical criterion to get my result, but am way off experimental data.

 

Vander Waals criteria states: dP/dv = 0 and d2P/dv2=0

 

My steps:

1. I solve the two partial derivative equations mentioned above to end up with a cubic equation in 'v' (molar volume). Solving this I get the critical molar volume

2. I put this back in either of the partial derivative conditions and calculate Temperature using iteration.

3. Now I know Tc and Vc. I put these values back in the equation of state to get my pressure.

 

Problem:

My Tc values are still within range or close to experimental data but the Pc values are way off.

 

Attached Plot:

Attached plot is of a binary mixture of methane(C1) and propane(C3). Tc and Pc plots are attached vs mole fraction of C1

 

[Art Montemayor]

You state that you "trying to calculate critical Temperature and Pressure of binary and ternary mixtures using an equation of state".  However, you have a problem:  you are "way off experimental data".

  

What is your question or query?   It can't be that your calculations are wrong or not the correct ones because you haven't furnished your calculations.  If your question is: "why are my calculations not giving correct results?", then you must furnish your calculations - assuming that you have already had them proofed and checked for accuracy and applicability.  Our Forum members can't do much without a copy of your work in order to determine if it is the correct method and if it is applied correctly.

[breizh]

http://www.scielo.br...322000000400040

 

https://www3.nd.edu/...kst/bas2000.pdf

 

Hi ,

you may consider the links above .

Breizh

Edited by breizh, 09 August 2018 - 07:05 AM.

[suman93]

Hi,

Thank you so much for the response.

Breizh:  The links you have given use  the Gibb's method of finding the critical T and P.  I am not following that method.

The goal is to try and use the inflection point property and find the critical points. Which states dP/dv = 0 and d²P/dv² = 0 [these are partial derivatives of pressure with respect to molar volume]

I am giving my method here: Will be grateful to anyone who can point put a bug and show me the correct way

1. I start with Peng Robinson EOS: 

      P = { RT/(v- }  -   {  a / (v² + 2vb - b²) }

2. I find expressions for dP/dv = 0 and d²P / dv²= 0.

 

3. I equate them to end up with a cubic equation:  v³ - (3b)v² - (3b²)v - (3b³) = 0
    This has 3 roots but one is real and rest two are complex conjugates.
    The real root gives me the critical volume

4. I then put back this value of 'v' in the equation of   dP/dv = 0.   T is the only unknown in here now. I solve it to get T.

5. Now I put this T back in equation of state to get pressure.

I am giving my code below.

 

% Getting Pseudo Critical properties
% Using Peng-Robinson Equation of State
% Binary mixture: n-methane + propane
% Author: Suman Chakraborty

clear, clc;

%% User Inputs
% Universal gas constant
% R = 8.31446 J/mol K; 83.144 cm^3 bar/mol K
R = 83.144;
% mole fraction liquid phase
x_C1 = 0:0.001:1;
x_C3 = 1 - x_C1;
% size of mole fraction array
N = size(x_C1);
N = N(1,2);

%% methane Data
% Critical Temperature (K)
Tc_C1 = 190.6;
% Critical Pressure (bar)
Pc_C1 = 46.1;
% Molecular weight (gm/mole)
MW_C1 = 16.0425;
% Acentric Factor (w)
w_C1 = 0.011; %0.010

%% propane Data
% Critical Temperature (K)
Tc_C3 = 369.9;
% Critical Pressure (bar)
Pc_C3 = 42.5;
% Molecular weight (gm/mole)
MW_C3 = 44.0956;
% Acentric Factor (w) 
w_C3 = 0.153;   %0.152

%% methane Pure component calculation
% a(Tc)  [cm^6 * bar/mol^2]
a_Tc_C1 = (0.45724 * R * R * Tc_C1 * Tc_C1)/Pc_C1;
% b(Tc)  [cm^3 / mole]
b_C1 = (0.0778 * R * Tc_C1) / Pc_C1;
% kappa
kappa_C1 = 0.37464 + (1.54226*w_C1) - (0.26992*w_C1*w_C1);

%% propane Pure component calculation
% a(Tc)  [cm^6 * bar/mol^2]
a_Tc_C3 = (0.45724 * R * R * Tc_C3 * Tc_C3)/Pc_C3;
% b(Tc)  [cm^3 / mole]
b_C3 = (0.0778 * R * Tc_C3) / Pc_C3;
% kappa
kappa_C3 = 0.37464 + (1.54226*w_C3) - (0.26992*w_C3*w_C3);

%% Main Loop Starts
for i=1:N
    % Application of mixing rules
    % Repulsion parameter of mixture - [process 1]
    %bm = (x_C1(i)*b_C1) + (x_C3(i)*b_C3);
    
    % Repulsion parameter of mixture - [process 2]
    term1 = (x_C1(i)*(b_C1^(1/3))) + (x_C3(i)*(b_C3^(1/3)));
    bm = term1^3;
    
    % Solving cubic equation
    a3 = 1; % coeff of v^3
    a2 = -(3*bm); % coeff of v^2
    a1 = -(3*bm*bm); % coeff of v
    a0 = -3*(bm^3);  % constant
    coeff = [a3 a2 a1 a0];

    % Z is compressibility factor
    soln_v = roots(coeff);
    % Roots real and >0
    v_real = [];
    for j=1:3
        %>B to get real fugacity, since fhi has a term log(Z-B) 
        if isreal(soln_v(j)) && soln_v(j) > 0  
          temp = real(soln_v(j));
          v_real = [v_real temp];
        end    
    end
    
    index = length(v_real);
    if index >= 1
       % Vapor Phase
       v_final = max(v_real);
    end
    
    % Guess Temperature [K]
    T_guess = 2000;

    % Convergence variable
    err = 1.0;
    eps = 0.01;
    
     while true
        % Reduced Temperature
        Tr_C1 = T_guess / Tc_C1;
        % Reduced Temperature
        Tr_C3 = T_guess / Tc_C3;
        % a(T) of C1 [cm^6 * bar/mol^2]
        term1 = 1 - (Tr_C1^0.5);
        term2 = 1 + (kappa_C1*term1);
        a_T_C1 = a_Tc_C1 * term2 * term2;
        % a(T) of C3 [cm^6 * bar/mol^2]
        term1 = 1 - (Tr_C3^0.5);
        term2 = 1 + (kappa_C3*term1);
        a_T_C3 = a_Tc_C3 * term2 * term2;
        % Attraction parameter of mixture
        % Assuming K_ij=0
        % 1-> C1 2-> C3
        a11 = (x_C1(i)*x_C1(i)) * a_T_C1;
        a12 = (x_C1(i)*x_C3(i)) * ((a_T_C1*a_T_C3)^0.5);
        a21 = (x_C3(i)*x_C1(i)) * ((a_T_C3*a_T_C1)^0.5);
        a22 = (x_C3(i)*x_C3(i)) * a_T_C3;
        am = a11 + a12 + a21 + a22;
        
        %LHS of condition
        term1 = (v_final-bm) ^ 2;
        LHS = (R*T_guess)/term1;
        %RHS of condition
        term1 = 2*am*(v_final+bm);
        term2 = (v_final*v_final) + (2*v_final*bm) - (bm*bm);
        RHS = term1/(term2*term2);
        
        err = (LHS-RHS)/LHS;
        err = abs(err);
        
        if err<eps
            break;
        end
        
        T_guess = T_guess - 0.5;
        
     end
     
     Tc(i) = T_guess;
     
     %calculating pressure
     term1 = (R*Tc(i))/(v_final-bm);
     term2 = (v_final*v_final) + (2*v_final*bm) - (bm*bm);
     Pc(i) = term1 - (am/term2);
    
end

 

[breizh]

http://pillars.che.p...ns_of_State.pdf

 

Hi.

Consider the document above in addition to my previous email.

 

Breizh