|
|
| (35 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]]). |
| ===Algorithm From Literature===
| | 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} |
|
| |
|
| ===Issues That Forced to Deviate From Literature=== | | a(T) &= 0.45724 \frac{{\bar R}^2{T_c}^2}{p_c} \alpha \left(T \right) \\ |
| *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.
| |
|
| |
|
| ==Limitations== | | \alpha &= \left( 1+\kappa \left( 1-\sqrt{\frac{T}{T_c}} \right) \right)^2 \\ |
| ===Errors===
| | |
| ===Possible Changes that can be done to remove limitations/errors===
| | \kappa &= 0.37464+1.54226\omega - 0.26992\omega^2 \\ |
| ===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= | | b &= \frac{0.0778\bar R T_c}{p_c} |
| <source lang="c">
| | \end{align} |
| #include <stdio.h>
| | </math> |
| #include <math.h>
| |
| #include <iostream>
| |
| using namespace std;
| |
|
| |
|
| | 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] \\ |
|
| |
|
| int noRoots(double Psat); //gives out number of roots (of V) at a particular Pressure
| | 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] |
| double Pcalc(double T, double V); // Calculates Pressure at Given T and V using PR-EOS
| | \end{align} |
| int PhiCalc(double PsatGuess); //Calculates fugacities at a given Pressure and Stores it in global variables
| | </math> |
| double phiDiff(double Pguess); //Finds fugacity difference at a particular P, uses function PhiCalc and is required for golden search algorithm
| | Clearly to evaluate these functions we need to be able to evaluate <math>\frac{da}{dT}</math> (checked, agrees with Sandler): |
| 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.
| |
| //Global Variables | |
| double Tc,Pc,Omega,R,Tre,b,a,RT,T,P[20000],V[20000],Phil,Phiv,Vl,Vv;
| |
| int j;
| |
|
| |
|
| int main()
| | :<math> |
| { | | \frac{da}{dT}= -0.45724 \frac{{\bar R}^{2} {T_c}^{\frac{3}{2}} }{p_c} \kappa \frac{\sqrt{\alpha} }{ \sqrt{T}} |
| 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[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;
| | ==Comparisons== |
| 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==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;
| |
| }
| |
| | |
| 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
| | [[Category:Development]] |
| 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==
| |
| ==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