Cantera  2.1.2
WaterPropsIAPWS.h
Go to the documentation of this file.
1 /**
2  * @file WaterPropsIAPWS.h
3  * Headers for a class for calculating the equation of state of water
4  * from the IAPWS 1995 Formulation based on the steam tables thermodynamic
6  */
7 /*
8  * Copyright (2005) Sandia Corporation. Under the terms of
9  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
10  * U.S. Government retains certain rights in this software.
11  */
12 #ifndef WATERPROPSIAPWS_H
13 #define WATERPROPSIAPWS_H
14
15 #include "WaterPropsIAPWSphi.h"
16 #include "cantera/base/config.h"
17
18 namespace Cantera
19 {
20 /**
21  * @name Names for the phase regions
22  *
23  * These constants are defined and used in the interface
24  * to describe the location of where we are in (T,rho) space.
25  *
26  * WATER_UNSTABLELIQUID indicates that we are in the unstable region, inside the
27  * spinodal curve where dpdrho < 0.0 amonst other properties. The difference
28  * between WATER_UNSTABLELIQUID and WATER_UNSTABLEGAS is that
29  * for WATER_UNSTABLELIQUID d2pdrho2 > 0 and dpdrho < 0.0
30  * for WATER_UNSTABLEGAS d2pdrho2 < 0 and dpdrho < 0.0
31  */
32 //@{
33 #define WATER_GAS 0
34 #define WATER_LIQUID 1
35 #define WATER_SUPERCRIT 2
36 #define WATER_UNSTABLELIQUID 3
37 #define WATER_UNSTABLEGAS 4
38 //@}
39
40 //! Class for calculating the equation of state of water.
41 /*!
42  * The reference is W. Wagner, A. Prub, "The IAPWS Formulation 1995 for the
43  * Thermodynamic Properties of Ordinary Water Substance for General and
44  * Scientific Use," J. Phys. Chem. Ref. Dat, 31, 387, 2002.
45  *
46  * This class provides a very complicated polynomial for the specific
47  * Helmholtz free energy of water, as a function of temperature and density.
48  *
49  * \f[
50  * \frac{M\hat{f}(\rho,T)}{R T} = \phi(\delta, \tau) =
51  * \phi^o(\delta, \tau) + \phi^r(\delta, \tau)
52  * \f]
53  *
54  * where
55  *
56  * \f[
57  * \delta = \rho / \rho_c \quad \mathrm{and} \quad \tau = T_c / T
58  * \f]
59  *
60  * The following constants are assumed
61  *
62  * \f[
63  * T_c = 647.096\mathrm{\;K}
64  * \f]
65  * \f[
66  * \rho_c = 322 \mathrm{\;kg\,m^{-3}}
67  * \f]
68  * \f[
69  * R/M = 0.46151805 \mathrm{\;kJ\,kg^{-1}\,K^{-1}}
70  * \f]
71  *
72  * The free energy is a unique single-valued function of the temperature and
73  * density over its entire range.
74  *
75  * Note, the base thermodynamic state for this class is the one used in the
76  * steam tables, i.e., the liquid at the triple point for water has the
77  * following properties:
78  *
79  * - u(273.16, rho) = 0.0
80  * - s(273.16, rho) = 0.0
81  * - psat(273.16) = 611.655 Pascal
82  * - rho(273.16, psat) = 999.793 kg m-3
83  *
84  * Therefore, to use this class within %Cantera, offsets to u() and s() must
85  * be used to put the water class onto the same basis as other thermodynamic
86  * quantities. For example, in the WaterSSTP class, these offsets are
87  * calculated in the following way. The thermodynamic base state for water is
88  * set to the NIST basis here by specifying constants EW_Offset and SW_Offset.
89  * These offsets are calculated on the fly so that the following properties
90  * hold:
91  *
92  * - Delta_Hfo_idealGas(298.15, 1bar) = -241.826 kJ/gmol
93  * - So_idealGas(298.15, 1bar) = 188.835 J/gmolK
94  *
95  * The offsets are calculated by actually computing the above quantities and
96  * then calculating the correction factor.
97  *
98  * This class provides an interface to the WaterPropsIAPWSphi class, which
99  * actually calculates the \f$\phi^o(\delta, \tau) \f$ and the
100  * \f$\phi^r(\delta, \tau) \f$ polynomials in dimensionless form.
101  *
102  * All thermodynamic results from this class are returned in dimensional form.
103  * This is because the gas constant (and molecular weight) used within this
104  * class is allowed to be potentially different than that used elsewhere in
105  * %Cantera. Therefore, everything has to be in dimensional units. Note,
106  * however, the thermodynamic basis is set to that used in the steam tables.
107  * (u = s = 0 for liquid water at the triple point).
108  *
109  * This class is not a %ThermoPhase. However, it does maintain an internal
110  * state of the object that is dependent on temperature and density. The
111  * internal state is characterized by an internally stored \f$\tau\f$ and a
112  * \f$\delta \f$ value, and an iState value, which indicates whether the
113  * point is a liquid, a gas, or a supercritical fluid. Along with that the
114  * \f$\tau\f$ and a \f$\delta \f$ values are polynomials of \f$\tau\f$ and
115  * a \f$\delta \f$ that are kept by the WaterPropsIAPWSphi class. Therefore,
116  * whenever \f$\tau\f$ or \f$\delta \f$ is changed, the function setState()
117  * must be called in order for the internal state to be kept up to date.
118  *
119  * The class is pretty straightforward. However, one function deserves
120  * mention. The density() function calculates the density that is consistent
121  * with a particular value of the temperature and pressure. It may therefore
122  * be multivalued or potentially there may be no answer from this function. It
123  * therefore takes a phase guess and a density guess as optional parameters.
124  * If no guesses are supplied to density(), a gas phase guess is assumed. This
125  * may or may not be what is wanted. Therefore, density() should usually at
126  * least be supplied with a phase guess so that it may manufacture an
127  * appropriate density guess. density() manufactures the initial density
128  * guess, nondimensionalizes everything, and then calls
129  * WaterPropsIAPWSphi::dfind(), which does the iterative calculation to find
130  * the density condition that matches the desired input pressure.
131  *
132  * The phase guess defines are located in the .h file. they are
133  *
134  * - WATER_GAS
135  * - WATER_LIQUID
136  * - WATER_SUPERCRIT
137  *
138  * There are only three functions which actually change the value of the
139  * internal state of this object after it's been instantiated
140  *
141  * - setState_TR(temperature, rho)
142  * - density(temperature, pressure, phase, rhoguess)
143  * - psat(temperature, waterState);
144  *
145  * The setState_TR() is the main function that sets the temperature and rho
146  * value. The density() function serves as a setState_TP() function, in that
147  * it sets internal state to a temperature and pressure. However, note that
148  * this is potentially multivalued. Therefore, we need to supply in addition a
149  * phase guess and a rho guess to the input temperature and pressure. The
150  * psat() function sets the internal state to the saturated liquid or
151  * saturated gas state, depending on the waterState parameter.
152  *
153  * Because the underlying object WaterPropsIAPWSphi is privately held, you can
154  * be sure that the underlying state of this object doesn't change except due
155  * to the three function calls listed above.
156  *
157  * @ingroup thermoprops
158  */
160 {
161 public:
162  //! Base constructor
163  WaterPropsIAPWS();
164
165  //! Copy constructor
166  WaterPropsIAPWS(const WaterPropsIAPWS& right);
167
168  //! assignment constructor
170
171  //! destructor
173
174  //! Set the internal state of the object wrt temperature and density
175  /*!
176  * @param temperature temperature (kelvin)
177  * @param rho density (kg m-3)
178  */
179  void setState_TR(doublereal temperature, doublereal rho);
180
181  //! Calculate the Helmholtz free energy in mks units of J kmol-1 K-1,
182  //! using the last temperature and density
183  doublereal helmholtzFE() const;
184
185  //! Calculate the Gibbs free energy in mks units of J kmol-1 K-1.
186  //! using the last temperature and density
187  doublereal Gibbs() const;
188
189  //! Calculate the enthalpy in mks units of J kmol-1
190  //! using the last temperature and density
191  doublereal enthalpy() const;
192
193  //! Calculate the internal energy in mks units of J kmol-1
194  doublereal intEnergy() const;
195
196  //! Calculate the entropy in mks units of J kmol-1 K-1
197  doublereal entropy() const;
198
199  //! Calculate the constant volume heat capacity in mks units of J kmol-1 K-1
200  //! at the last temperature and density
201  doublereal cv() const;
202
203  //! Calculate the constant pressure heat capacity in mks units of J kmol-1 K-1
204  //! at the last temperature and density
205  doublereal cp() const;
206
207  //! Calculate the molar volume (kmol m-3)
208  //! at the last temperature and density
209  doublereal molarVolume() const;
210
211  //! Calculates the pressure (Pascals), given the current value of the
212  //! temperature and density.
213  /*!
214  * The density is an independent variable in the underlying equation of state
215  *
216  * @return
217  * returns the pressure (Pascal)
218  */
219  doublereal pressure() const;
220
221  //! Calculates the density given the temperature and the pressure,
222  //! and a guess at the density. Sets the internal state.
223  /*!
224  * Note, below T_c, this is a multivalued function.
225  *
226  * The density() function calculates the density that is consistent with
227  * a particular value of the temperature and pressure. It may therefore be
228  * multivalued or potentially there may be no answer from this function.
229  * It therefore takes a phase guess and a density guess as optional
230  * parameters. If no guesses are supplied to density(), a gas phase guess
231  * is assumed. This may or may not be what is wanted. Therefore, density()
232  * should usually at least be supplied with a phase guess so that it may
233  * manufacture an appropriate density guess. density() manufactures the
234  * initial density guess, nondimensionalizes everything, and then calls
235  * WaterPropsIAPWSphi::dfind(), which does the iterative calculation to
236  * find the density condition that matches the desired input pressure.
237  *
238  * @param temperature: Kelvin
239  * @param pressure : Pressure in Pascals (Newton/m**2)
240  * @param phase : guessed phase of water
241  * : -1: no guessed phase
242  * @param rhoguess : guessed density of the water
243  * : -1.0 no guessed density
244  * @return
245  * Returns the density. If an error is encountered in the calculation
246  * the value of -1.0 is returned.
247  */
248  doublereal density(doublereal temperature, doublereal pressure,
249  int phase = -1, doublereal rhoguess = -1.0);
250
251  //! Calculates the density given the temperature and the pressure,
252  //! and a guess at the density, while not changing the internal state
253  /*!
254  * Note, below T_c, this is a multivalued function.
255  *
256  * The density() function calculates the density that is consistent with a
257  * particular value of the temperature and pressure. It may therefore be
258  * multivalued or potentially there may be no answer from this function.
259  * It therefore takes a phase guess and a density guess as optional
260  * parameters. If no guesses are supplied to density(), a gas phase guess
261  * is assumed. This may or may not be what is wanted. Therefore, density()
262  * should usually at least be supplied with a phase guess so that it may
263  * manufacture an appropriate density guess. density() manufactures the
264  * initial density guess, nondimensionalizes everything, and then calls
265  * WaterPropsIAPWSphi::dfind(), which does the iterative calculation to
266  * find the density condition that matches the desired input pressure.
267  *
268  * @param pressure : Pressure in Pascals (Newton/m**2)
269  * @param phase : guessed phase of water
270  * : -1: no guessed phase
271  * @param rhoguess : guessed density of the water
272  * : -1.0 no guessed density
273  * @return
274  * Returns the density. If an error is encountered in the calculation
275  * the value of -1.0 is returned.
276  */
277  doublereal density_const(doublereal pressure, int phase = -1, doublereal rhoguess = -1.0) const;
278
279  //! Returns the density (kg m-3)
280  /*!
281  * The density is an independent variable in the underlying equation of state
282  *
283  * @return Returns the density (kg m-3)
284  */
285  doublereal density() const;
286
287  //! Returns the temperature (Kelvin)
288  /*!
289  * @return Returns the internally stored temperature
290  */
291  doublereal temperature() const;
292
293  //! Returns the coefficient of thermal expansion.
294  /*!
295  * alpha = d (ln V) / dT at constant P.
296  *
297  * @return
298  * Returns the coefficient of thermal expansion
299  */
300  doublereal coeffThermExp() const;
301
302  //! Returns the isochoric pressure derivative wrt temperature
303  /*!
304  * beta = M / (rho * Rgas) (d (pressure) / dT) at constant rho
305  *
306  * Note for ideal gases this is equal to one.
307  *
308  * beta = delta (phi0_d() + phiR_d())
309  * - tau delta (phi0_dt() + phiR_dt())
310  */
311  doublereal coeffPresExp() const;
312
313  //! Returns the coefficient of isothermal compressibility for the
314  //! state of the object
315  /*!
316  * kappa = - d (ln V) / dP at constant T.
317  *
318  * units - 1/Pascal
319  *
320  * @return
321  * returns the isothermal compressibility
322  */
323  doublereal isothermalCompressibility() const;
324
325  //! Returns the value of dp / drho at constant T for the
326  //! state of the object
327  /*!
328  * units - Joules / kg
329  *
330  * @return
331  * returns dpdrho
332  */
333  doublereal dpdrho() const;
334
335  //! This function returns an estimated value for the saturation pressure.
336  /*!
337  * It does this via a polynomial fit of the vapor pressure curve.
338  * units = (Pascals)
339  *
340  * @param temperature Input temperature (Kelvin)
341  *
342  * @return
343  * Returns the estimated saturation pressure
344  */
345  doublereal psat_est(doublereal temperature) const;
346
347  //! This function returns the saturation pressure given the temperature as
348  //! an input parameter, and sets the internal state to the saturated
349  //! conditions.
350  /*!
351  * Note this function will return the saturation pressure, given the
352  * temperature. It will then set the state of the system to the
353  * saturation condition. The input parameter waterState is used to either
354  * specify the liquid state or the gas state at the desired temperature
355  * and saturated pressure.
356  *
357  * If the input temperature, T, is above T_c, this routine will set the
358  * internal state to T and the pressure to P_c. Then, return P_c.
359  *
360  * @param temperature input temperature (kelvin)
361  * @param waterState integer specifying the water state
362  *
363  * @return Returns the saturation pressure. units = Pascal
364  */
365  doublereal psat(doublereal temperature, int waterState = WATER_LIQUID);
366
367  //! Return the value of the density at the water spinodal point (on the liquid side)
368  //! for the current temperature.
369  /*!
370  * @return returns the density with units of kg m-3
371  */
372  doublereal densSpinodalWater() const;
373
374  //! Return the value of the density at the water spinodal point (on the gas side)
375  //! for the current temperature.
376  /*!
377  * @return returns the density with units of kg m-3
378  */
379  doublereal densSpinodalSteam() const;
380
381  //! Returns the Phase State flag for the current state of the object
382  /*!
383  * @param checkState If true, this function does a complete check to see
384  * where in parameters space we are
385  *
386  * There are three values:
387  * - WATER_GAS below the critical temperature but below the critical density
388  * - WATER_LIQUID below the critical temperature but above the critical density
389  * - WATER_SUPERCRIT above the critical temperature
390  */
391  int phaseState(bool checkState = false) const ;
392
393  //! Returns the critical temperature of water (Kelvin)
394  /*!
395  * This is hard coded to the value 647.096 Kelvin
396  */
397  doublereal Tcrit() const {
398  return 647.096;
399  }
400
401  //! Returns the critical pressure of water (22.064E6 Pa)
402  /*!
403  * This is hard coded to the value of 22.064E6 pascals
404  */
405  doublereal Pcrit() const {
406  return 22.064E6;
407  }
408
409  //! Return the critical density of water (kg m-3)
410  /*!
411  * This is equal to 322 kg m-3.
412  */
413  doublereal Rhocrit() const {
414  return 322.;
415  }
416
417 private:
418  //! Calculate the dimensionless temp and rho and store internally.
419  /*!
420  * @param temperature input temperature (kelvin)
421  * @param rho density in kg m-3
422  */
423  void calcDim(doublereal temperature, doublereal rho);
424
425  //! Utility routine in the calculation of the saturation pressure
426  /*!
427  * Calculate the Gibbs free energy in mks units of J kmol-1 K-1.
428  *
429  * @param temperature temperature (kelvin)
430  * @param pressure pressure (Pascal)
431  * @param densLiq Output density of liquid
432  * @param densGas output Density of gas
433  * @param delGRT output delGRT
434  */
435  void corr(doublereal temperature, doublereal pressure, doublereal& densLiq,
436  doublereal& densGas, doublereal& delGRT);
437
438  //! Utility routine in the calculation of the saturation pressure
439  /*!
440  * @param temperature temperature (kelvin)
441  * @param pressure pressure (Pascal)
442  * @param densLiq Output density of liquid
443  * @param densGas output Density of gas
444  * @param pcorr output corrected pressure
445  */
446  void corr1(doublereal temperature, doublereal pressure, doublereal& densLiq,
447  doublereal& densGas, doublereal& pcorr);
448
449  //! pointer to the underlying object that does the calculations.
451
452  //! Dimensionless temperature
453  /*!
454  * tau = T_C / T
455  */
456  doublereal tau;
457
458  //! Dimensionless density
459  /*!
460  * delta = rho / rho_c
461  */
462  mutable doublereal delta;
463
464  //! Current state of the system
465  mutable int iState;
466 };
467
468 }
469 #endif
doublereal densSpinodalSteam() const
Return the value of the density at the water spinodal point (on the gas side) for the current tempera...
doublereal isothermalCompressibility() const
Returns the coefficient of isothermal compressibility for the state of the object.
doublereal psat_est(doublereal temperature) const
This function returns an estimated value for the saturation pressure.
doublereal delta
Dimensionless density.
doublereal Tcrit() const
Returns the critical temperature of water (Kelvin)
void setState_TR(doublereal temperature, doublereal rho)
Set the internal state of the object wrt temperature and density.
Header for Lowest level of the classes which support a real water model (see class WaterPropsIAPWS an...
doublereal cv() const
Calculate the constant volume heat capacity in mks units of J kmol-1 K-1 at the last temperature and ...
int phaseState(bool checkState=false) const
Returns the Phase State flag for the current state of the object.
Class for calculating the equation of state of water.
doublereal density_const(doublereal pressure, int phase=-1, doublereal rhoguess=-1.0) const
Calculates the density given the temperature and the pressure, and a guess at the density...
WaterPropsIAPWS()
Base constructor.
doublereal Gibbs() const
Calculate the Gibbs free energy in mks units of J kmol-1 K-1.
doublereal density() const
Returns the density (kg m-3)
void corr(doublereal temperature, doublereal pressure, doublereal &densLiq, doublereal &densGas, doublereal &delGRT)
Utility routine in the calculation of the saturation pressure.
doublereal coeffThermExp() const
Returns the coefficient of thermal expansion.
doublereal densSpinodalWater() const
Return the value of the density at the water spinodal point (on the liquid side) for the current temp...
WaterPropsIAPWS & operator=(const WaterPropsIAPWS &right)
assignment constructor
void calcDim(doublereal temperature, doublereal rho)
Calculate the dimensionless temp and rho and store internally.
doublereal entropy() const
Calculate the entropy in mks units of J kmol-1 K-1.
int iState
Current state of the system.
doublereal dpdrho() const
Returns the value of dp / drho at constant T for the state of the object.
doublereal molarVolume() const
Calculate the molar volume (kmol m-3) at the last temperature and density.
doublereal enthalpy() const
Calculate the enthalpy in mks units of J kmol-1 using the last temperature and density.
doublereal Rhocrit() const
Return the critical density of water (kg m-3)
doublereal intEnergy() const
Calculate the internal energy in mks units of J kmol-1.
doublereal tau
Dimensionless temperature.
doublereal pressure() const
Calculates the pressure (Pascals), given the current value of the temperature and density...
Low level class for the real description of water.
void corr1(doublereal temperature, doublereal pressure, doublereal &densLiq, doublereal &densGas, doublereal &pcorr)
Utility routine in the calculation of the saturation pressure.
doublereal Pcrit() const
Returns the critical pressure of water (22.064E6 Pa)
doublereal cp() const
Calculate the constant pressure heat capacity in mks units of J kmol-1 K-1 at the last temperature an...
doublereal helmholtzFE() const
Calculate the Helmholtz free energy in mks units of J kmol-1 K-1, using the last temperature and dens...
doublereal temperature() const
Returns the temperature (Kelvin)
WaterPropsIAPWSphi * m_phi
pointer to the underlying object that does the calculations.
doublereal psat(doublereal temperature, int waterState=WATER_LIQUID)
This function returns the saturation pressure given the temperature as an input parameter, and sets the internal state to the saturated conditions.
doublereal coeffPresExp() const
Returns the isochoric pressure derivative wrt temperature.