View source for FPROPS

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

<small>38 high-accuracy Helmholtz fluids:</small> {{fpropsfluid|acetone}},
{{fpropsfluid|ammonia}},
{{fpropsfluid|butane}},
{{fpropsfluid|butene}},
{{fpropsfluid|carbondioxide}},
{{fpropsfluid|carbonmonoxide}},
{{fpropsfluid|carbonylsulfide}},
{{fpropsfluid|cisbutene}},
{{fpropsfluid|decane}},
{{fpropsfluid|dimethylether}},
{{fpropsfluid|ethane}},
{{fpropsfluid|ethanol}},
{{fpropsfluid|hydrogen}},
{{fpropsfluid|hydrogensulfide}},
{{fpropsfluid|isobutane}},
{{fpropsfluid|isobutene}},
{{fpropsfluid|isohexane}},
{{fpropsfluid|isopentane}},
{{fpropsfluid|krypton}},
{{fpropsfluid|methane}},
{{fpropsfluid|neopentane}},
{{fpropsfluid|nitrogen}},
{{fpropsfluid|nitrousoxide}},
{{fpropsfluid|nonane}},
{{fpropsfluid|parahydrogen}},
{{fpropsfluid|propane}},
{{fpropsfluid|r116}},
{{fpropsfluid|r134a}},
{{fpropsfluid|r141b}},
{{fpropsfluid|r142b}},
{{fpropsfluid|r218}},
{{fpropsfluid|r245fa}},
{{fpropsfluid|r41}},
{{fpropsfluid|sulfurdioxide}},
{{fpropsfluid|toluene}},
{{fpropsfluid|transbutene}},
{{fpropsfluid|water}},
{{fpropsfluid|xenon}}

<small>408 cubic EOS fluids:</small> {{fpropsfluid|_rpp}}
</div>
'''FPROPS''' is a free open-source C-based library for high-accuracy evaluation of thermodynamic properties for a number of pure substances. It makes use of published data for the Helmholtz fundamental equation for those substances. It has been developed by [[User:Jpye|John Pye]] and [[#Contributors|others]], 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 on the right. The properties that can be calculated are internal energy <math>u</math>, entropy <math>s</math>, pressure <math>p</math>, enthalpy <math>h</math> and Helmholtz energy <math>a</math>, 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.cfm REFPROP], [http://www.mech.kyushu-u.ac.jp/~heat/propath/ PROPATH], EES, [http://www.thermo.ruhr-uni-bochum.de/en/prof-w-wagner/software/fluidcal.html FLUIDCAL], SteamTab and others, but is free open-source software, licensed under the [http://www.gnu.org/licenses/gpl.html GPL]. Similar open-source libraries include [[freesteam]] and [http://coolprop.sourceforge.net CoolProp]. Contributions are most welcome!

Enhancement to FPROPS feature set and list of supported fluids has been going on as part of ASCEND's involvement in [[GSOC]].

FPROPS is included in current versions of ASCEND, but you can also browse the source code at {{srcdir|models/johnpye/fprops/}} (current trunk) <!-- or {{srcbranchdir|fprops2|models/johnpye/fprops/}}-->. Currently (Mar 2014), [[FPROPS/Viscosity|viscosity]] and [[FPROPS/Thermal conductivity|thermal conductivity]] correlations are being implemented in FPROPS.

== Implementation ==

The calculation routines implement calculation of reduced Helmholtz energy <math>\phi = a / R T</math> as a function of reduced density <math>\delta = \rho / \rho_c</math> and ''inverse'' reduced temperature <math>\tau = T_c / T</math>, by dividing it into ideal <math>\phi^0 \left(\tau,\delta\right)</math> and 'real' or 'residual' <math>\phi^r \left(\tau,\delta \right)</math> parts, as described in IAPWS-95<ref name=iapws95>http://www.iapws.org/relguide/IAPWS95-Rev.pdf</ref>.

=== 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.

An expression for <math>c_p^0</math> can then be used to calculate ideal enthalpy, ideal entropy and ideal helmholtz energy as follows <ref name="lemmon20fluids">E W Lemmon and R Span, 2006. Short Fundamental Equations of State for 20 Industrial Fluids, ''J. Chem. Eng. Data'' '''51''', {{doi|10.1021/je050186n}}</ref>:

:<math>h^0 = h_0^0 + \int_{T_0}^{T}{c_p^0 dT}</math>

:<math>s^0 = s_0^0 + \int_{T_0}^{T}{\frac{c_p^0}{T}dT} - R \ln \left(\frac{\rho T}{\rho_0 T_0}\right)</math> 

:<math>a^0 = h^0  - RT - T s^0</math>

Combining these, an expression for calculating the nondimensionalised ideal component of Helmholtz energy <math>a^0 \left(T , \rho \right)</math> from <math>c_p^0</math> can be found, which incorporates additional correlation constants relating to the reference state <math>\left(T_0, \rho_0, h_0^0, s_0^0 \right)</math>: <ref>Span & Wagner (1996), "A new equation of state for carbon dioxide covering from the triple point temperature to 1100 K at pressures up to 800 MPa" ''J. Phys. Chem. Ref. Data'', '''25''' p 1509 ff, http://link.aip.org/link/?jpr/25/1509%26agg=NIST, {{doi|10.1063/1.555991}}</ref>

:<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>.

In the reverse direction, using the relations for the Helmholtz fundamental equation, an important result is

:<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). 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 are, in that case, determined using specified values of enthalpy <math>h^0_{ref}</math> and entropy <math>s^0_{ref}</math> at a reference state <math>\left(T_{ref},\rho_{ref} \right)</math>

In FPROPS, then, we permit <math>\phi^0</math> expressions of this form:

:<math>\phi^0 = \ln \left(\delta \right) + c + m \tau + \sum_{i=0}^{\mathrm{np-1}} {\begin{cases}a_i \ln \tau &\mbox{if }p_i =0 \\ a_i \tau^p &\mbox{otherwise}\end{cases}} + \sum_0^{\mathrm{ne}-1} {n_i \ln \left[1 - \exp\left(-\gamma_i \tau\right)\right]}</math>

This corresponds, after the above transformation, to:

:<math>\frac{c_p^0}{R} = 1 + \sum_{i=0}^{\mathrm{np}-1} {c_i T^{t_i}} + \sum_{i=0}^{\mathrm{ne}-1} \frac{b_i x_i^2 \exp\left(-x_i\right)}{[1-\exp\left(-x_i\right)]^2} \mbox{  where }x_i = \frac{\beta_i}{T}</math>

The conversion from <math>\phi^0\left(\tau,\delta\right)</math> terms to <math>c_p^0</math> terms is as follows (<math>\tau = T^{*} / T</math>): 

{| border="1" cellpadding="3" cellspacing="0"
! type !! term in <math>\phi^0</math> !! term in <math>c_p^0</math> !! conversion
|-
| power terms || <math>\begin{cases} a \ln \tau, & \mbox{if } p = 0 \\
a \tau^p, & \mbox {otherwise} \end{cases}</math>
|| <math>c T^{t}</math> || <math>c = \begin{cases}a, & \mbox{if } p = 0 \\ - a p (p-1) \left(T^{*}\right)^p & \mbox{otherwise}\end{cases}</math><br><math>t = -p</math>
|-
| Planck-Einstein terms || <math>n_i \ln \left[1 - e^{-\gamma_i \tau}\right]</math> || <math>\frac{b_i x_i^2 \exp\left(-x_i\right)}{[1-\exp\left(-x_i\right)]^2}</math> where <math>x_i = \frac{\beta_i}{T}</math> || <math>b = n</math><br><math>\beta = T^{*}\gamma</math>
|-
| 'lead' terms || <math>\ln\left(\delta\right) + c + m \tau</math> || <math>1</math> || none required
|}

Note the factor '1' in the lead terms is not usually specified separately in publications; our code compensates for that when reading in <math>c_p^0</math> data; you don't need to subtract 1 from the provided constant term.

If a publication provides <math>\phi^0</math>, then the above allows calculation of <math>c_p^0</math>.

If on the other hand a publication provides <math>c_p^0</math> then we need to use the reference state properties <math>\left(T_0,\rho_0,h_0^0, s_0^0\right)</math>, to define <math>c</math> and <math>m</math>, as follows:

:<math>m = \frac{h_0^0}{R T_c}</math>
:<math>c = -\frac{s_0^0}{R} - 1 - \ln\delta_0 + ln\tau_0</math>

where <math>\tau_0 =  T^{*}/T_0</math> and <math>\delta_0 = \rho_0 / \rho^{*}</math>.

However, if we are using this ideal correlation together with a 'real' part, then these constants <math>c</math> and <math>m</math> will instead take values determined from reference state properties <math>\left(T_0,\rho_0,h_0, s_0\right)</math>, meaning that these expressions will be offset as a result of <math>h_0^0 - h_0</math> and <math>s_0^0 - s_0</math> variations.

=== Residual part ===

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

{| border="1" cellpadding="3" cellspacing="0"
! type !! expression !! correlation constants
|-
| power terms || <math>

\begin{cases} 
  a_i \tau^{t_i} \delta^{d_i} \exp\left(-\delta^{l_i}\right), & \mbox{if }l \ne 0 \\
  a_i \tau^{t_i} \delta^{d_i},  & \mbox{otherwise} 
\end{cases}


</math> || <math>a, t, d, l</math>
|-
| gaussian terms || <math>n_i \tau^{t_i} \delta^{d_i} \exp\left[-\alpha_i(\delta - \epsilon_i)^2 - \beta_i(\tau - \gamma_i)^2\right]</math> || <math> n, t, d, \alpha, \beta, \gamma, \epsilon</math>
|-
| critical terms || <math>n_i \Delta^{b_i} \delta \psi \mbox{  where } \begin{cases}
  \Delta = \theta^2 + B_i \left[ \left( \delta - 1 \right)^2 \right]^{a_i} \\
  \theta = \left( 1- \tau \right) + A_i \left[ \left( \delta - 1 \right)^2 \right]^{\frac{1}{2 \beta_i}} \\
  \psi = \exp \left[-C_i \left(\delta-1 \right)^2 - D_i \left( \tau - 1 \right)^2 \right]
\end{cases}</math>
|| <math>n, a, b, \beta, A, B, C, D</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> and hence the specific Helmholtz energy <math>a = R T \phi \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<ref name=iapws95></ref> as well as by Tillner Roth et al<ref>Tillner-Roth, Harms-Watzenberg and Baehr (1993), "Eine neue Fundamentalgleichung für Ammoniak",
''DKV-Tagungsbericht'', '''20''', 167-181</ref>.

=== Saturated mixture properties ===

The Helmholtz function <math>a \left(\rho,T \right)</math>, given above, applies for all regions outside the saturation region. The saturation region, on the other hand, is not calculated directly from the raw Helmholtz correlation. Instead within the saturation region, the stable equilibrium of phases must be calculated such that the liquid phase vapour phases have equal Gibbs free energy <math>g = a + p / \rho</math>. This requirement is equivalent to the Maxwell equal area rule<ref>http://en.wikipedia.org/wiki/Van_der_Waals_equation</ref><ref name=iapws95></ref><ref>Eubank & Hall (2004), "Equal area rule and algorithm for determining phase compositions", ''AIChE Journal'', '''41''', {{doi|10.1002/aic.690410419}}</ref>. The condition can be expressed as a set of three simultaneous equations in the variables <math>p_{sat}, \rho_f, \rho_g</math>:

:<math>
\begin{align}
p_{sat} & = p \left(\rho_f, T \right) \\
p_{sat} & = p \left(\rho_g, T \right) \\
g\left(\rho_f,T \right) &= g\left(\rho_g,T \right)
\end{align}
</math>

This can be further reduced to :

:<math>
\begin{align}
p_{sat} & = R T \rho_f \left(1 + \delta_f \phi^r_\delta \left(\tau,\delta_f \right) \right) \\
p_{sat} & = R T \rho_g \left(1 + \delta_g \phi^r_\delta \left(\tau,\delta_g, \right) \right) \\
\frac{p_{sat}}{R T} \left( \frac{1}{\rho_g} - \frac{1}{\rho_f} \right) & = \phi \left(\tau, \delta_f \right) - \phi \left(\tau, \delta_g \right)
 = \ln \left( \frac{\rho_f}{\rho_g} \right) + \phi^r \left(\tau, \delta_f \right) - \phi^r \left(\tau, \delta_g  \right)
\end{align}
</math>

Calculation of the phase equilibrium conditions requires iterative solution of these three equations; one suitable approach is that of Akasaka<ref name=akasaka2008>Akasaka (2008) "A Reliable and Useful Method to Determine the Saturation State from Helmholtz Energy Equations of State", ''Journal of Thermal Science and Technology'', '''3''', 442-451, {{doi|10.1299/jtst.3.442}}</ref>, which first eliminates <math>p_{sat}</math> from the above equations, then iterates only on <math>\rho_f</math> and <math>\rho_g</math>. The approach seems to be robust, with tests currently performed with a sample of five fluids including water, carbon-dioxide and nitrogen. Our code is in {{src|models/johnpye/fprops/sat.c}}.

The Akasaka approach is markedly simpler to implement that the complicated algorithm used in REFPROP (which is explained by Akasaka, and also by Soave<ref>Soave (1986), "Direct calculation of pure-compound vapour pressures through cubic equations of state ", ''Fluid Phase Equilibria'', '''31''', 203-207, {{doi|10.1016/0378-3812(86)90013-0}}</ref>). A limitation of the Soave approach is that it does not converge for temperatures <math>T \gtrapprox 0.99 T_c</math>, and hence requires yet another, even more complicated, algorithm to be implemented for that region (including need for the very tedious-to-calculate third partial derivatives <math>\phi^r_{\delta \delta \delta}</math>!).

=== Saturation Tsat(p) ===

To calculation the saturation temperature as a function of pressure, a different iteration scheme from the above for <math>p_{sat}(T)</math> is needed. <ref>Long X. Nghiem and Yau-Kun Li, 1985, ''Application of the tangent plane criterion to saturation pressure and temperature computations'', Fluid Phase Equilibria, {{doi|10.1016/0378-3812(85)90059-7}}</ref><ref>M Michelson, 1984, ''Saturation point calculations'', Fluid Phase Equilibria, {{doi|10.1016/0378-3812(85)90005-6}}</ref>

At present, we have implemented <tt>fprops_sat_p</tt> in {{src|models/johnpye/fprops/sat.c}} which provides a solution via simple Newton iteration that makes repeated calls to <tt>fprops_sat_T</tt>. The iteration makes use of the Clapeyron equation to calculate the gradients,
:<math>\left(\frac{\partial p}{\partial T} \right)_{sat} = \frac{h_{fg}}{T v_{fg}}</math>
There is probably a more efficient approach, but this is adequate for the moment.

=== Thermodynamic property derivatives ===
{{experimental}}

The IAPWS has made a release<ref>IAPWS, 2008, ''Revised Advisory Note No. 3 Thermodynamic Derivatives from IAPWS Formulation'', http://www.iapws.org</ref> with information about how partial derivatives of thermodynamic properties with respect to other properties can be calculated in general. 

FPROPS includes an initial implementation of these derivative calculations in {{src|models/johnpye/fprops/derivs.c}}, but they have not yet been thoroughly tested.

=== Properties in terms of (p,h) ===

{{experimental}}
It seems to be very helpful to have property correlations calculable in terms of the independent variables (p,h). This seems to aid stability of calculation in larger systems, compared to the (v,T) property pair. This would appear to be in part due to the very small, steep subcooled region when using (v,T) coordinates.

An initial implementation of an iterative solver for properties in terms of (p,h) has been implemented for FPROPS in {{src|models/johnpye/fprops/solve_ph.c}}. Currently, it solves for all regions except for pressures below the triple-point pressure. We hope to overcome this problem soon.

=== Properties in terms of (T,x) ===

For convenience some simple functions for calculating properties in terms of <math>(T,x)</math> have been written. The implementation is in {{src|models/johnpye/fprops/solve_Tx.c}}. Any input with <math>x</math> outside <math>[0,1]</math> or <math>T</math> outside <math>\left[ T_t, T_c \right]</math> will return an error code. Otherwise, the corresponding rho (<math>\rho</math>) value will be returned such that <math>(T,\rho)</math> is equivalent to the desired state.

A python demonstration of this code is in {{src|models/johnpye/fprops/python/solve_Tx.py}}, which gives the following output:

<pre>
john@thunder:~/ascend/models/johnpye/fprops/python$ python solve_Tx.py 
solved state
T = 300.000000, x = 0.500000
--> rho = 0.309802, h = 395476.135716, s = 381991.898474, u = 1491.770280

checking reverse solve:
--> x = 0.500000
</pre>

=== Transport properties ===

Initial implementations of [[FPROPS/Viscosity|dynamic viscosity]] and [[FPROPS/Thermal conductivity|thermal conductivity]] have been created for carbon dioxide. Providing the fluid data supports this, you should be able to access these properties from Python, C and ASCEND. In Python,

<source lang="a4c">
from fprops import *
D = fluid("carbondioxide")
S = D.set_Trho(220,2.5)
print "mu =",S.mu
print "lam =",S.lam
</source>

Both properties are returned in base SI units (Pa*s for viscosity, and W/m·K for thermal conductvity).

== 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.cfm 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&lt;/tdd&gt;) then run the test script &lt;tt&gt;./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 (should be) embedded at the end of the corresponding fluid file, for example {{src|models/johnpye/fprops/fluids/hydrogen.c}}.

During development, some additional testing routines are in the python subdirectory. These include

* <tt>satcvgc.py</tt> to test converge of saturation region routines <tt>fprops_sat_T</tt> and <tt>fprops_sat_p</tt> are working.
* <tt>solve_ph_array</tt> to test convergence of the <tt>fprops_solve_ph</tt> routine, including plotting of regions where the routine is failing.

== 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/fluids/hydrogen.c}}.

After adding the data structures in "component".c file, a corresponding "component".h is required. 

An example of which is shown in {{src|models/johnpye/fprops/fluids/hydrogen.h}}

The new component has to be declared in SConstruct and asc_helmholtz.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 <math>h</math> and <math>s</math>, these formulae tell you the offset to apply to the constant and linear terms in <math>\phi^0(\tau,\delta)</math>.

== Using within ASCEND ==

''Note that the interface for accessing FPROPS is being changed significantly in 2012, in the 'fprops' branch code and work by Sean Muratet.''

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>, to calculate p(T,rho)
* <tt>helmholtz_u</tt>, to calculate u(T,rho)
* <tt>helmholtz_h</tt>, to calculate h(T,rho)
* <tt>helmholtz_s</tt>, to calculate s(T,rho)
* <tt>helmholtz_a</tt>, to calculate a(T,rho)
* <tt>helmholtz_cp</tt>, to calculate cp(T,rho)
* <tt>helmholtz_cv</tt>, to calculate cv(T,rho)
* <tt>helmholtz_w</tt>, to calculate speed of sound for given (T,rho)
* <tt>helmholtz_phsx_vT</tt> is also available which returns p, h, s and x in for specified (v,T).
* <tt>helmholtz_Tvsx_ph</tt> returns T, v, s and x for specified (p,h). Recommended for energy system models.

'''Note'''. When solving properties in 'reverse' mode (eg solving for <span class="texhtml">&rho;,''T''</span> when <span class="texhtml">''p'',''h''</span> 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. Also, use <tt>OPTION iterationlimit 200</tt> for larger models.

A first system model using FPROPS is available in {{src|models/johnpye/fprops/rankine_fprops.a4c}}. The key part is
<source lang="a4c">
cd IS_A fluid;
cd.component :== 'water';
(* ...definitions of p, h, T, v, s, x omitted here... *)

calc_ph: helmholtz_Tvsx_ph(
    p, h : INPUT;
    T, v, s, x : OUTPUT;
    cd : DATA
);
</source>

Some details on the implementation of the 'wrapper' code for connecting FPROPS with ASCEND is given in [[writing ASCEND external relations in C]].

== Using from C ==

Using the FPROPS code from C should also be very straightforward once the relevant header files have been examined. Warning! We are still working on the API and it is subject to change, especially in the light of current [[GSOC2011]] work.

A very simple test of selecting a fluid, then solving for (p,h) values and then returning results, is given in {{src|models/johnpye/fprops/test_ph.c}} (''Note: this file is currently missing and needs to be restored.''). An extract of that file is shown below:

<source lang=c>
double p = 1e5;
double h = 3500e3;const char *fluid = "water";
double T = 0,rho = 0;
const HelmholtzData *D = fprops_fluid(fluid);
int res = fprops_solve_ph(p,h,&T, &rho, 0, D);
if(res)printf("ERROR %d!)\n",res);
printf("Result: T = %f K, rho = %f kg/m3\n",T,rho);
double p1, h1;
p1 = helmholtz_p(T,rho, D);
h1 = helmholtz_h(T,rho, D);
printf("Double check: h = %f kJ/kg\n",h1/1e3);
</source>

Other test code written in C is shown in {{src|models/johnpye/fprops/test.c}}, which may be helpful for seeing some other features in use.

At this stage no external documentation for the use of the C code has been written, other that what is on this page.

== Using from Python ==

FPROPS can also be used from the Python language. There are a number of sample programs in {{srcdir|models/johnpye/fprops/python}}. The API is still evolving as we seek the cleanest ways to report errors back to the user, and catch floating-point errors, etc.

== Further work ==
{{task}}

At the conclusion of [[GSOC2012]] with [[Sean Muratet]], we completed implementation of a cubic EOS (Peng Robinson) as well as heavily-rewritten internal fluids database and data structures. Work continues in Sean's branch as well as the {{srcbranchdir|fprops2|models/johnpye/fprops}} branch.

At the end of [[GSOC2011]] with [[Richard Towers]], we concluded the following steps are required using his {{srcbranchdir|richard|}} branch as the basis. Further work is now (Aug 2011) goin on in the {{srcbranchdir|fprops2|models/johnpye/fprops}} branch.

# <s>checking/evaluatingthe proposed data structure (is the design there right?) and modifying if required</s>
# fixing up broken ASCEND bindings
# <s>adding more fluid data, eg from critical point/omega data from RPP</s>
# adding test code for PR EOS-based evaluations
# fixing up the PR EOS evaluation code until the test code works
# merging to trunk
# <s>extending python bindings to raise exceptions instead of returning error flags (nicer API design)</s>
# performance testing, profiling, etc.

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

* Adding incompressible fluids EOS
* Adding additional cubic EOS
* Adding option for interpolation of saturation curve (as seen in CoolProp)
* Extend test suite with more comparison of RPP data and cubic EOS code
* Revisit derivatives calculation for new EOSes.
* Expose derivatives code to ASCEND and Python; extend ASCEND to support optimisation with [[external relations]].
* <s>add support for '''calculation of the saturation curve''' using Maxwell phase-equilibrium condition</s>
* <s>add iterative solver for properties in terms of (p,h)</s>
* <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.
* <s>move specification of the '''zero point''' out of IdealData and into HelmholtzData. Third law of thermodynamics? How to avoid correlations with negative ''h'' or ''s'' values?</s>
* <s>implement a '''pre-calculation routine''' where for example coefficients of <math>\phi^0</math> can be prepared?</s>
* implement support for thermophysical properties '''conductivity, viscosity, surface tension,''' and possibly others.
* <s>implement support for calculation of the speed of sound.</s>
* <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?</s>
* <s>remove p_c from HelmholtzData structure, because that value can be calculated using p(T_c,rho_c).</s>
* add support for solving in terms of (p,h) where p < p_t, the triple point pressure.
* possibly some additional coefficients needs in order to support correlations such as for Methanol<ref>K de Reuck, 1993, ''Methanol, International Thermodynamic Tables of the Fluid State Vol. 12 (1993)'', Blackwell/IUPAC, http://www.iupac.org/objID/Source/sou65427568843800267934799</ref>
* <s>add properties of R-134a (see {{doi|10.1063/1.555958}})</s>

Further down the track we would like to support:

* reading property correlations from a '''database''' (SQLite?) or text file.
* <font color=green>some simpler equations of state such as MBWR or Peng-Robinson.</font>(Assigned to Ankit, see [[PengRobinson EOS in FPROPS]])
* lots more fluids
* '''mixing models''' for multi-phase, multi-component systems (or else incorporate FPROPS into exist code that supports this).
* [http://www.gerg.info/publications/tm/tm15_04.pdf natural gas properties?]
* [http://www.if.uidaho.edu/~vutgikar/Fall2007/Hydrogen/AIR.pdf air properties] {{doi|10.1063/1.1285884}}

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, subject to the restrictions of the GPL license.

== References ==

<references/>

== Contributors ==

The people who have contributed code and/or implementable advice for FPROPS:
* [[User:Jpye|John Pye]]
* [[User:Smuratet|Sean Muratet]]
* [[User:Kchittur|Krishnan Chittur]]
* Hong Ke
* [[Richard Towers]]
* [[User:Rshrikanth|Shrikanth Ranganadham ]]
* [[User:Ankitml|Ankit Mittal]]

== See Also ==

* [[Combined-cycle gas turbine]]
* [[Thermodynamics with ASCEND]]
* [[freesteam]]

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