|
|
| (36 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]]). |
| ===The general Algorithm===
| | See also [[FPROPS]]. |
| *Guess a Pressure.
| |
| *Find roots of volume using PR-EOS, if roots<3 guess P again.
| |
| *Min V is the liquid density at that P, and the max V is vapour.
| |
| *get Zl and Zv.
| |
| *Get fugacities of both phases.
| |
| *Pnew = Pold*fl/fv
| |
|
| |
|
| ===The current Algorithm=== | | Comments and suggestions are welcome |
| *Calculate PV data using PR_EOS from 0.001 to 6666 cm3/mol.
| | ==Overview== |
| *Get the range of Pressure in which volume has 3 roots.
| | 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>: |
| *Give this range as lower and upper bounds to a golden search algorithm, minimizing the difference in fugacities.
| | :<math> |
| *To calculate fugacities a non iterative analytical expression was derived for PR-EOS.
| | 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} |
|
| |
|
| ==Limitations== | | a(T) &= 0.45724 \frac{{\bar R}^2{T_c}^2}{p_c} \alpha \left(T \right) \\ |
| ===Errors===
| |
| ===Possible Changes that can be done to remove limitations/errors===
| |
| ===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= | | \alpha &= \left( 1+\kappa \left( 1-\sqrt{\frac{T}{T_c}} \right) \right)^2 \\ |
| <source lang="c">
| | |
| #include <stdio.h>
| | \kappa &= 0.37464+1.54226\omega - 0.26992\omega^2 \\ |
| #include <math.h>
| |
| #include <iostream>
| |
| using namespace std;
| |
|
| |
|
| | 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> |
|
| |
|
| int noRoots(double Psat); //gives out number of roots (of V) at a particular Pressure
| | ==Departure Functions== |
| double Pcalc(double T, double V); // Calculates Pressure at Given T and V using PR-EOS
| | 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. |
| int PhiCalc(double PsatGuess); //Calculates fugacities at a given Pressure and Stores it in global variables
| | The departure functions of the Peng-Robinson equation of state are as follows: |
| double phiDiff(double Pguess); //Finds fugacity difference at a particular P, uses function PhiCalc and is required for golden search algorithm
| | :<math> |
| 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.
| | \begin{align} |
| //Global Variables
| | 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] \\ |
| double Tc,Pc,Omega,R,Tre,b,a,RT,T,P[20000],V[20000],Phil,Phiv,Vl,Vv;
| |
| int j;
| |
|
| |
|
| int main()
| | 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} |
| Tc = 408.2; // Kelvin
| | </math> |
| Pc = 36.5; // bar
| | Clearly to evaluate these functions we need to be able to evaluate <math>\frac{da}{dT}</math> (checked, agrees with Sandler): |
| Omega = 0.183; // acentric factor
| |
| R = 83.14; //bar[cm3]/k/mol
| |
| T = 300.1;
| |
| Tre = T/Tc;
| |
| RT = R*T;
| |
|
| |
| 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;
| | :<math> |
| pEnd = pStart;
| | \frac{da}{dT}= -0.45724 \frac{{\bar R}^{2} {T_c}^{\frac{3}{2}} }{p_c} \kappa \frac{\sqrt{\alpha} }{ \sqrt{T}} |
| deltaP = 1.0;
| | </math> |
| for(int n = 0;n<2;n++)
| |
| {
| |
| for(double pp = pEnd;pp<60.0;pp=pp+deltaP)
| |
| {
| |
| int nor;
| |
| nor = noRoots(pp);
| |
| if (nor==1)
| |
| {
| |
| pEnd = pp - deltaP;
| |
| 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;
| | ==Comparisons== |
| }
| |
| | |
| 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)
| | [[Category:Development]] |
| {
| |
| 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==
| |
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