## Featured Articles

Check out the latest featured articles.

## New Article

Product Viscosity vs. Shear

## Featured File

Vertical Tank Selection

## New Blog Entry

Gas Turbine Re-Ratings- posted in Ankur's blog

4

# Calculating Mixture Critical Temperature And Pressure Using Equation O

critical temperature equation of state critical pressure peng-robinson critical state

4 replies to this topic
|

### #1 suman93

suman93

Brand New Member

• Members
• 2 posts

Posted 07 August 2018 - 07:58 PM

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

### #2 Art Montemayor

Art Montemayor

Gold Member

• 5,596 posts

Posted 08 August 2018 - 01:47 PM

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.

### #3 breizh

breizh

Gold Member

• ChE Plus Subscriber
• 3,959 posts

Posted 09 August 2018 - 07:02 AM

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.

### #4 suman93

suman93

Brand New Member

• Members
• 2 posts

Posted 09 August 2018 - 08:13 PM

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 d2P/dv2 = 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

P = { RT/(v- }  -   {  a / (v2 + 2vb - b2) }

2. I find expressions for dP/dv = 0 and d2P / dv2= 0.

3. I equate them to end up with a cubic equation:  v3 - (3b)v2 - (3b2)v - (3b3) = 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

```

### #5 breizh

breizh

Gold Member

• ChE Plus Subscriber
• 3,959 posts

Posted 12 August 2018 - 02:26 AM

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

Hi.

Consider the document above in addition to my previous email.

Breizh