From 871ab9f64cb9739627c31fa91d1bae1d7f71ab36 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Fri, 16 Mar 2007 17:22:38 +0000 Subject: [PATCH] WaterSSTP update Worked on doxygen documentation Changed the XML definition Worked on making it fall in line with other ThermoPhase objects. --- Cantera/src/IdealGasPhase.h | 2 +- Cantera/src/ThermoPhase.h | 4 - Cantera/src/thermo/DebyeHuckel.h | 4 +- Cantera/src/thermo/WaterPropsIAPWS.cpp | 157 +++++--------- Cantera/src/thermo/WaterPropsIAPWS.h | 198 ++++++++++++------ Cantera/src/thermo/WaterPropsIAPWSphi.h | 16 +- Cantera/src/thermo/WaterSSTP.cpp | 47 +++-- Cantera/src/thermo/WaterSSTP.h | 112 ++++++++-- .../cathermo/testWaterTP/output_blessed.txt | 2 - .../cathermo/testWaterTP/testWaterSSTP.cpp | 15 ++ .../cathermo/testWaterTP/waterTPphase.xml | 2 +- 11 files changed, 332 insertions(+), 227 deletions(-) diff --git a/Cantera/src/IdealGasPhase.h b/Cantera/src/IdealGasPhase.h index 6a2b9fdd0..4589ac146 100644 --- a/Cantera/src/IdealGasPhase.h +++ b/Cantera/src/IdealGasPhase.h @@ -278,7 +278,7 @@ namespace Cantera { * XML_Node * const xs = xc->findNameID("phase", "silane"); * IdealGasPhase *silaneGas = new IdealGasPhase(*xs); * @endcode - + * *
*

XML Example

*
diff --git a/Cantera/src/ThermoPhase.h b/Cantera/src/ThermoPhase.h index 08abc04b5..46edd00fd 100755 --- a/Cantera/src/ThermoPhase.h +++ b/Cantera/src/ThermoPhase.h @@ -1459,7 +1459,3 @@ namespace Cantera { #endif - - - - diff --git a/Cantera/src/thermo/DebyeHuckel.h b/Cantera/src/thermo/DebyeHuckel.h index d8783fc50..76fd90faa 100644 --- a/Cantera/src/thermo/DebyeHuckel.h +++ b/Cantera/src/thermo/DebyeHuckel.h @@ -523,7 +523,7 @@ namespace Cantera { * or * * @code - * char iFile[80]; + * char iFile[80], file_ID[80]; * strcpy(iFile, "DH_NaCl.xml"); * sprintf(file_ID,"%s#NaCl_electrolyte", iFile); * XML_Node *xm = get_XML_NameID("phase", file_ID, 0); @@ -533,7 +533,7 @@ namespace Cantera { * or by the following call to importPhase(): * * @code - * char iFile[80]; + * char iFile[80], file_ID[80]; * strcpy(iFile, "DH_NaCl.xml"); * sprintf(file_ID,"%s#NaCl_electrolyte", iFile); * XML_Node *xm = get_XML_NameID("phase", file_ID, 0); diff --git a/Cantera/src/thermo/WaterPropsIAPWS.cpp b/Cantera/src/thermo/WaterPropsIAPWS.cpp index af06a784a..33ffcd083 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWS.cpp @@ -25,6 +25,8 @@ static const double M_water = 18.015268; // kg kmol-1 /* * Note, this is the Rgas value quoted in the paper. For consistency * we have to use that value and not the updated value + * + * The Ratio of R/M = 0.46151805 kJ kg-1 K-1 , which is Eqn. (6.3) in the paper. */ //static const double Rgas = 8.314472E3; // Joules kmol-1 K-1 static const double Rgas = 8.314371E3; // Joules kmol-1 K-1 @@ -67,15 +69,18 @@ WaterPropsIAPWS::~WaterPropsIAPWS() { void WaterPropsIAPWS::calcDim(double temperature, double rho) { tau = T_c / temperature; delta = rho / Rho_c; -} - -/** - * Calculate the Helmholtz Free energy in dimensionless units - * - */ -double WaterPropsIAPWS::helmholtzFE_RT() const{ - double retn = m_phi->phi(tau, delta); - return (retn); + /* + * Determine the internal state + */ + if (temperature > T_c) { + iState = WATER_SUPERCRIT; + } else { + if (delta < 1.0) { + iState = WATER_GAS; + } else { + iState = WATER_LIQUID; + } + } } /* @@ -84,18 +89,17 @@ double WaterPropsIAPWS::helmholtzFE_RT() const{ */ double WaterPropsIAPWS::helmholtzFE(double temperature, double rho) { setState(temperature, rho); - double retn = helmholtzFE_RT(); + double retn = m_phi->phi(tau, delta); double RT = Rgas * temperature; return (retn * RT); } double WaterPropsIAPWS::helmholtzFE() const{ - double retn = helmholtzFE_RT(); + double retn = m_phi->phi(tau, delta); double temperature = T_c/tau; double RT = Rgas * temperature; return (retn * RT); } - /* * Calculate the pressure (Pascals), given the temperature and density * Temperature: kelvin @@ -103,25 +107,16 @@ double WaterPropsIAPWS::helmholtzFE() const{ */ double WaterPropsIAPWS::pressure(double temperature, double rho) { calcDim(temperature, rho); - double retn = pressureM_rhoRT(); + double retn = m_phi->pressureM_rhoRT(tau, delta); return (retn * rho * Rgas * temperature/M_water); } double WaterPropsIAPWS::pressure() const{ - double retn = pressureM_rhoRT(); + double retn = m_phi->pressureM_rhoRT(tau, delta); double rho = delta * Rho_c; double temperature = T_c / tau; return (retn * rho * Rgas * temperature/M_water); } -/* - * Calculates the pressure in dimensionless form - * pM/(rhoRT) at the currently stored tau and delta values - */ -double WaterPropsIAPWS::pressureM_rhoRT() const { - double retn = m_phi->pressureM_rhoRT(tau, delta); - return retn; -} - /* * Calculates the density given the temperature and the pressure, * and a guess at the density. Note, below T_c, this is a @@ -150,7 +145,8 @@ density(double temperature, double pressure, int phase, double rhoguess) { rhoguess = pressure * M_water / (Rgas * temperature); } else { /* - * Provide a guess about the liquid density + * Provide a guess about the liquid density that is + * relatively high -> convergnce from above seems robust. */ rhoguess = 1000.; } @@ -166,7 +162,7 @@ density(double temperature, double pressure, int phase, double rhoguess) { } double p_red = pressure * M_water / (Rgas * temperature * Rho_c); deltaGuess = rhoguess / Rho_c; - calcDim(temperature, rhoguess); + setState(temperature, rhoguess); double delta_retn = m_phi->dfind(p_red, tau, deltaGuess); double density_retn; if (delta_retn >0.0) { @@ -177,17 +173,11 @@ density(double temperature, double pressure, int phase, double rhoguess) { */ density_retn = delta_retn * Rho_c; /* - * Determine the internal state + * Set the internal state -> this may be + * a duplication. However, let's just be sure. */ - if (temperature > T_c) { - iState = WATER_SUPERCRIT; - } else { - if (delta_retn < 1.0) { - iState = WATER_GAS; - } else { - iState = WATER_LIQUID; - } - } + setState(temperature, density_retn); + } else { density_retn = -1.0; @@ -327,34 +317,25 @@ isothermalCompressibility(double temperature, double pressure) { return retn; } -/** - * Calculate the Gibbs Free energy in dimensionless units - * - */ -double WaterPropsIAPWS:: -Gibbs_RT() const{ - double gRT = m_phi->gibbs_RT(); - return gRT; -} -/** +/* * Calculate the Gibbs free energy in mks units of * J kmol-1 K-1. */ double WaterPropsIAPWS:: Gibbs(double temperature, double rho) { setState(temperature, rho); - double gRT = Gibbs_RT(); + double gRT = m_phi->gibbs_RT(); return (gRT * Rgas * temperature); } double WaterPropsIAPWS:: Gibbs() const { - double gRT = Gibbs_RT(); + double gRT = m_phi->gibbs_RT(); double temperature = T_c/tau; return (gRT * Rgas * temperature); } -/** +/* * Calculate the Gibbs free energy in mks units of * J kmol-1 K-1. */ @@ -368,7 +349,7 @@ corr(double temperature, double pressure, double &densLiq, exit(-1); } setState(temperature, densLiq); - double gibbsLiqRT = Gibbs_RT(); + double gibbsLiqRT = m_phi->gibbs_RT(); densGas = density(temperature, pressure, WATER_GAS, densGas); if (densGas <= 0.0) { @@ -376,7 +357,7 @@ corr(double temperature, double pressure, double &densLiq, exit(-1); } setState(temperature, densGas); - double gibbsGasRT = Gibbs_RT(); + double gibbsGasRT = m_phi->gibbs_RT(); delGRT = gibbsLiqRT - gibbsGasRT; } @@ -443,133 +424,91 @@ setState(double temperature, double rho) { m_phi->tdpolycalc(tau, delta); } -/** - * Calculate the enthalpy in dimensionless units - * - */ -double WaterPropsIAPWS:: -enthalpy_RT() const{ - double hRT = m_phi->enthalpy_RT(); - return hRT; -} -/** + +/* * Calculate the enthalpy in mks units of * J kmol-1 K-1. */ double WaterPropsIAPWS:: enthalpy(double temperature, double rho) { setState(temperature, rho); - double hRT = enthalpy_RT(); + double hRT = m_phi->enthalpy_RT(); return (hRT * Rgas * temperature); } double WaterPropsIAPWS:: enthalpy() const { double temperature = T_c/tau; - double hRT = enthalpy_RT(); + double hRT = m_phi->enthalpy_RT(); return (hRT * Rgas * temperature); } -/** - * Calculate the internal Energy in dimensionless units - * - */ -double WaterPropsIAPWS:: -intEnergy_RT() const { - double uRT = m_phi->intEnergy_RT(); - return uRT; -} -/** +/* * Calculate the internal Energy in mks units of * J kmol-1 K-1. */ double WaterPropsIAPWS:: intEnergy(double temperature, double rho) { setState(temperature, rho); - double uRT = intEnergy_RT(); + double uRT = m_phi->intEnergy_RT(); return (uRT * Rgas * temperature); } double WaterPropsIAPWS:: intEnergy() const{ double temperature = T_c / tau; - double uRT = intEnergy_RT(); + double uRT = m_phi->intEnergy_RT(); return (uRT * Rgas * temperature); } -/** - * Calculate the enthalpy in dimensionless units - * - */ -double WaterPropsIAPWS:: -entropy_R() const { - double sR = m_phi->entropy_R(); - return sR; -} - -/** +/* * Calculate the enthalpy in mks units of * J kmol-1 K-1. */ double WaterPropsIAPWS:: entropy(double temperature, double rho) { setState(temperature, rho); - double sR = entropy_R(); + double sR = m_phi->entropy_R(); return (sR * Rgas); } -/** +/* * Calculate the enthalpy in mks units of * J kmol-1 K-1. */ double WaterPropsIAPWS:: entropy() const { - double sR = entropy_R(); + double sR = m_phi->entropy_R(); return (sR * Rgas); } -/** - * Calculate the dimensionless Heat capacity at constant volume - */ -double WaterPropsIAPWS:: -cv_R() const { - double cvR = m_phi->cv_R(); - return cvR; -} -/** +/* * Calculate heat capacity at constant volume * J kmol-1 K-1. */ double WaterPropsIAPWS:: cv(double temperature, double rho) { setState(temperature, rho); - double cvR = cv_R(); + double cvR = m_phi->cv_R(); return (cvR * Rgas); } -/** - * Calculate the dimensionless Heat capacity at constant pressure - */ -double WaterPropsIAPWS:: -cp_R() const { - double cpR = m_phi->cp_R(); - return cpR; -} -/** +/* * Calculate heat capacity at constant pressure * J kmol-1 K-1. */ double WaterPropsIAPWS:: cp(double temperature, double rho) { setState(temperature, rho); - double cpR = cp_R(); + double cpR = m_phi->cp_R(); return (cpR * Rgas); } + double WaterPropsIAPWS:: cp() const { - double cpR = cp_R(); + double cpR = m_phi->cp_R(); return (cpR * Rgas); } diff --git a/Cantera/src/thermo/WaterPropsIAPWS.h b/Cantera/src/thermo/WaterPropsIAPWS.h index 114a9d0b9..c508cb0fc 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.h +++ b/Cantera/src/thermo/WaterPropsIAPWS.h @@ -1,6 +1,6 @@ /** * @file WaterPropsIAPWS.h - * + * Definitions for a class for calculating the equation of state of water. */ /* * Copywrite (2005) Sandia Corporation. Under the terms of @@ -28,9 +28,41 @@ #define WATER_SUPERCRIT 2 //@} +//! Class for calculating the equation of state of water. /*! - * Class for calculating the properties of water. * + * The reference is W. Wagner, A. Prub, "The IAPWS Formulation 1995 for the Themodynamic + * Properties of Ordinary Water Substance for General and Scientific Use," + * J. Phys. Chem. Ref. Dat, 31, 387, 2002. + * + * This class provides a very complicated polynomial for the specific helmholtz free + * energy of water, as a function of temperature and density. + * + * \f[ + * \frac{M\hat{f}(\rho,T)}{R T} = \phi(\delta, \tau) = + * \phi^o(\delta, \tau) + \phi^r(\delta, \tau) + * \f] + * + * where + * + * \f[ + * \delta = \rho / \rho_c \mbox{\qquad and \qquad} \tau = T_c / T + * \f] + * + * The following constants are assumed + * + * \f[ + * T_c = 647.096\mbox{\ K} + * \f] + * \f[ + * \rho_c = 322 \mbox{\ kg\ m$^{-3}$} + * \f] + * \f[ + * R/M = 0.46151805 \mbox{\ kJ\ kg$^{-1}$\ K$^{-1}$} + * \f] + * + * The free energy is a unique single-valued function of the temperature and density + * over its entire range. * * Note, the base thermodynamic state for this class is the one * used in the steam tables, i.e., the liquid at the triple point @@ -41,6 +73,59 @@ * - psat(273.16) = 611.655 Pascal * - rho(273.16, psat) = 999.793 kg m-3 * + * Therefore, to use this class within %Cantera, offsets to u() and s() must be used + * to put the water class onto the same basis as other thermodynamic quantities. + * For example, in the WaterSSTP class, these offsets are calculated in the following way. + * The thermodynamic base state for water is set to the NIST basis here + * by specifying constants EW_Offset and SW_Offset. These offsets are + * calculated on the fly so that the following properties hold: + * + * - Delta_Hfo_idealGas(298.15, 1bar) = -241.826 kJ/gmol + * - So_idealGas(298.15, 1bar) = 188.835 J/gmolK + * + * The offsets are calculated by actually computing the above quantities and then + * calculating the correction factor. + * + * This class provides an interface to the #WaterPropsIAPWSphi class, which actually + * calculates the \f$ \phi^o(\delta, \tau) \f$ and the \f$ \phi^r(\delta, \tau) \f$ + * polynomials in dimensionless form. + * + * All thermodynamic results from this class are returned in dimensional form. This + * is because the gas constant (and molecular weight) used within this class is allowed to be potentially + * different than that used elsewhere in %Cantera. Therefore, everything has to be + * in dimensional units. Note, however, the thermodynamic basis is set to that used + * in the steam tables. (u = s = 0 for liquid water at the triple point). + * + * This class is not a %ThermoPhase. However, it does maintain an internal state of + * the object that is dependent on temperature and density. The internal state + * is characterized by an internally storred \f$ \tau\f$ and a \f$ \delta \f$ value, + * and an iState value, which indicates whether the point is a liquid, a gas, + * or a supercritical fluid. + * Along with that the \f$ \tau\f$ and a \f$ \delta \f$ values are polynomials of + * \f$ \tau\f$ and a \f$ \delta \f$ that are kept by the #WaterPropsIAPWSphi class. + * Therefore, whenever \f$ \tau\f$ or \f$ \delta \f$ is changed, the function setState() + * must be called in order for the internal state to be kept up to date. + * + * The class is pretty straightfoward. However, one function deserves mention. + * the #density() function calculates the density that is consistent with + * a particular value of the temperature and pressure. It may therefore be + * multivalued or potentially there may be no answer from this function. It therefore + * takes a phase guess and a density guess as optional parameters. If no guesses are + * supplied to density(), a gas phase guess is assumed. This may or may not be what + * is wanted. Therefore, density() should usually at leat be supplied with a phase + * guess so that it may manufacture an appropriate density guess. + * #density() manufactures the initial density guess, nondimensionalizes everything, + * and then calls #WaterPropsIAPWSphi::dfind(), which does the iterative calculation + * to find the density condition that matches the desired input pressure. + * + * The phase guess defines are located in the .h file. they are + * + * - WATER_GAS + * - WATER_LIQUID + * - WATER_SUPERCRIT + * + * @ingroup thermoprops + * */ class WaterPropsIAPWS { public: @@ -174,10 +259,20 @@ public: double pressure() const; //! Calculates the density given the temperature and the pressure, - //! and a guess at the density. + //! and a guess at the density. Sets the internal state. /*! * Note, below T_c, this is a multivalued function. * + * The #density() function calculates the density that is consistent with + * a particular value of the temperature and pressure. It may therefore be + * multivalued or potentially there may be no answer from this function. It therefore + * takes a phase guess and a density guess as optional parameters. If no guesses are + * supplied to density(), a gas phase guess is assumed. This may or may not be what + * is wanted. Therefore, density() should usually at leat be supplied with a phase + * guess so that it may manufacture an appropriate density guess. + * #density() manufactures the initial density guess, nondimensionalizes everything, + * and then calls #WaterPropsIAPWSphi::dfind(), which does the iterative calculation + * to find the density condition that matches the desired input pressure. * * @param temperature: Kelvin * @param pressure : Pressure in Pascals (Newton/m**2) @@ -186,7 +281,8 @@ public: * @param rhoguess : guessed density of the water * : -1.0 no guessed density * @return - * Returns the density + * Returns the density. If an error is encountered in the calculation + * the value of -1.0 is returned. */ double density(double temperature, double pressure, int phase = -1, double rhoguess = -1.0); @@ -197,19 +293,6 @@ public: */ double density() const; - - //! This function returns an estimated value for the saturation pressure. - /*! - * It does this via a polynomial fit of the vapor pressure curve. - * units = (Pascals) - * - * @param temperature Input temperature (Kelvin) - * - * @return - * Returns the estimated saturation pressure - */ - double psat_est(double temperature); - //! Returns the coefficient of thermal expansion as a function of temperature and pressure. /*! * alpha = d (ln V) / dT at constant P. @@ -235,31 +318,19 @@ public: * returns the isothermal compressibility */ double isothermalCompressibility(double temperature, double pressure); - - - //! Utility routine in the calculation of the saturation pressure + + //! This function returns an estimated value for the saturation pressure. /*! - * @param temperature temperature (kelvin) - * @param pressure pressure (Pascal) - * @param densLiq Output density of liquid - * @param densGas output Density of gas - * @param delGRT output delGRT - */ - void corr(double temperature, double pressure, double &densLiq, - double &densGas, double &delGRT); + * It does this via a polynomial fit of the vapor pressure curve. + * units = (Pascals) + * + * @param temperature Input temperature (Kelvin) + * + * @return + * Returns the estimated saturation pressure + */ + double psat_est(double temperature); - //! Utility routine in the calculation of the saturation pressure - /*! - * @param temperature temperature (kelvin) - * @param pressure pressure (Pascal) - * @param densLiq Output density of liquid - * @param densGas output Density of gas - * @param pcorr output corrected pressure - */ - void corr1(double temperature, double pressure, double &densLiq, - double &densGas, double &pcorr); - - //! This function returns the saturation pressure given the //! temperature as an input parameter. /*! @@ -297,35 +368,29 @@ private: */ void calcDim(double temperature, double rho); - /* - * Dimensionless versions of thermo functions. Note these are - * private, because R value is specific to the class. We only - * show the dimensional functions in the interface. + //! Utility routine in the calculation of the saturation pressure + /*! + * @param temperature temperature (kelvin) + * @param pressure pressure (Pascal) + * @param densLiq Output density of liquid + * @param densGas output Density of gas + * @param delGRT output delGRT */ - double helmholtzFE_RT() const; + void corr(double temperature, double pressure, double &densLiq, + double &densGas, double &delGRT); - //! Returns the dimensionless gibbs free energy - double Gibbs_RT() const; + //! Utility routine in the calculation of the saturation pressure + /*! + * @param temperature temperature (kelvin) + * @param pressure pressure (Pascal) + * @param densLiq Output density of liquid + * @param densGas output Density of gas + * @param pcorr output corrected pressure + */ + void corr1(double temperature, double pressure, double &densLiq, + double &densGas, double &pcorr); - //! Returns the dimensionless enthalpy - double enthalpy_RT() const; - - //! Returns the dimensionless internal energy - double intEnergy_RT() const; - - //! Returns the dimensionless entropy - double entropy_R() const; - - //! Returns the dimensionless heat capacity at constant volume - double cv_R() const; - - //! Returns the dimensionless heat capacity at constant pressure - double cp_R() const; - - //! Return the current dimensionless pressure - double pressureM_rhoRT() const; - -protected: +private: //! pointer to the underlying object that does the calculations. WaterPropsIAPWSphi *m_phi; @@ -346,4 +411,3 @@ protected: int iState; }; #endif - diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.h b/Cantera/src/thermo/WaterPropsIAPWSphi.h index a6f65d563..0022afd96 100644 --- a/Cantera/src/thermo/WaterPropsIAPWSphi.h +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.h @@ -1,7 +1,7 @@ /** * @file WaterPropsIAPWSphi.h - * * Lowest level of the classes which support a real water model. + * This class calculates dimensionless quantitites. */ /* * Copywrite (2006) Sandia Corporation. Under the terms of @@ -18,6 +18,10 @@ /*! * the WaterPropsIAPSWSphi class support low level calls for * the real description of water. + * + * The reference is W. Wagner, A. Prub, "The IAPWS Formulation 1995 for the Themodynamic + * Properties of Ordinary Water Substance for General and Scientific Use," + * J. Phys. Chem. Ref. Dat, 31, 387, 2002. * * Units Note: This class works with reduced units exclusively. */ @@ -80,7 +84,6 @@ public: //! Internal check # 2 void check2(); - //! Calculate the dimensionless pressure at tau and delta; /*! * @@ -138,9 +141,10 @@ public: */ double cp_R() const; - /** - * Calculates internal polynomials in tau and delta. This - * routine is used to store the internal state of tau and delta + + //! Calculates internal polynomials in tau and delta. + /*! + * This routine is used to store the internal state of tau and delta * for later use by the other routines in the class. * * @param tau Dimensionless temperature = T_c/T @@ -186,7 +190,7 @@ private: */ void intCheck(double tau, double delta); -protected: +private: //! Value of internally calculated polynomials of powers of TAU double TAUp[52]; diff --git a/Cantera/src/thermo/WaterSSTP.cpp b/Cantera/src/thermo/WaterSSTP.cpp index d50037383..2d6ee9d0e 100644 --- a/Cantera/src/thermo/WaterSSTP.cpp +++ b/Cantera/src/thermo/WaterSSTP.cpp @@ -1,7 +1,7 @@ /** * @file WaterSSTP.cpp - * Definitions for the Object WaterSSTP, which creates a - * single species ThermoPhase object for real liquid water. + * Declarations for the object WaterSSTP, which creates a + * single species %ThermoPhase object for real liquid water. */ /* * Copywrite (2006) Sandia Corporation. Under the terms of @@ -25,25 +25,25 @@ namespace Cantera { WaterSSTP::WaterSSTP() : SingleSpeciesTP(), m_sub(0), - m_subflag(0), m_mw(0.0), EW_Offset(0.0), SW_Offset(0.0), m_verbose(0), + m_ready(false), m_allowGasPhase(false) { - constructPhase(); + //constructPhase(); } WaterSSTP::WaterSSTP(std::string inputFile, std::string id) : SingleSpeciesTP(), m_sub(0), - m_subflag(0), m_mw(0.0), EW_Offset(0.0), SW_Offset(0.0), m_verbose(0), + m_ready(false), m_allowGasPhase(false) { constructPhaseFile(inputFile, id); @@ -53,11 +53,11 @@ namespace Cantera { WaterSSTP::WaterSSTP(XML_Node& phaseRoot, std::string id) : SingleSpeciesTP(), m_sub(0), - m_subflag(0), m_mw(0.0), EW_Offset(0.0), SW_Offset(0.0), m_verbose(0), + m_ready(false), m_allowGasPhase(false) { constructPhaseXML(phaseRoot, id) ; @@ -68,14 +68,14 @@ namespace Cantera { WaterSSTP::WaterSSTP(const WaterSSTP &b) : SingleSpeciesTP(b), m_sub(0), - m_subflag(b.m_subflag), m_mw(b.m_mw), EW_Offset(b.EW_Offset), SW_Offset(b.SW_Offset), m_verbose(b.m_verbose), + m_ready(false), m_allowGasPhase(b.m_allowGasPhase) { - m_sub = new WaterPropsIAPWS(*(b.m_sub)); + m_sub = new WaterPropsIAPWS(*(b.m_sub)); /* * Use the assignment operator to do the brunt * of the work for the copy construtor. @@ -89,9 +89,9 @@ namespace Cantera { WaterSSTP& WaterSSTP::operator=(const WaterSSTP&b) { if (&b == this) return *this; m_sub->operator=(*(b.m_sub)); - m_subflag = b.m_subflag; m_mw = b.m_mw; m_verbose = b.m_verbose; + m_ready = b.m_ready; m_allowGasPhase = b.m_allowGasPhase; return *this; } @@ -109,7 +109,7 @@ namespace Cantera { void WaterSSTP::constructPhase() { - throw CanteraError("constructPhaseXML", "unimplemented"); + throw CanteraError("WaterSSTP::constructPhase()", "unimplemented"); } @@ -183,10 +183,16 @@ namespace Cantera { void WaterSSTP::initThermo() { + SingleSpeciesTP::initThermo(); } void WaterSSTP:: initThermoXML(XML_Node& phaseNode, std::string id) { + + /* + * Do initializations that don't depend on knowing the XML file + */ + initThermo(); if (m_sub) delete m_sub; m_sub = new WaterPropsIAPWS(); if (m_sub == 0) { @@ -238,7 +244,7 @@ namespace Cantera { } s = entropy_mole(); s -= GasConstant * log(oneBar/presLow); - printf("s = %g\n", s); + //printf("s = %g\n", s); doublereal h = enthalpy_mole(); if (h != -241.826E6) { @@ -246,7 +252,7 @@ namespace Cantera { } h = enthalpy_mole(); - printf("h = %g\n", h); + //printf("h = %g\n", h); /* @@ -264,15 +270,16 @@ namespace Cantera { delete m_spthermo; m_spthermo = 0; } + + /* + * Set the flag to say we are ready to calculate stuff + */ + m_ready = true; } void WaterSSTP:: setParametersFromXML(const XML_Node& eosdata) { - eosdata._require("model","PureFluid"); - m_subflag = atoi(eosdata["fluid_type"].c_str()); - if (m_subflag < 0) - throw CanteraError("WaterSSTP::setParametersFromXML", - "missing or negative substance flag"); + eosdata._require("model","PureLiquidWater"); } /* @@ -315,6 +322,9 @@ namespace Cantera { double dens = density(); doublereal g = m_sub->Gibbs(T, dens); *grt = (g + EW_Offset - SW_Offset*T) / (GasConstant * T); + if (!m_ready) { + throw CanteraError("waterSSTP::", "Phase not ready"); + } } /* @@ -326,6 +336,9 @@ namespace Cantera { double dens = density(); doublereal g = m_sub->Gibbs(T, dens); *gss = (g + EW_Offset - SW_Offset*T); + if (!m_ready) { + throw CanteraError("waterSSTP::", "Phase not ready"); + } } void WaterSSTP::getCp_R(doublereal* cpr) const { diff --git a/Cantera/src/thermo/WaterSSTP.h b/Cantera/src/thermo/WaterSSTP.h index 49beb3dfb..d4567e2fa 100644 --- a/Cantera/src/thermo/WaterSSTP.h +++ b/Cantera/src/thermo/WaterSSTP.h @@ -24,19 +24,37 @@ namespace Cantera { //! Class for single-component water. This is designed to cover just the //! liquid part of water. /*! + * The reference is W. Wagner, A. Prub, "The IAPWS Formulation 1995 for the Themodynamic + * Properties of Ordinary Water Substance for General and Scientific Use," + * J. Phys. Chem. Ref. Dat, 31, 387, 2002. * + *
+ *

Specification of Species Standard %State Properties

+ *
+ * + * The offsets used in the steam tables are different than NIST's. + * They assume u_liq(TP) = 0.0, s_liq(TP) = 0.0, where TP is the + * triple point conditions: + * + * - u(273.16, rho) = 0.0 + * - s(273.16, rho) = 0.0 + * - psat(273.16) = 611.655 Pascal + * - rho(273.16, psat) = 999.793 kg m-3 + * + * These "steam table" assumptions are used by the WaterPropsIAPWS class. + * Therefore, offsets must be calculated to make the thermodynamic + * properties calculated within this class to be consistent with + * thermo properties within Cantera. * - * Notes: - * Base state for thermodynamic properties: - * * The thermodynamic base state for water is set to the NIST basis here - * by specifying constants EW_Offset and SW_Offset. These offsets are + * by specifying constants, #EW_Offset and #SW_Offset, one for energy + * quantities and one for entropy quantities. The offsets are * specified so that the following properties hold: * - * Delta_Hfo_gas(298.15) = -241.826 kJ/gmol - * So_gas(298.15, 1bar) = 188.835 J/gmolK + * - Delta_Hfo_idealgas(298.15) = -241.826 kJ/gmol + * - So_idealgas(298.15, 1bar) = 188.835 J/gmolK * - * (http://webbook.nist.gov) + * ref -> (http://webbook.nist.gov) * * The "o" here refers to a hypothetical ideal gas state. The way * we achieve this in practice is to evaluate at a very low pressure @@ -47,10 +65,69 @@ namespace Cantera { * * So(1bar) = S(P0) + RT ln(1bar/P0) * - * The offsets used in the steam tables are different than NIST's. - * They assume u_liq(TP) = 0.0, s_liq(TP) = 0.0, where TP is the - * triple point conditions. + *
+ *

%Application within %Kinetics Managers

+ *
* + * This is unimplemented. + * + *
+ *

Instantiation of the Class

+ *
+ * + * The constructor for this phase is NOT located in the default ThermoFactory + * for %Cantera. However, a new %WaterSSTP object may be created by + * the following code snippets, combined with an XML file given in the + * XML example section. + * + * @code + * WaterSSTP *w = new WaterSSTP("waterSSTPphase.xml",""); + * @endcode + * + * or + * + * @code + * char iFile[80], file_ID[80]; + * strcpy(iFile, "waterSSTPphase.xml"); + * sprintf(file_ID,"%s#water", iFile); + * XML_Node *xm = get_XML_NameID("phase", file_ID, 0); + * WaterSSTP *w = new WaterSSTP(*xm); + * @endcode + * + * or by the following call to importPhase(): + * + * @code + * char iFile[80], file_ID[80]; + * strcpy(iFile, "waterSSTPphase.xml"); + * sprintf(file_ID,"%s#water", iFile); + * XML_Node *xm = get_XML_NameID("phase", file_ID, 0); + * WaterSSTP water; + * importPhase(*xm, &water); + * @endcode + * + *
+ *

XML Example

+ *
+ * + * An example of an XML Element named phase setting up a WaterSSTP object with + * id "water" is given below. + * + * @verbatim + + + O H + H2O + + 300.0 + 101325.0 + + + + + @endverbatim + * + * Note the model "PureLiquidWater" indicates the usage of the WaterSSTP object. + * * @ingroup thermoprops * */ @@ -406,14 +483,8 @@ namespace Cantera { */ virtual void setParametersFromXML(const XML_Node& eosdata); - protected: - - void Set(int n, double x, double y) const; - void setTPXState() const; - void check(doublereal v = 0.0) const; - void reportTPXError() const; - protected: + /** * @internal * This internal routine must be overwritten because @@ -422,8 +493,11 @@ namespace Cantera { void _updateThermo() const; private: + //! Pointer to the WaterPropsIAPWS that calculates the real properties + //! of water. mutable WaterPropsIAPWS *m_sub; - int m_subflag; + + //! Molecular weight of Water -> Cantera assumption doublereal m_mw; /** @@ -442,6 +516,8 @@ namespace Cantera { bool m_verbose; + bool m_ready; + /** * Since this phase represents a liquid phase, it's an error to * return a gas-phase answer. However, if the below is true, then diff --git a/test_problems/cathermo/testWaterTP/output_blessed.txt b/test_problems/cathermo/testWaterTP/output_blessed.txt index 507913150..92b8ff263 100644 --- a/test_problems/cathermo/testWaterTP/output_blessed.txt +++ b/test_problems/cathermo/testWaterTP/output_blessed.txt @@ -1,5 +1,3 @@ -s = 188835 -h = -2.41826e+08 psat(273.16) = 611.655 Comparisons to NIST: (see http://webbook.nist.gov): diff --git a/test_problems/cathermo/testWaterTP/testWaterSSTP.cpp b/test_problems/cathermo/testWaterTP/testWaterSSTP.cpp index cec3e0fc6..42a2d3ba5 100644 --- a/test_problems/cathermo/testWaterTP/testWaterSSTP.cpp +++ b/test_problems/cathermo/testWaterTP/testWaterSSTP.cpp @@ -4,6 +4,7 @@ #include "stdio.h" #include "math.h" #include "WaterSSTP.h" +#include "importCTML.h" #include using namespace std; using namespace Cantera; @@ -23,7 +24,21 @@ int main () { double pres; try { WaterSSTP *w = new WaterSSTP("waterTPphase.xml",""); + delete w; + char iFile[80], file_ID[80]; + strcpy(iFile, "waterTPphase.xml"); + sprintf(file_ID,"%s#water", iFile); + XML_Node *xm = get_XML_NameID("phase", file_ID, 0); + w = new WaterSSTP(*xm); + delete w; + + strcpy(iFile, "waterTPphase.xml"); + sprintf(file_ID,"%s#water", iFile); + xm = get_XML_NameID("phase", file_ID, 0); + w = new WaterSSTP(); + importPhase(*xm, w); + /* * Print out the triple point conditions diff --git a/test_problems/cathermo/testWaterTP/waterTPphase.xml b/test_problems/cathermo/testWaterTP/waterTPphase.xml index 827bc2948..a58b5a2ba 100644 --- a/test_problems/cathermo/testWaterTP/waterTPphase.xml +++ b/test_problems/cathermo/testWaterTP/waterTPphase.xml @@ -10,7 +10,7 @@ 300.0 101325.0 - +