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