|
|
| (25 intermediate revisions by 2 users not shown) |
| Line 1: |
Line 1: |
| {{task}} | | {{task}} |
| Work on this is still going on as a part of GSoC 2010, Project [[User:Ankitml|Ankit]] | | 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]]). |
| | See also [[FPROPS]]. |
|
| |
|
| Comments and suggestions are welcome | | Comments and suggestions are welcome |
| ==Algorithm== | | ==Overview== |
| ===Algorithm From Literature===
| | 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>: |
| *Guess a Pressure.
| | :<math> |
| *Find roots of volume using PR-EOS, if roots<3 guess P again.
| | p =\frac{{\bar R} T}{{\bar v}-b}-\frac{a(T)}{{\bar v}({\bar v}+b)+b({\bar v}-b)} |
| *Min V is the molar liquid volume, and the max V is vapour molar volume.
| | </math> |
| *get Zl and Zv.
| | where |
| *Get fugacities of both phases.
| | :<math>\begin{align} |
| *Pnew = Pold*fl/fv
| |
|
| |
|
| ===The current Algorithm=== | | a(T) &= 0.45724 \frac{{\bar R}^2{T_c}^2}{p_c} \alpha \left(T \right) \\ |
| *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=== | | \alpha &= \left( 1+\kappa \left( 1-\sqrt{\frac{T}{T_c}} \right) \right)^2 \\ |
| *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.
| | \kappa &= 0.37464+1.54226\omega - 0.26992\omega^2 \\ |
| *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== | | b &= \frac{0.0778\bar R T_c}{p_c} |
| ===Errors===
| | \end{align} |
| *Not working at some temperatures. [It should work on entire range from triple point to Critical Temperature]
| | </math> |
|
| |
|
| [Solution] - Currently trying got find out the reason for this behaviour. - [http://ascendwiki.cheme.cmu.edu/User:Ankitml Ankit ]
| | 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> |
|
| |
|
| *Not predicting Saturation Volumes for Liquids correctly.
| | ==Departure Functions== |
| [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.
| | 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] \\ |
|
| |
|
| ===Changes to be done to integrate with ASCEND=== | | 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] |
| *Data (Pc,Tc,Omega) to be taken from data files and only component name should be required in the code.
| | \end{align} |
| *No main function
| | </math> |
| *ASCEND wrapper
| | Clearly to evaluate these functions we need to be able to evaluate <math>\frac{da}{dT}</math> (checked, agrees with Sandler): |
| *Code should be in C, remove C++ functions.
| |
|
| |
|
| =Code=
| | :<math> |
| <source lang="c">
| | \frac{da}{dT}= -0.45724 \frac{{\bar R}^{2} {T_c}^{\frac{3}{2}} }{p_c} \kappa \frac{\sqrt{\alpha} }{ \sqrt{T}} |
| #include <stdio.h>
| | </math> |
| #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 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