FPROPS

2010-05-21T16:20:13Z

Ankitml:

<div style="float:right;border:solid 1pt green;background-color=#EEFFEE;width:18em;padding:12px">
'''Fluids supported by FPROPS'''

* {{src|models/johnpye/fprops/ammonia.c}}
* {{src|models/johnpye/fprops/hydrogen.c}}
* {{src|models/johnpye/fprops/nitrogen.c}}
* {{src|models/johnpye/fprops/water.c}}
* {{src|models/johnpye/fprops/carbondioxide.c}}

* Some more, yet to be merged: [http://ascendcode.cheme.cmu.edu/viewvc.cgi/code/branches/hongke/models/johnpye/fprops/ hongke:models/johnpye/fprops/]
</div>
'''FPROPS''' is a free open-source C-based library for evaluating thermodynamic properties of a number of pure substances via published data for the Helmholtz fundamental equation for those substances. It has been developed by [[User:Jpye|John Pye]] and will happily function as standalone code, but is also provided with [[external library]] code for ASCEND so that it can be used to access these accurate property correlations from within a [[MODEL]].

Currently FPROPS supports calculation of the properties of the pure substances shown at right. The properties that can be calculated are internal energy, entropy, pressure, enthalpy and Helmholtz energy, as well as various partial derivatives of these with respect to temperature and density.

FPROPS reproduces a limited subset of the functionality of commercial programs such as [http://www.nist.gov/srd/nist23.htm REFPROP], EES, FLUIDCAL, [[freesteam]], SteamTab, and others, but is open source software and licensed under the [http://www.gnu.org/licenses/gpl.html GPL]. Contributions are most welcome!

'''IMPORTANT NOTE:''' FPROPS does not yet implement calculation of the saturation curve, so it gives incorrect property evaluations anywhere within the saturation region. This is because the code for evaluating the 'Maxwell criterion' is not yet completely implemented.

FPROPS is included with ASCEND as of version 0.9.5.116, or you can browse the code at {{srcdir|models/johnpye/fprops/}} or access it [[VersionManagement|via subversion]] (there have been changes since the 0.9.5.116 release).

== Implementation ==

The calculation routines implement calculation of reduced Helmholtz energy φ by dividing it into ideal (<math>\phi^0</math>) and residual (<math>\phi^r</math>) parts, as described in IAPWS-95<a href="#_note-iapws95" title="">[1]</a>.

=== Ideal part ===

Calculation of the '''ideal part''' of the reduced Helmholtz energy is derived from the expression for the zero-pressure limit of reduced specific heat capacity at constant pressure, <math>c_p^0 / R</math>, where <math>R</math> is the specific gas constant. This quantity <math>{c_p^0}(T)/R</math> can be expressed as a series composed of two different types of terms:

* power terms, <math>c_i T^{t_i}</math>
* exponential terms, <math>\frac{b_i x_i^2 exp(-x_i)}{[1-exp(-x_i)]^2}</math> where <math>x_i = \frac{\beta_i}{T}</math>

An expression for calculating the ideal component of Helmholtz energy <math>a^0 \left(T , \rho \right)</math> from <math>c_p^0</math> is given by Span et al<a href="#_note-0" title="">[1]</a>:

:<math>\frac{a^0}{RT} = \phi^0 = \frac{h_0^0 \tau}{R T_c} - \frac{s_0^0}{R} - 1 + \ln {\frac{\delta / \delta_0}{\tau / \tau_0}} - \frac{\tau}{R} \int_{\tau_0}^{\tau}{\frac{c_p^0}{\tau^2}d\tau} + \frac{1}{R} \int_{\tau_0}^{\tau}{\frac{c_p^0}{\tau} d \tau}</math>

where the '''reduced temperature''' is <math>\tau = \frac{T^{*}}{T}</math> and the '''reduced density''' is <math>\delta = \frac{\rho}{\rho^{*}}</math>.

On the other hand, using the relations for the Helmholtz fundamental equation, one can see that

:<math>\frac{c_p^0}{R} = 1 -\tau^2 \phi^0_{\tau \tau}</math>

so coefficients of <math>c_p^0</math> can be calculated from coefficients of <math>\phi^0</math> (providing, we assume, that the fundamental equation has smooth first and second (and possibly third) derivatives, and that when terms of <math>{c_p^0}(T)</math> are integrated they don't 'interact'). One can go the other way, and determine <math>\phi^0</math> from <math>c_p^0</math>, with the exception of the integration constants, which would in that case be determined using specified values of enthalpy <math>h_ref</math> and entropy <math>s_ref</math> at a reference state <math>\left(T_{ref},\rho_{ref} \right)</math>

=== Residual part ===

For the '''residual part''' of the reduced Helmholtz energy, FPROPS currently supports three different kinds of terms:

* power terms, <math>a_i \tau^{t_i}\delta^{d_i}</math> (these are HelmoltzPowTerm with ''l''''i'' = 0)
* exponential terms, <math>a_i \tau^{t_i} \delta^{d_i} exp(-\delta^{l_i})</math> (these are HelmholtzPowTerm with <math>l_i \ne 0</math>)
* gaussian terms, <math>n_i \tau^{t_i} \delta^{d_i} exp[-\alpha_i(\delta - \epsilon_i)^2 - \beta_i(\tau - \gamma_i)^2]</math>
* critical terms, <math>n_i \Delta^{b_i} \delta \psi</math>, where <math>\Delta = \theta^2 + B_i \left[ \left( \delta - 1 \right)^2 \right]^{a_i}</math>, and <math>\theta = \left( 1- \tau \right) + A_i \left[ \left( \delta - 1 \right)^2 \right]^{\frac{1}{2 \beta_i}}</math> and <math>\psi = exp \left[-C_i \left(\delta-1 \right)^2 - D_i \left( \tau - 1 \right)^2 \right]</math>

These different terms are used in various combinations across a number of high-accuracy publications of thermodynamic property data. Overall, the residual series permits calculation of 'real' reduced Helmholtz energy <math>\phi \left(\tau,\delta \right) = \phi^0 \left(\tau,\delta \right) + \phi^r \left( \tau,\delta \right)</math>.

From the expression for <math>\phi</math>, we can then derive expressions for <math>p</math>, <math>h</math>, <math>u</math>, <math>s</math> and their partial derivatives with respect to density and temperature. These expressions are summarised in IAPWS-95 as well as by Tillner Roth et al<a href="#_note-0" title="">[1]</a>.

The source code for FPROPS is available via the ASCEND [[VersionManagement|subversion repository]], in the {{srcdir|models/johnpye/fprops}} directory. There is some further information on the implementation on the page about [[writing ASCEND external relations in C]].

=== Maxwell criterion ===

The Helmholtz function above applies for all regions outside the saturation region. The saturation region is determined by a set of equations call the Maxwell phase equilibrium criteria or Maxwell equal-area rule <a href="#_note-iapws95" title="">[1]</a><a href="#_note-0" title="">[1]</a>

Calculation of the phase equilibrium conditions requires iterative solution of three equations; a suitable algorithm for efficiently performing this has not yet been incorporated into FPROPS; contributions are welcome.

The Maxwell criterion is stated in IAPWS-95 as being the preferred approach for determining the location of the saturation line, but it could be possible that we can use equality of Gibbs Free Energy to work this out. Need to investigate this and work out a suitable approach.

== Testing the code ==

FPROPS includes a suite of self tests for each fluid that has been implemented. Some of the test data comes from the publications from which the correlations have been taken; other test data comes from using the commercial package [http://www.nist.gov/srd/nist23.htm REFPROP] to calculate some sample data points, and comparing the FPROPS output with that data.

To run the FPROPS tests for a given fluid, <tt>cd</tt> to the FPROPS directory (eg <tt>cd ~/ascend/models/johnpye/fprops</tdd>) then run the test script <tt>./test.py ''fluidname''</tt>, for example:

<source lang="a4c">./test.py hydrogen
./test.py nitrogen</source>
''Requires GNU scientific library(libgsl0-dev) to be installed in the system''.

The output will tell you if all calculated values have conformed to expected error margins.

Test code for each fluid is embedded at the end of the corresponding fluid file, for example {{src|models/johnpye/fprops/hydrogen.c}}.

== Adding new fluids ==

Currently, new fluids must be added by declaring constant data structures in C code.

An example of a fluid declaration is shown in {{src|models/johnpye/fprops/hydrogen.c}}.

The data structures which are to be filled are declared in {{src|models/johnpye/fprops/ideal.h}} and {{src|models/johnpye/fprops/helmholtz.h}}.

'''Note''': when adding new fluids it is important to observe that different papers give different 'zero points' for enthalpy and entropy for the substances in question, and the zero points vary from paper to paper. These Helmholtz equations of state can be adjusted for different 'zero points' by adjusting the constant and linear parameters ('<tt>c</tt>' and '<tt>m</tt>') in the IdealData object that is provided in your data structures. The equation to fix these values is

:<math>\frac{\Delta h}{RT} = \tau \Delta \phi_{\tau}^0 = T^{*} \Delta a_2^0</math>

:<math>\frac{\Delta s}{R} = \tau \Delta \phi_{\tau}^0 - \Delta \phi^0 = - \Delta a_1^0</math>

where we assume

:<math>\phi^0(\tau,\delta) = a_1^0 + a_2^0 \tau + \mathrm{other\ terms\ in\ \tau, \delta}</math>

In other words, to correct an offset in ''h'' and ''s'', these formulae tell you the offset to apply to the constant and linear terms in φ0(τ,δ).

== Using within ASCEND ==

To use FPROPS from ASCEND, see the example code in {{src|models/johnpye/fprops/fprops_test.a4c}}. Essentially you must declare a [[DATA]] instance into which you must place the name of the substance for which you wish to calculate properties. Then you can use the functions <tt>helmholtz_p</tt>, <tt>helmholtz_u</tt>, <tt>helmholtz_h</tt>, <tt>helmholtz_p</tt>, <tt>helmholtz_s</tt>, <tt>helmholtz_a</tt> to calculate the property you need.

'''Note'''. When solving properties in 'reverse' mode (eg solving for ρ,''T'' when ''p'',''h'' are know), you will get significantly better results by switching the QRSlv solver parameter for 'convergence test' (convopt) to 'RELNOM_SCALE'. This significantly improves ASCEND's ability to solve these problems.

== Using from C ==

Using the FPROPS code from C should also be very straightforward; example code is shown in {{src|models/johnpye/fprops/test.c}}.

== Further work ==
{{task}}

FPROPS is still under development! We aim to implement a few more features including

* add support for '''calculation of the saturation curve''' using Maxwell phase-equilibrium condition (Assign to Ankit)
* <s>support for Water using IAPWS-95 correlation (partially complete, but needs support for another type of residual function term).</s>
* <s>refactor the 'HelmholtzExpTerm' which is actually a double of 'HelmholtzGausTerm' with epsilon = 1.</s>
* <s>add support for calculation of partial derivatives dp/dT, dp/drho, dh/dT, dh/drho, du/dT, du/drho,</s> ds/dT, ds/drho, da/dT, da/drho <s>to improve convergence of ASCEND models that calculate p, h, u,</s> s, a etc.
* investigate apparent disagreement between cp0 for '''ammonia''' in FPROPS and values from REFPROP (is it consistent with Tillner-Roth?)
* investigate apparent very small discrepancy between REFPROP and IAPWS95 in value of speed of sound near critical point.
* move specification of the '''zero point''' out of IdealData and into HelmholtzData.
* implement a '''pre-calculation routine''' where for example coefficients of <math>\phi^0</math> can be prepared?
* implement support for thermophysical properties '''conductivity, viscosity, surface tension,''' and possibly others.
* <s>implement support for calculation of the speed of sound.</s>
* resolve the confusion between T* and T_crit in the data structures. Probably all correlations will be using T_crit -- but is that sure?

Further down the track we would like to support:

* reading property correlations from a '''database''' (SQLite?) or text file. (Assigned to Shrikanth)
* some simpler equations of state such as MBWR or Peng-Robinson.
* more fluids of general interest.
* '''mixing models''' for multi-phase, multi-component systems (or else incorporate FPROPS into exist code that supports this).

We'd be very pleased to be able to collaborate on this work with others; the code has been made deliberately able to stand alone; it has no dependency on any other code from ASCEND, and is suitable for incorporation into other open source programs.

== References ==

== See Also ==

* [[Thermodynamics with ASCEND]]

[[Category:Proposed]]
[[Category:Development]]
[[Category:Documentation]]