|
|
| Line 54: |
Line 54: |
| \end{align} | | \end{align} |
| </math> | | </math> |
| So, for given values of <math>R</math>, <math>T_c</math>, <math>T_r</math>, <math>P_c</math> and <math>\omega</math> we can solve for <math>Z</math>. What next??
| | For any given pair of P,V and T we can solve for the other. Other thermodynamic properties can be obtained from the departure functions. |
| ==Algorithm==
| |
| ===Algorithm From Literature===
| |
| *Guess a Pressure.
| |
| *Find roots of volume using PR-EOS, if roots<3 guess P again.
| |
| *Min V is the molar liquid volume, and the max V is vapour molar volume.
| |
| *get Zl and Zv.
| |
| *Get fugacities of both phases.
| |
| *Pnew = Pold*fl/fv
| |
| | |
| ===The current Algorithm===
| |
| *Calculate PV data using PR_EOS from 0.001 to 6666 cm3/mol.
| |
| *Get the range of Pressure in which volume has 3 roots.
| |
| *Give this range as lower and upper bounds to a golden search algorithm, minimizing the difference in fugacities.
| |
| *To calculate fugacities a non iterative analytical expression was derived for PR-EOS.
| |
| | |
| ===Issues That Forced to Deviate From Literature===
| |
| *Find a non iterative way to solve a cubic Equation for volume. The normal iterative one makes the algorithm very slow, as it is required to solve the volume cubic many number of times to arrive at the solution. So right now generating PV data for a range with a constant increment in V.
| |
| *This also helps in guessing the first P in the range of 3 roots, which is required as the first step.
| |
| *The golden search algorithm is just experimental over the standard Pnew = Pold*Fl/Fv iteration, which requires a good guess whereas Golden search requires a bracket in which solution lies. It is yet to be tested which one is better.
| |
| | |
| ==Limitations==
| |
| ===Errors===
| |
| *Not working at some temperatures. [It should work on entire range from triple point to Critical Temperature]
| |
| | |
| [Solution] - Currently trying got find out the reason for this behaviour. - [http://ascendwiki.cheme.cmu.edu/User:Ankitml Ankit ]
| |
| | |
| *Not predicting Saturation Volumes for Liquids correctly.
| |
| [Solution] - Its a problem with Peng Robinson EOS, a modification to PR_EOS was proposed by G. Heyen in the paper, "Liquid and Vapor Properties from a cubic Equation of state", in 2nd International Conference on Phase Equilibria and Fluid Properties, 1980.
| |
| | |
| ===Changes to be done to integrate with ASCEND===
| |
| *Data (Pc,Tc,Omega) to be taken from data files and only component name should be required in the code.
| |
| *No main function
| |
| *ASCEND wrapper
| |
| *Code should be in C, remove C++ functions.
| |
| | |
| =Code=
| |
| <source lang="c">
| |
| #include <stdio.h>
| |
| #include <math.h>
| |
| #include <iostream>
| |
| using namespace std;
| |
| | |
| //Global Variables
| |
| double Tc,Pc,Omega,R,Tre,b,a,RT,T,P[20000],V[20000],Phil,Phiv,Vl,Vv;
| |
| int j;
| |
| | |
| int noRoots(double Psat); //gives out number of roots (of V) at a particular Pressure
| |
| 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
| |
| 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()
| |
| {
| |
| | |
| Tc = 304.1282; // Kelvin
| |
| Pc = 73.773; // bar
| |
| Omega = 0.22394; // acentric factor
| |
| R = 83.14; //bar[cm3]/k/mol
| |
| T = 260;
| |
| RT = R*T;
| |
| Tre = T/Tc;
| |
|
| |
| double p[20000],v[20000];
| |
| //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<20000;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;
| |
| pEnd = pStart;
| |
| deltaP = 1.0;
| |
| for(int n = 0;n<2;n++)
| |
| {
| |
| for(double pp = pEnd;pp<60.0;pp=pp+deltaP)
| |
| {
| |
| int nor;
| |
| nor = noRoots(pp);
| |
| if (nor<3)
| |
| {
| |
| pEnd = pp - deltaP;
| |
| pp = 5000;
| |
| }
| |
| }
| |
| deltaP = deltaP/10;
| |
| }
| |
| | |
| | |
| | |
| 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);
| |
| 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';
| |
| cout<<T<<'\t'<<sol<<'\t'<<Vl<<'\t'<<Vv<<'\t'<<fabs(Phil-Phiv)<<'\n';
| |
| | |
| return 0;
| |
| }
| |
| | |
| int noRoots(double Psat)
| |
| {
| |
| double Vroot[3];
| |
| Vroot[0] = -99; Vroot[1] = -99; Vroot[2] = -99;
| |
| 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)
| |
| {
| |
| double P,m,tmp;
| |
| 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)
| |
| { 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;
| |
| 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==
| |
| | |
| * Maybe this is of interest: http://pubs.acs.org/doi/abs/10.1021/ed066p54.2 ? [[User:Jpye|Jpye]] 06:39, 5 August 2010 (UTC)
| |
|
| |
|
| ==Comparisons== | | ==Comparisons== |
|
| |
|
| [[Category:Development]] | | [[Category:Development]] |
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 one of several equations of state that can be expressed in the form
- <math>
\begin{align}
P=\frac{RT}{V-b}-\frac{a}{V^2+ubV+wb^2}
\end{align}
</math>
or
- <math>
\begin{align}
Z^3-(1+B^*-uB^*)Z^2+(A^*+wB^{*2}-uB^*-uB^{*2})Z-A^*B^*-wB^{*2}-wB^{*3}=0
\end{align}
</math>
Where
- <math>
\begin{align}
A^*=\frac{aP}{R^2T^2}
\end{align}
</math>
- <math>
\begin{align}
B^*=\frac{bP}{RT}
\end{align}
</math>
In the Peng-Robinson values for <math>u,</math> <math>w</math>, <math>b</math> and <math>a</math> are set as
- <math>
\begin{align}
u=2
\end{align}
</math>
- <math>
\begin{align}
w=-1
\end{align}
</math>
- <math>
\begin{align}
b=\frac{0.7780RT_c}{P_c}
\end{align}
</math>
- <math>
\begin{align}
a=\frac{0.45724R^2T_c^2}{P_c}[1+f\omega(1-T_r^{\frac{1}{2}})]^2
\end{align}
</math>
- <math>
\begin{align}
f\omega=0.37464+1.54226\omega-0.26992\omega^2
\end{align}
</math>
For any given pair of P,V and T we can solve for the other. Other thermodynamic properties can be obtained from the departure functions.
Comparisons