PengRobinson EOS in FPROPS: Difference between revisions

From ASCEND
Jump to navigation Jump to search
← Older edit
Ankitml (talk | contribs)
33 edits
No edit summary
 
(41 intermediate revisions by 3 users not shown)
Line 1: Line 1:
=Documentation=
{{task}}
==Algorithm==
Work on this is went on as a part of GSoC 2010 (Project [[User:Ankitml|Ankit]]) and continues as part of GSoC 2011 ([[User:richardTowers|Richard Towers]]).
==Limitations==
See also [[FPROPS]].
*limitations
*Possible Changes that can be done to remove
*Errors


=Code=
Comments and suggestions are welcome
<source lang="c">
==Overview==
#include <stdio.h>
The Peng-Robinson EOS is a cubic equation of state in that it contains volume terms to the third power. It is usually expressed to give pressure in terms of temperature and molar volume <math>{\bar v}</math>:
#include <math.h>
:<math>
#include <iostream>
p =\frac{{\bar R} T}{{\bar v}-b}-\frac{a(T)}{{\bar v}({\bar v}+b)+b({\bar v}-b)}
using namespace std;
</math>
where
:<math>\begin{align}


//Global Variables
a(T) &= 0.45724  \frac{{\bar R}^2{T_c}^2}{p_c} \alpha \left(T \right) \\
double Tc,Pc,Omega,R,Tre,b,a,RT,T,P[25000],V[25000],Phil,Phiv,Vl,Vv;
int j;


int noRoots(double Psat); //gives out number of roots (of V) at a particular Pressure
\alpha &= \left( 1+\kappa \left( 1-\sqrt{\frac{T}{T_c}} \right) \right)^2 \\
double Pcalc(double T, double V);  // Calculates Pressure at Given T and V using PR-EOS
int PhiCalc(double PsatGuess); //Calculates fugacities at a given Pressure and Stores it in global variables
\kappa &= 0.37464+1.54226\omega - 0.26992\omega^2 \\
double phiDiff(double Pguess); //Finds fugacity difference at a particular P, uses function PhiCalc and is required for golden search algorithm
double goldSer(double a,double c,double b); //Golden search algorithm, function to be minimized is phiDiff and a is the lower bound, b is upper bound, c is the center.


int main()
b &= \frac{0.0778\bar R T_c}{p_c}
{
\end{align}
        Tc = 408.2;  // Kelvin
</math>
        Pc = 36.5;  // bar
        Omega = 0.183; // acentric factor
        R = 83.14; //bar[cm3]/k/mol
T = 300.1;
        Tre = T/Tc;
RT = R*T;
double p[25000],v[25000];
//following for loop calculates P V data till v=6500cm3/mol and stores it in P,V arrays only if P is within bounds of 60 and -160
j = 0;
        for (int i=0;i<25000;i++)
        {
                v[i] = (i + 0.001)/3;
                p[i] = Pcalc(T,v[i]);
                if(p[i]<60.0 && p[i]>-160.0)
                {
                        P[j] = p[i];
                        V[j] = v[i];
                        j++;
                }
        }
//finding the bounds for P
        double pStart,deltaP;
        pStart = 0.01;
        deltaP = 1.0;
        for(int n = 0;n<3;n++)
        {
                for(double pp = pStart;pp<60.0;pp=pp+deltaP)
                {
                        int nor;
                        nor = noRoots(pp);
                        if (nor==3)
                        {
                                pStart = pp - deltaP;
                                pp = 59.0;
                        }
                }
                deltaP = deltaP/10;
        }
        pStart = pStart + 0.01;


        double pEnd;
It is sometimes more convenient to express the equation as a cubic polynomial in terms of compressibility factor <math>Z</math>
        pEnd = pStart;
:<math>
        deltaP = 1.0;
Z^3+(-1+B)Z^2+(A-3B^2-2B)Z-(AB-B^2-B^3)=0
        for(int n = 0;n<2;n++)
</math>
        {
in which
                for(double pp = pEnd;pp<60.0;pp=pp+deltaP)
:<math>
                {
\begin{align}
                        int nor;
A &= \frac{a \left(T \right) p}{({\bar R} T)^2} \\
                        nor = noRoots(pp);
B &= \frac{b p}{{\bar R} T} \\
                        if (nor==1)
Z &= \frac{p {\bar v}}{{\bar R} T}
                        {
\end{align}
                                pEnd = pp - deltaP;
</math>
                                pp = 59.0;
                        }
                }
                deltaP = deltaP/10;
        }
cout<<"At Temp = "<<T<<" K"<<'\n';
int sandbox;
double Pa,Pb,Pc,sol,phi,resphi;
phi = (1 + sqrt(5))/2;
        resphi = 2 - phi;
Pa = pStart;
Pb = pEnd;
Pc = (Pb - Pa)*resphi + Pa;
sol = goldSer(Pa,Pc,Pb);
cout<<"Psat = "<<sol<<" bar"<<'\n';
sandbox=PhiCalc(sol);
cout<<"Molar Volume (L) = "<<Vl<<" cm3/mol"<<'\n';
cout<<"Molar Volume (V) = "<<Vv<<" cm3/mol"<<'\n';
cout<<"Difference in fugacities of vapour and liquid = "<<fabs(Phil-Phiv)<<'\n';


return 0;
==Departure Functions==
}
Departure functions represent the departure of the ''real'' properties from the ''ideal'' properties - i.e the properties of a fluid at zero pressure or infinite molar volume.
The departure functions of the Peng-Robinson equation of state are as follows:
:<math>
\begin{align}
H_{m}-H_{m}^{\text{ideal}}&={\bar R} T(Z-1)+\frac{T\left(\frac{da}{dT}\right)-a}{2\sqrt{2}b}\ln\left[\frac{Z+(1+\sqrt{2})B}{Z+(1-\sqrt{2})B}\right] \\


int noRoots(double Psat)
S_{m}-S_{m}^{\text{ideal}}&={\bar R} \ln (Z-B)+\frac{\frac{da}{dT}}{2\sqrt{2}b}\ln\left[\frac{Z+(1+\sqrt{2})B}{Z+(1-\sqrt{2})B}\right]
{
\end{align}
        double Vroot[3];
</math>
        Vroot[0] = -99; Vroot[1] = -99; Vroot[2] = -99;
Clearly to evaluate these functions we need to be able to evaluate <math>\frac{da}{dT}</math> (checked, agrees with Sandler):
        int rootCount = 0;
        for (int i = 0; i < j-1;i++)
        {
                if((Psat<=P[i+1] && Psat>=P[i]) || (Psat>=P[i+1] && Psat<=P[i]))
                {
                        Vroot[rootCount] = V[i]+(Psat-P[i])*(V[i+1]-V[i])/(P[i+1]-P[i]);
                        rootCount++;
                }
        }
        return rootCount;
}


double Pcalc(double T, double V)
:<math>
{
\frac{da}{dT}= -0.45724 \frac{{\bar R}^{2} {T_c}^{\frac{3}{2}} }{p_c} \kappa \frac{\sqrt{\alpha} }{ \sqrt{T}}
        double P,m,tmp;
</math>
        b = 0.0778*R*Tc/Pc;
        m = 0.37464+1.54226*Omega-0.26992*pow(Omega,2);
        tmp = 1+m*(1-sqrt(Tre));
        a = 0.45724*R*R*Tc*(Tc/Pc)*pow(tmp,2);
        P = R*T/(V-b) - a/(V*(V+b)+b*(V-b));
        return P;
}


int PhiCalc(double PsatGuess)
==Comparisons==
{      double Vroot[3],A,B;
        B = PsatGuess*b/RT;
        A = a*PsatGuess/pow(RT,2);
        Vroot[0] = -99; Vroot[1] = -99; Vroot[2] = -99;
        int rootCount = 0;
        for (int i = 0; i < j-1;i++)
        {
                if((PsatGuess<=P[i+1] && PsatGuess>=P[i]) || (PsatGuess>=P[i+1] && PsatGuess<=P[i]))
                {
                        Vroot[rootCount] =V[i]+(PsatGuess-P[i])*(V[i+1]-V[i])/(P[i+1]-P[i]);
                        rootCount++;
                }
        }
 
        double Zv,Zl;
        double term_a,term_b,term_c;
 
        Vl = Vroot[0];  //liquid volume is the smaller root
        Vv = Vroot[2];  // vapour volume is the largest one
 
        Zl = PsatGuess*Vl/(R*T);  //the corresponding Z
        Zv = PsatGuess*Vv/(R*T);


        term_a = (1 + sqrt(2))*B;
[[Category:Development]]
        term_b = (1 - sqrt(2))*B;
        term_c = 2*sqrt(2)*B;
 
        Phil = (Zl-1) - log(Zl-B) - A/term_c*log((Zl+term_a)/(Zl+term_b));
        Phiv = (Zv-1) - log(Zv-B) - A/term_c*log((Zv+term_a)/(Zv+term_b));
return 0;
}
 
double phiDiff(double Pguess)
{
double phiD;
int sandbox;
sandbox = PhiCalc(Pguess);
phiD = fabs(Phil - Phiv);
return phiD;
}
 
double goldSer(double a,double c,double b)
{
        double d,phi,resphi;
phi = (1 + sqrt(5))/2;
        resphi = 2 - phi;
        if (fabs(a-b)<0.000001)
                return ((a+b)/2);
 
        d = c + resphi*(b - c);
 
        if (phiDiff(d)<phiDiff(c))
                return goldSer(c,d,b);
        else
                return goldSer(d,c,a);
}
</source>
 
==Possible Improvements in the code==
==Comparisons==

Latest revision as of 23:38, 13 January 2013

This article is about planned development or proposed functionality. Comments welcome.

Work on this is went on as a part of GSoC 2010 (Project Ankit) and continues as part of GSoC 2011 (Richard Towers). See also FPROPS.

Comments and suggestions are welcome

Overview

The Peng-Robinson EOS is a cubic equation of state in that it contains volume terms to the third power. It is usually expressed to give pressure in terms of temperature and molar volume <math>{\bar v}</math>:

<math>

p =\frac{{\bar R} T}{{\bar v}-b}-\frac{a(T)}{{\bar v}({\bar v}+b)+b({\bar v}-b)} </math> where

<math>\begin{align}

a(T) &= 0.45724 \frac{{\bar R}^2{T_c}^2}{p_c} \alpha \left(T \right) \\

\alpha &= \left( 1+\kappa \left( 1-\sqrt{\frac{T}{T_c}} \right) \right)^2 \\

\kappa &= 0.37464+1.54226\omega - 0.26992\omega^2 \\

b &= \frac{0.0778\bar R T_c}{p_c} \end{align} </math>

It is sometimes more convenient to express the equation as a cubic polynomial in terms of compressibility factor <math>Z</math>

<math>

Z^3+(-1+B)Z^2+(A-3B^2-2B)Z-(AB-B^2-B^3)=0 </math> in which

<math>

\begin{align} A &= \frac{a \left(T \right) p}{({\bar R} T)^2} \\ B &= \frac{b p}{{\bar R} T} \\ Z &= \frac{p {\bar v}}{{\bar R} T} \end{align} </math>

Departure Functions

Departure functions represent the departure of the real properties from the ideal properties - i.e the properties of a fluid at zero pressure or infinite molar volume. The departure functions of the Peng-Robinson equation of state are as follows:

<math>

\begin{align} H_{m}-H_{m}^{\text{ideal}}&={\bar R} T(Z-1)+\frac{T\left(\frac{da}{dT}\right)-a}{2\sqrt{2}b}\ln\left[\frac{Z+(1+\sqrt{2})B}{Z+(1-\sqrt{2})B}\right] \\

S_{m}-S_{m}^{\text{ideal}}&={\bar R} \ln (Z-B)+\frac{\frac{da}{dT}}{2\sqrt{2}b}\ln\left[\frac{Z+(1+\sqrt{2})B}{Z+(1-\sqrt{2})B}\right] \end{align} </math> Clearly to evaluate these functions we need to be able to evaluate <math>\frac{da}{dT}</math> (checked, agrees with Sandler):

<math>

\frac{da}{dT}= -0.45724 \frac{{\bar R}^{2} {T_c}^{\frac{3}{2}} }{p_c} \kappa \frac{\sqrt{\alpha} }{ \sqrt{T}} </math>

Comparisons

Retrieved from "https://ascend4.org/index.php?title=PengRobinson_EOS_in_FPROPS&oldid=4102"