User:Smuratet
Sean Muratet is a GSOC 2012 student from Huntsville, AL, USA. He is currently pursuing a Master of Science degree in Modeling and Simulation at Center for Modeling, Simulation, and Analysis of the University of Alabama in Huntsville. Mr. Muratet graduated from the University of Alabama in Huntsville with a BSE in Chemical Engineering in 2011. He has previously worked as a Graduate Research Assistant at the Center for Modeling, Simulation, and Analysis at UAH.
GSOC 2012 Goals
- Complete and test Peng Robinson EOS code currently in the FPROPS 2 branch.
- Integrate PR code into Ascend
- Code additional cubic EOS('s) into Ascend, and possible MWBR
- Continue the work of Richard_Towers on an XML based component file format
Progress Reports
6 July 2012
To this point, the main effort has been to determine the failure conditions of the current Peng Robinson code. A preliminary investigation dealing strictly with the case of water found that known accurate values of density at temperatures within the 2 phase region return negative pressure values. My first line involved tracing possible unit and conversion errors within the code. This turned out to be a dead end. At this point, I am fairly certain of the accuracy and appropriateness of the units within the PR code, at least as it applies to the saturation and pressure functions. Other function such as enthalpy and entropy have not been fully vetted at this point, but that will come soon. After more debugging, it became apparent that densities outside a certain range returned negative pressure values for a given temperature. This brings up an interesting numerical issue. The PR code takes a density and temperature as arguments to return pressure. This seems to be the tail wagging the dog in a sense. Obviously, at a certain temperature and density, the thermodynamic state of the system is set (Two intensive variables=system defined). However, in the two phase region, more than one density can map to the same pressure. In addition, one density root of the PR equation (the middle value of the three) is not physically significant.
This puts the onus on the user to ensure an appropriate density value is entered, consistent with the PR EOS. This may not be the case in many instances for someone not very familiar with the PR EOS in conjunction with a given fluid. It could be difficult to know a priori how to proceed. I will briefly try to provide examples of common situations which could prove troublesome:
1. Taking the fluid water at 298.15 K and density 997. This is generally known to be a situation where a two phase region is possible. With the PR EOS, this density would return a negative pressure. This could mean one of a few things. Either the fluid is in tension (a physical possibility in general but not here) or we are outside the range of physical significance of the PR EOS. Based on my own investigations, the maximum density must be less than 1/b or in this case approximately 951.752 where the PR EOS exhibits asymptotic behavior. There seems to be an additional asymptote around 2297, beyond which a positive pressure may be found.
2. Under the same temperature conditions, a density of 40 is a root of pressure 1.27 bar. However, this is a pressure root which is not physically significant. The two other roots 12.37 and 850.66 are the physically significant roots. However, without some mechanism to stop it, the user could enter this density and a pressure would be returned which is neither correct nor proper.
3. In general, a cubic will have three roots. For the case of the PR EOS a solution will have three real density roots or one real and two complex roots. If three real roots are found, it seems to indicate a two phase region although there is some disagreement as to whether the phases are in equilibrium as a consequence. If two complex solutions and one real is found, a single phase is indicated with density of the real root. However, a condition arises that is unexpected (to me at least). There is a case when three real roots are returned, but the system should be single phase. For example with carbon dioxide at 300 K and a density of 555, three real roots are returned. one is 555, the other two are -636 and -889.66. Just below at 554, two complex and one real roots (554) are returned. This point at which the transition is made seems to vary with temperature. The pressure from PR EOS is 237 bar (.237 GPa) at 555. This pressure does seem to be well within the range of possibility in terms of a pressure which could be reached. I don't currently know the implications of this behavior, or whether to ignore it or not. It is a situation though which would have to be dealt with in any code to try and find a saturation condition. We then can't simply say if three real roots are returned from a single pressure that the system is two phase and that the highest and lowest roots are the liquid and gas phases, respectively. We could be in the situation just described, and return values of -889.66 and 555 as vapor and liquid densities. Right now I believe it is sufficient to ignore roots which are negative or larger than 1/b. Update: These values were obtained in error on my part. This may not be applicable
Generally speaking, the PR EOS does not describe water very well at all. Still, the same problems might occur for other fluids which are better described. The end user must have some knowledge of equations of state generally to use them. We cannot hope to protect people from the poor results which will come from misuse, while at the same time providing software which is useful. I hope to write code which is a balance between the two principles. At this time, I think that means throwing an error for a density input at or above that which is 1/b. There may be a need to write a more elegant routine, but this should be the beginning.
15 August 2012 The Peng-Robinson EOS has been implemented within FPROPS and is now available from within ascend. The functions of specific volume and temperature implemented include pressure, enthalpy, entropy, internal energy, gibbs energy, helmholtz energy,constant pressure heat capacity, constant volume heat capacity, and fluid sonic velocity. Specifics can be found in the following sections. Many of these functions had been previously coded, but were not working in that a saturation condition could not be found via the implemented code through the PR EOS. I have solved this problem by writing a saturation routine that is specific to cubic eos's, but may well function in other cases. The routine requires an initial guess for vapor phase density provided by the chouaieb correlation. A pressure calculation is made from the PR EOS based on this density, and the liquid phase root is found. The fugacity coefficients of each phase are compared and if not equivalent i.e. <math>\left(\phi_{l}=\phi_{v}\right)</math> then values are successively substituted until equivalence is obtained. Near the critical point, and in some other problem cases, this method alone will not converge. In these cases, Brent's method is used to find the maximum and minimum of <math>\left(\frac{dP}{d\bar{v}}\right)_{T}</math> at the requested temperature. Within these bounds, we can be assured three roots exist. The average of the two pressures (maximum and minimum) is used to continue successive substitution. Using this method, the saturation routine will converge fairly close to the critical point <math>\left(.999984T_{c}\right)</math> in tested cases.
Overview
The following is some explanation of Peng-Robinson functions implemented within my branch (with some information taken from : The Peng-Robinson EOS is a cubic equation of state in that it contains volume terms to the third power. It is expressed to give pressure in terms of temperature and specific 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=\frac{p {\bar v}}{{\bar R} T}</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}{{\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>
Additional Functions
From the above functions, many other properties can be determined by their respective definitions
- <math>
\begin{align} Internal~Energy\,~U \\ U=H-PV \\
Helmholtz~Energy\,~A \\ A=U-TS \\
Gibbs~Energy\,~G \\ G=H-TS \\ \end{align} </math> We can also determine the heat capacities and fluid sonic velocity (speed of sound in a fluid). Currently, the PR EOS in FPROPS uses the rpp.c data file for components by default. This file contains constant pressure heat capacity correlation information. This will provide the ideal portion of the constant pressure heat capacity. The real portion can be determined from the definition:
- <math>
C_{p}=\left(\frac{\partial H}{\partial T}\right)_{P}=\frac{T}{2\sqrt{2}b}\left[\frac{d^{2}a}{dT^{2}}\right]\ln\left[\frac{Z+2.414B}{Z-0.414B}\right]+\left[\frac{R(M-N)^{2}}{M^{2}-2A(Z+B)}\right]-R+\overset{i}{C_{p}} </math> where:
- <math>
\begin{align} N=\frac{B}{Rb}\frac{da}{dT} \\ M=\frac{Z^{2}+2BZ-B^{2}}{Z-B} \\ \left[\frac{d^{2}a}{dT}\right]=\left[\frac{a(T_{c})\kappa}{2TT_{c}^{.5}}\right]\left[\frac{\alpha^{.5}}{T^{.5}}+\frac{\kappa}{T_{c}^{.5}}\right] \\ \end{align} </math>