From ebdc0d014ba3e74018aca4c9aea6995d93d9dc9b Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Mon, 12 Mar 2007 01:12:57 +0000 Subject: [PATCH] Fixed some errors in DebyeHuckel that weren't covered by test suite (partial molar enthalpy and heat capacity) Doxygen update - added private functions and added some water property files. --- Cantera/src/ConstDensityThermo.h | 2 +- Cantera/src/Kinetics.cpp | 26 +- Cantera/src/Kinetics.h | 10 +- Cantera/src/NasaThermo.h | 14 ++ Cantera/src/Phase.h | 91 ++++--- Cantera/src/ThermoPhase.h | 11 +- Cantera/src/thermo/DebyeHuckel.cpp | 39 +-- Cantera/src/thermo/DebyeHuckel.h | 41 ++-- Cantera/src/thermo/SingleSpeciesTP.h | 6 +- Cantera/src/thermo/VPStandardStateTP.h | 6 +- Cantera/src/thermo/WaterPropsIAPWS.cpp | 10 +- Cantera/src/thermo/WaterPropsIAPWS.h | 284 ++++++++++++++++------ Cantera/src/thermo/WaterPropsIAPWSphi.cpp | 18 +- Cantera/src/thermo/WaterPropsIAPWSphi.h | 111 ++++++++- Cantera/src/units.h | 20 +- tools/doc/Cantera.cfg.in | 5 +- 16 files changed, 506 insertions(+), 188 deletions(-) diff --git a/Cantera/src/ConstDensityThermo.h b/Cantera/src/ConstDensityThermo.h index f1bcd4673..51df0d8ac 100755 --- a/Cantera/src/ConstDensityThermo.h +++ b/Cantera/src/ConstDensityThermo.h @@ -443,9 +443,9 @@ namespace Cantera { //! Current pressure (Pa) doublereal m_press; - private: + //! Function to update the reference state thermo functions void _updateThermo() const; }; } diff --git a/Cantera/src/Kinetics.cpp b/Cantera/src/Kinetics.cpp index ff0941865..a0b3f5e88 100644 --- a/Cantera/src/Kinetics.cpp +++ b/Cantera/src/Kinetics.cpp @@ -214,19 +214,17 @@ namespace Cantera { m_phaseindex[m_thermo.back()->id()] = nPhases(); } - /** - * err(): - * - * Private function of the class Kinetics, indicating that a function - * inherited from the base class hasn't had a definition assigned to it - */ - void Kinetics::err(std::string m) const { - throw CanteraError("Kinetics::" + m, - "The default Base class method was called, when " - "the inherited class's method should " - "have been called"); - } - - + + //! Private function of the class Kinetics, indicating that a function + //! inherited from the base class hasn't had a definition assigned to it + /*! + * @param m String message + */ + void Kinetics::err(std::string m) const { + throw CanteraError("Kinetics::" + m, + "The default Base class method was called, when " + "the inherited class's method should " + "have been called"); + } } diff --git a/Cantera/src/Kinetics.h b/Cantera/src/Kinetics.h index 092d40f76..676399011 100755 --- a/Cantera/src/Kinetics.h +++ b/Cantera/src/Kinetics.h @@ -936,8 +936,14 @@ namespace Cantera { private: - std::vector m_dummygroups; - void err(std::string m) const; + //! Vector of group lists + std::vector m_dummygroups; + + //! Function for unhandled situations + /*! + * @param m String error message + */ + void err(std::string m) const; }; diff --git a/Cantera/src/NasaThermo.h b/Cantera/src/NasaThermo.h index 093e6f8e8..c1d481f31 100755 --- a/Cantera/src/NasaThermo.h +++ b/Cantera/src/NasaThermo.h @@ -468,10 +468,20 @@ namespace Cantera { private: //! see SpeciesThermoFactory.cpp for the definition + /*! + * @param name string name of species + * @param tmid Mid temperature, between the two temperature regions + * @param clow coefficients for lower temperature region + * @param chigh coefficients for higher temperature region + */ void checkContinuity(std::string name, double tmid, const doublereal* clow, doublereal* chigh); //! for internal use by checkContinuity + /*! + * @param t temperature + * @param c coefficient array + */ doublereal enthalpy_RT(double t, const doublereal* c) { return c[0] + 0.5*c[1]*t + OneThird*c[2]*t*t + 0.25*c[3]*t*t*t + 0.2*c[4]*t*t*t*t @@ -479,6 +489,10 @@ namespace Cantera { } //! for internal use by checkContinuity + /*! + * @param t temperature + * @param c coefficient array + */ doublereal entropy_R(double t, const doublereal* c) { return c[0]*log(t) + c[1]*t + 0.5*c[2]*t*t + OneThird*c[3]*t*t*t + 0.25*c[4]*t*t*t*t diff --git a/Cantera/src/Phase.h b/Cantera/src/Phase.h index 0ea2cacf1..edb964386 100755 --- a/Cantera/src/Phase.h +++ b/Cantera/src/Phase.h @@ -71,10 +71,10 @@ namespace Cantera { * The XML_Node for the phase contains all of the input data used * to set up the model for the phase, during its initialization. */ - XML_Node& xml() { return *m_xml; } + XML_Node& xml() { return *m_xml; } //! Return the string id for the phase - std::string id() const { return m_id; } + std::string id() const { return m_id; } //! Set the string id for the phase /*! @@ -92,7 +92,7 @@ namespace Cantera { void setName(std::string nm) { m_name = nm; } //! Returns the index of the phase - int index() const { return m_index; } + int index() const { return m_index; } //! Sets the index of the phase /*! @@ -329,10 +329,10 @@ namespace Cantera { */ doublereal massFraction(std::string name) const; - /** - * Charge density [C/m^3]. - */ - doublereal chargeDensity() const; + /** + * Charge density [C/m^3]. + */ + doublereal chargeDensity() const; /// Returns the number of spatial dimensions (1, 2, or 3) int nDim() {return m_ndim;} @@ -346,47 +346,64 @@ namespace Cantera { */ void setNDim(int ndim) {m_ndim = ndim;} - /** - * Finished adding species, prepare to use them for calculation - * of mixture properties. - */ - virtual void freezeSpecies(); + /** + * Finished adding species, prepare to use them for calculation + * of mixture properties. + */ + virtual void freezeSpecies(); - virtual bool ready() const; + virtual bool ready() const; protected: - /** - * m_kk = Number of species in the phase. @internal m_kk is a - * member of both the State and Constituents classes. - * Therefore, to avoid multiple inheritance problems, we need - * to restate it in here, so that the declarations in the two - * base classes become hidden. - */ - int m_kk; - /** - * m_ndim is the dimensionality of the phase. Volumetric - * phases have dimensionality 3 and surface phases have - * dimensionality 2. - */ - int m_ndim; - /** - * m_index is the index of the phase - * - */ - int m_index; + /** + * m_kk = Number of species in the phase. @internal m_kk is a + * member of both the State and Constituents classes. + * Therefore, to avoid multiple inheritance problems, we need + * to restate it in here, so that the declarations in the two + * base classes become hidden. + */ + int m_kk; + /** + * m_ndim is the dimensionality of the phase. Volumetric + * phases have dimensionality 3 and surface phases have + * dimensionality 2. + */ + int m_ndim; + /** + * m_index is the index of the phase + * + */ + int m_index; private: - vector_fp m_data; - XML_Node* m_xml; - std::string m_id; - std::string m_name; + //! This stores the initial state of the system + /*! + * @deprecated + * This doesn't seem to be used much anymore. + */ + vector_fp m_data; + + //! Pointer to the XML node containing the XML info for this phase + XML_Node* m_xml; + + //! ID of the phase. + /*! + * This is the value of the ID attribute of the XML phase node. + */ + std::string m_id; + + //! Name of the phase. + /*! + * Initially, this is the value of the ID attribute of the XML phase node. + */ + std::string m_name; }; //! typedef for the base Phase class - typedef Phase phase_t; + typedef Phase phase_t; } #endif diff --git a/Cantera/src/ThermoPhase.h b/Cantera/src/ThermoPhase.h index ee4dff675..148d92247 100755 --- a/Cantera/src/ThermoPhase.h +++ b/Cantera/src/ThermoPhase.h @@ -1394,14 +1394,19 @@ namespace Cantera { private: - doublereal err(std::string msg) const; + //! Error function that gets called for unhandled cases + /*! + * @param msg String containing the message. + */ + doublereal err(std::string msg) const; }; //! typedef for the ThermoPhase class - typedef ThermoPhase thermophase_t; + typedef ThermoPhase thermophase_t; + //! typedef for the ThermoPhase class - typedef ThermoPhase thermo_t; + typedef ThermoPhase thermo_t; } #endif diff --git a/Cantera/src/thermo/DebyeHuckel.cpp b/Cantera/src/thermo/DebyeHuckel.cpp index fa85788cb..348cfe5d9 100644 --- a/Cantera/src/thermo/DebyeHuckel.cpp +++ b/Cantera/src/thermo/DebyeHuckel.cpp @@ -1987,7 +1987,7 @@ namespace Cantera { } - /** + /* * _activityWaterHelgesonFixedForm() * * Formula for the log of the activity of the water @@ -2300,6 +2300,7 @@ namespace Cantera { void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const { double z_k, coeff, tmp, y, yp1, sigma, tmpLn; int k; + // First we store dAdT explicitly here double dAdT = dA_DebyedT_TP(); if (dAdT == 0.0) { for (k = 0; k < m_kk; k++) { @@ -2318,13 +2319,17 @@ namespace Cantera { double sqrtI = sqrt(m_IionicMolality); double numdAdTTmp = dAdT * sqrtI; double denomTmp = m_B_Debye * sqrtI; + double d_lnActivitySolvent_dT = 0; switch (m_formDH) { case DHFORM_DILUTE_LIMIT: - for (int k = 0; k < m_kk; k++) { + for (int k = 1; k < m_kk; k++) { m_dlnActCoeffMolaldT[k] = m_lnActCoeffMolal[k] * dAdT / m_A_Debye; } + d_lnActivitySolvent_dT = 2.0 / 3.0 * dAdT * m_Mnaught * + m_IionicMolality * sqrt(m_IionicMolality); + m_dlnActCoeffMolaldT[m_indexSolvent] = d_lnActivitySolvent_dT; break; case DHFORM_BDOT_AK: @@ -2365,7 +2370,8 @@ namespace Cantera { sigma = 0.0; } m_dlnActCoeffMolaldT[m_indexSolvent] = - 2.0 /3.0 * dAdT * m_Mnaught * m_IionicMolality * sqrtI * sigma; + 2.0 /3.0 * dAdT * m_Mnaught * + m_IionicMolality * sqrtI * sigma; break; case DHFORM_BETAIJ: @@ -2384,8 +2390,7 @@ namespace Cantera { } else { sigma = 0.0; } - m_dlnActCoeffMolaldT[m_indexSolvent] = - (xmolSolvent - 1.0)/xmolSolvent + + m_dlnActCoeffMolaldT[m_indexSolvent] = 2.0 /3.0 * dAdT * m_Mnaught * m_IionicMolality * sqrtI * sigma; break; @@ -2405,8 +2410,7 @@ namespace Cantera { } sigma = 1.0 / ( 1.0 + denomTmp); - m_dlnActCoeffMolaldT[m_indexSolvent] = - (xmolSolvent - 1.0)/xmolSolvent + + m_dlnActCoeffMolaldT[m_indexSolvent] = 2.0 /3.0 * dAdT * m_Mnaught * m_IionicMolality * sqrtI * sigma; break; @@ -2416,6 +2420,8 @@ namespace Cantera { exit(-1); break; } + + } /* @@ -2521,8 +2527,7 @@ namespace Cantera { } else { sigma = 0.0; } - m_d2lnActCoeffMolaldT2[m_indexSolvent] = - (xmolSolvent - 1.0)/xmolSolvent + + m_d2lnActCoeffMolaldT2[m_indexSolvent] = 2.0 /3.0 * d2AdT2 * m_Mnaught * m_IionicMolality * sqrtI * sigma; break; @@ -2533,17 +2538,16 @@ namespace Cantera { for (int k = 0; k < m_kk; k++) { if (k != m_indexSolvent) { z_k = m_speciesCharge[k]; - m_dlnActCoeffMolaldT[k] = + m_d2lnActCoeffMolaldT2[k] = - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp) - 2.0 * z_k * z_k * d2AdT2 * tmpLn / (m_B_Debye * m_Aionic[0]); - m_dlnActCoeffMolaldT[k] /= 3.0; + m_d2lnActCoeffMolaldT2[k] /= 3.0; } } sigma = 1.0 / ( 1.0 + denomTmp); - m_dlnActCoeffMolaldT[m_indexSolvent] = - (xmolSolvent - 1.0)/xmolSolvent + + m_d2lnActCoeffMolaldT2[m_indexSolvent] = 2.0 /3.0 * d2AdT2 * m_Mnaught * m_IionicMolality * sqrtI * sigma; break; @@ -2641,7 +2645,8 @@ namespace Cantera { sigma = 0.0; } m_dlnActCoeffMolaldP[m_indexSolvent] = - 2.0 /3.0 * dAdP * m_Mnaught * m_IionicMolality * sqrtI * sigma; + 2.0 /3.0 * dAdP * m_Mnaught * + m_IionicMolality * sqrtI * sigma; break; case DHFORM_BETAIJ: @@ -2660,8 +2665,7 @@ namespace Cantera { } else { sigma = 0.0; } - m_dlnActCoeffMolaldP[m_indexSolvent] = - (xmolSolvent - 1.0)/xmolSolvent + + m_dlnActCoeffMolaldP[m_indexSolvent] = 2.0 /3.0 * dAdP * m_Mnaught * m_IionicMolality * sqrtI * sigma; break; @@ -2681,8 +2685,7 @@ namespace Cantera { } sigma = 1.0 / ( 1.0 + denomTmp); - m_dlnActCoeffMolaldP[m_indexSolvent] = - (xmolSolvent - 1.0)/xmolSolvent + + m_dlnActCoeffMolaldP[m_indexSolvent] = 2.0 /3.0 * dAdP * m_Mnaught * m_IionicMolality * sqrtI * sigma; break; diff --git a/Cantera/src/thermo/DebyeHuckel.h b/Cantera/src/thermo/DebyeHuckel.h index 0440137bc..d8783fc50 100644 --- a/Cantera/src/thermo/DebyeHuckel.h +++ b/Cantera/src/thermo/DebyeHuckel.h @@ -53,7 +53,9 @@ namespace Cantera { * The concentrations of the ionic species are assumed to obey the electroneutrality * condition. * - * Specification of Species Standard %State Properties + *
+ *

Specification of Species Standard %State Properties

+ *
* * The standard states are on the unit molality basis. Therefore, in the * documentation below, the normal \f$ o \f$ superscript is replaced with @@ -454,10 +456,10 @@ namespace Cantera { * Currently, \f$ B_{Debye} \f$ is a constant in the model, specified either by a default * water value, or through the input file. This may have to be looked at, in the future. * - * *
*

%Application within %Kinetics Managers

*
+ * * For the time being, we have set the standard concentration for all species in * this phase equal to the default concentration of the solvent at 298 K and 1 atm. * This means that the @@ -507,7 +509,7 @@ namespace Cantera { * Note, this treatment may be modified in the future, as events dictate. * *
- * Instantiation of the Class + *

Instantiation of the Class

*
* * The constructor for this phase is NOT located in the default ThermoFactory @@ -540,7 +542,7 @@ namespace Cantera { * @endcode * *
- * XML Example + *

XML Example

*
* * The phase model name for this is called StoichSubstance. It must be supplied @@ -597,8 +599,6 @@ namespace Cantera { @endverbatim * - * The model attribute, "StoichSubstanceSSTP", on the thermo element identifies the phase as - * being a StoichSubstanceSSTP object. * */ class DebyeHuckel : public MolalityVPSSTP { @@ -994,16 +994,19 @@ namespace Cantera { //! in the mixture. Units (J/kmol) /*! * For this phase, the partial molar enthalpies are equal to the - * pure species enthalpies + * standard state enthalpies modified by the derivative of the + * molality-based activity coefficent wrt temperature + * * \f[ - * \bar h_k(T,P) = \hat h^{ref}_k(T) + (P - P_{ref}) \hat V^0_k + * \bar h_k(T,P) = h^{\triangle}_k(T,P) - R T^2 \frac{d \ln(\gamma_k^\triangle)}{dT} * \f] - * The reference-state pure-species enthalpies, - * \f$ \hat h^{ref}_k(T) \f$, - * at the reference pressure,\f$ P_{ref} \f$, - * are computed by the species thermodynamic - * property manager. They are polynomial functions of temperature. - * @see SpeciesThermo + * The solvent partial molar enthalpy is equal to + * \f[ + * \bar h_o(T,P) = h^{o}_o(T,P) - R T^2 \frac{d \ln(a_o}{dT} + * \f] + * + * The temperature dependence of the activity coefficients currently + * only occurs through the temperature dependence of the Debye constant. * * @param hbar Output vector of species partial molar enthalpies. * Length: m_kk. units are J/kmol. @@ -1562,6 +1565,8 @@ namespace Cantera { //! salt-out modifications. /*! * Returns the calculated activity coefficients. + * + * @param IionicMolality Value of the ionic molality (sqrt(gmol/kg)) */ double _nonpolarActCoeff(double IionicMolality) const; @@ -1572,6 +1577,12 @@ namespace Cantera { * NaCl brine. It's to be used with extreme caution. */ double _osmoticCoeffHelgesonFixedForm() const; + + //! Formula for the log of the water activity that occurs in the GWB. + /*! + * It is originally from Helgeson for a variable + * NaCl brine. It's to be used with extreme caution. + */ double _lnactivityWaterHelgesonFixedForm() const; @@ -1868,6 +1879,8 @@ namespace Cantera { * * We assume that the activity coefficients are current in this routine * + * + * * The solvent activity coefficient is on the molality scale. It's derivative is too. */ void s_update_dlnMolalityActCoeff_dT() const; diff --git a/Cantera/src/thermo/SingleSpeciesTP.h b/Cantera/src/thermo/SingleSpeciesTP.h index 7df40963c..ba583222f 100644 --- a/Cantera/src/thermo/SingleSpeciesTP.h +++ b/Cantera/src/thermo/SingleSpeciesTP.h @@ -703,8 +703,12 @@ namespace Cantera { void _updateThermo() const; private: - doublereal err(std::string msg) const; + //! Error return for unhandled cases + /*! + * @param msg String message + */ + doublereal err(std::string msg) const; }; } diff --git a/Cantera/src/thermo/VPStandardStateTP.h b/Cantera/src/thermo/VPStandardStateTP.h index a40da3245..6dc1e6cea 100644 --- a/Cantera/src/thermo/VPStandardStateTP.h +++ b/Cantera/src/thermo/VPStandardStateTP.h @@ -572,12 +572,12 @@ namespace Cantera { */ mutable vector_fp m_Vss; - - + private: + //! VPStandardStateTP has its own err routine /*! - * VPStandardStateTP has its own err routine + * @param msg Error message string */ doublereal err(std::string msg) const; diff --git a/Cantera/src/thermo/WaterPropsIAPWS.cpp b/Cantera/src/thermo/WaterPropsIAPWS.cpp index 35c685eaf..92b34721e 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWS.cpp @@ -78,7 +78,7 @@ double WaterPropsIAPWS::helmholtzFE_RT() const{ return (retn); } -/** +/* * Calculate the Helmholtz free energy in mks units of * J kmol-1 K-1. */ @@ -96,7 +96,7 @@ double WaterPropsIAPWS::helmholtzFE() const{ } -/** +/* * Calculate the pressure (Pascals), given the temperature and density * Temperature: kelvin * rho: density in kg m-3 @@ -113,7 +113,7 @@ double WaterPropsIAPWS::pressure() const{ return (retn * rho * Rgas * temperature); } -/** +/* * Calculates the pressure in dimensionless form * p/(rhoRT) at the currently stored tau and delta values */ @@ -199,7 +199,7 @@ double WaterPropsIAPWS::density() const { return (delta * Rho_c); } -/** +/* * psat_est provides a rough estimate of the saturation * pressure given the temperature. This is used as an initial * guess for refining the pressure. @@ -245,7 +245,7 @@ double WaterPropsIAPWS::psat_est(double temperature) { return ps; } -/** +/* * Returns the coefficient of thermal expansion as a function * of temperature and pressure. * alpha = d (ln V) / dT at constant P. diff --git a/Cantera/src/thermo/WaterPropsIAPWS.h b/Cantera/src/thermo/WaterPropsIAPWS.h index fe63aff77..66ba81a36 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.h +++ b/Cantera/src/thermo/WaterPropsIAPWS.h @@ -16,15 +16,19 @@ #include "WaterPropsIAPWSphi.h" -/* - * These constants are defined and used in the interphase - * to describe desired phases. - */ +/** + * @name Names for the phase regions + * + * These constants are defined and used in the interface + * to describe desired phases. + */ +//@{ #define WATER_GAS 0 #define WATER_LIQUID 1 #define WATER_SUPERCRIT 2 +//@} -/** +/*! * Class for calculating the properties of water. * * @@ -32,149 +36,264 @@ * used in the steam tables, i.e., the liquid at the triple point * for water has the following properties: * - * 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 + * - 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 * */ class WaterPropsIAPWS { public: + + //! Base constructor WaterPropsIAPWS(); - WaterPropsIAPWS(const WaterPropsIAPWS &b); - WaterPropsIAPWS & operator=(const WaterPropsIAPWS &b); + + //! Copy constructor + WaterPropsIAPWS(const WaterPropsIAPWS &); + + //! assignment constructor + WaterPropsIAPWS & operator=(const WaterPropsIAPWS &); + + //! destructor ~WaterPropsIAPWS(); - + //! Set the internal state of the object wrt temperature and density + /*! + * @param temperature temperature (kelvin) + * @param rho density (kg m-3) + */ void setState(double temperature, double rho); - /** - * Calculate the Helmholtz free energy in mks units of - * J kmol-1 K-1. + + //! Calculate the Helmholtz free energy in mks units of J kmol-1 K-1. + /*! + * @param temperature temperature (kelvin) + * @param rho density (kg m-3) */ double helmholtzFE(double temperature, double rho); + + //! Calculate the Helmholtz free energy in mks units of J kmol-1 K-1, + //! using the last temperature and density double helmholtzFE() const; - /** - * Calculate the Gibbs free energy in mks units of - * J kmol-1 K-1. + + //! Calculate the Gibbs free energy in mks units of J kmol-1 K-1. + /*! + * @param temperature temperature (kelvin) + * @param rho density (kg m-3) */ double Gibbs(double temperature, double rho); + + //! Calculate the Gibbs free energy in mks units of J kmol-1 K-1. + //! using the last temperature and density double Gibbs() const; - /** - * Calculate the enthalpy in mks units of - * J kmol-1 + + //! Calculate the enthalpy in mks units of J kmol-1 + /*! + * @param temperature temperature (kelvin) + * @param rho density (kg m-3) */ double enthalpy(double temperature, double rho); + + //! Calculate the enthalpy in mks units of J kmol-1 + //! using the last temperature and density double enthalpy() const; - /** - * Calculate the internal energy in mks units of - * J kmol-1 + //! Calculate the internal energy in mks units of J kmol-1 + /*! + * @param temperature temperature (kelvin) + * @param rho density (kg m-3) */ double intEnergy(double temperature, double rho); - double intEnergy() const; - /** - * Calculate the entropy in mks units of - * J kmol-1 K-1 + //! Calculate the internal energy in mks units of J kmol-1 + //! at the last internal energy + double intEnergy() const; + + //! Calculate the entropy in mks units of J kmol-1 K-1 + /*! + * @param temperature temperature (kelvin) + * @param rho density (kg m-3) */ double entropy(double temperature, double rho); + + //! Calculate the entropy in mks units of J kmol-1 K-1 + //! at the last temperature and density double entropy() const; - /** - * Calculate the constant volume heat capacity - * in mks units of J kmol-1 K-1 + + //! Calculate the constant volume heat capacity in mks units of J kmol-1 K-1 + /*! + * @param temperature temperature (kelvin) + * @param rho density (kg m-3) */ double cv(double temperature, double rho); + + //! Calculate the constant volume heat capacity in mks units of J kmol-1 K-1 + //! at the last temperature and density double cv() const; - /** - * Calculate the constant pressure heat capacity - * in mks units of J kmol-1 K-1 + //! Calculate the constant pressure heat capacity in mks units of J kmol-1 K-1 + /*! + * @param temperature temperature (kelvin) + * @param rho density (kg m-3) */ double cp(double temperature, double rho); + + //! Calculate the constant pressure heat capacity in mks units of J kmol-1 K-1 + //! at the last temperature and density double cp() const; + //! Calculate the molar volume (kmol m-3) + /*! + * @param temperature temperature (kelvin) + * @param rho density (kg m-3) + */ double molarVolume(double temperature, double rho); - double molarVolume() const; - /** - * Calculate the pressure (Pascals), given the temperature and density - * Temperature: kelvin - * rho: density in kg m-3 + //! Calculate the molar volume (kmol m-3) + //! at the last temperature and density + double molarVolume() const; + + //! Calculate the pressure (Pascals), given the temperature and density + /*! + * @param temperature input temperature kelvin + * @param rho density in kg m-3 + * + * @return + * returns the pressure (Pascal) */ double pressure(double temperature, double rho); + + //! Calculates the pressure (Pascals), given the current value of the + //! temperature and density. + /*! + * The density is an independent variable in the underlying equation of state + */ double pressure() const; - /* - * Calculates the density given the temperature and the pressure, - * and a guess at the density. Note, below T_c, this is a - * multivalued function. + //! Calculates the density given the temperature and the pressure, + //! and a guess at the density. + /*! + * Note, below T_c, this is a multivalued function. * - * parameters: - * temperature: Kelvin - * pressure : Pressure in Pascals (Newton/m**2) - * phase : guessed phase of water - * : -1: no guessed phase - * rhoguess : guessed density of the water - * : -1.0 no guessed density + * + * @param temperature: Kelvin + * @param pressure : Pressure in Pascals (Newton/m**2) + * @param phase : guessed phase of water + * : -1: no guessed phase + * @param rhoguess : guessed density of the water + * : -1.0 no guessed density + * @return + * Returns the density */ double density(double temperature, double pressure, int phase = -1, double rhoguess = -1.0); + + //! Returns the density (kg m-3) + /*! + * The density is an independent variable in the underlying equation of state + */ 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. + + //! 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. + //! Returns the coefficient of thermal expansion as a function of temperature and pressure. + /*! * alpha = d (ln V) / dT at constant P. * - * + * @param temperature Input temperature (Kelvin) + * @param pressure Input pressure (Pa) + * @return + * Returns the coefficient of thermal expansion */ double coeffThermExp(double temperature, double pressure); - /** - * Returns the coefficient of isothermal compressibility as a function - * of temperature and pressure. + + //! Returns the coefficient of isothermal compressibility as a function + //! of temperature and pressure. + /*! * kappa = - d (ln V) / dP at constant T. * * units - 1/Pascal + * + * @param temperature Input temperature (Kelvin) + * @param pressure Input pressure (Pa) + * @return + * returns the isothermal compressibility */ double isothermalCompressibility(double temperature, double pressure); - /** - * Utility routine in the calculation of the saturation pressure - */ + + //! 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 + */ void corr(double temperature, double pressure, double &densLiq, double &densGas, double &delGRT); + + //! 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. - * units = Pascal + + //! This function returns the saturation pressure given the + //! temperature as an input parameter. + /*! + * @param temperature input temperature (kelvin) + * @return + * Returns the saturation pressure + * units = Pascal */ double psat(double temperature); + //! Returns the critical temperature of water (Kelvin) + /*! + * This is hard coded to the value 647.096 Kelvin + */ double Tcrit() { return 647.096;} + + //! Returns the critical pressure of water (22.064E6 Pa) + /*! + * This is hard coded to the value of 22.064E6 pascals + */ double Pcrit() { return 22.064E6;} + //! Return the critical density of water (kg m-3) + /*! + * This is equal to 322 kg m-3. + */ double Rhocrit() { return 322.;} - private: /** * Calculate the dimensionless temp and rho and store internally. + * + * @param temperature input temperature (kelvin) + * @param rho density in kg m-3 */ void calcDim(double temperature, double rho); @@ -184,19 +303,46 @@ private: * show the dimensional functions in the interface. */ double helmholtzFE_RT() const; + + //! Returns the dimensionless gibbs free energy double Gibbs_RT() const; + + //! 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 pressure_rhoRT() const; protected: + + //! pointer to the underlying object that does the calculations. WaterPropsIAPWSphi *m_phi; + //! Dimensionless temperature + /*! + * tau = T_C / T + */ double tau; + + //! Dimensionless density + /*! + * delta = rho / rho_c + */ double delta; + + //! Current state of the system int iState; }; #endif diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.cpp b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp index bc2d6e1f3..e0d764e3f 100644 --- a/Cantera/src/thermo/WaterPropsIAPWSphi.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp @@ -363,7 +363,7 @@ WaterPropsIAPWSphi::WaterPropsIAPWSphi() : { } -/** +/* * intCheck() calculates all of the functions at a one point and * prints out the result. It's used for conducting the internal * check. @@ -624,7 +624,7 @@ double WaterPropsIAPWSphi::phi_d(double tau, double delta) { return retn; } -/** +/* * Calculate the dimensionless pressure at tau and delta; * * p/(rhoRT) = delta * phi_d() @@ -854,7 +854,7 @@ double WaterPropsIAPWSphi::phiR_t() const { return val; } -/** +/* * Calculate the dPhidtau function, which is basically the derivative * of helmholtz free energy wrt tau * Eqn. (6.4) @@ -867,7 +867,7 @@ double WaterPropsIAPWSphi::phi_t(double tau, double delta) { return retn; } -/** +/* * Calculate d2_phi0/dtau2 */ double WaterPropsIAPWSphi::phi0_tt() const { @@ -882,7 +882,7 @@ double WaterPropsIAPWSphi::phi0_tt() const { return retn; } -/** +/* * Calculate Eqn. 6.6 for dphiRdtau, the second derivative residual part of the * dimensionless Helmholtz free energy wrt temperature * @@ -1190,7 +1190,7 @@ double WaterPropsIAPWSphi::enthalpy_RT() const { return hRT; } -/** +/* * Calculate the dimensionless entropy s/R. */ double WaterPropsIAPWSphi::entropy_R() const { @@ -1203,7 +1203,7 @@ double WaterPropsIAPWSphi::entropy_R() const { return sR; } -/** +/* * Calculate the dimensionless internal energy, u/RT. */ double WaterPropsIAPWSphi::intEnergy_RT() const { @@ -1214,7 +1214,7 @@ double WaterPropsIAPWSphi::intEnergy_RT() const { return uR; } -/** +/* * Calculate the dimensionless constant volume Heat Capacity, Cv/R */ double WaterPropsIAPWSphi::cv_R() const { @@ -1225,7 +1225,7 @@ double WaterPropsIAPWSphi::cv_R() const { return cvR; } -/** +/* * Calculate the dimensionless constant pressure Heat Capacity, Cp/R */ double WaterPropsIAPWSphi::cp_R() const { diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.h b/Cantera/src/thermo/WaterPropsIAPWSphi.h index 9a9eeb463..a8e5653d5 100644 --- a/Cantera/src/thermo/WaterPropsIAPWSphi.h +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.h @@ -1,5 +1,7 @@ /** * @file WaterPropsIAPWSphi.h + * + * Lowest level of the classes which support a real water model. */ /* * Copywrite (2006) Sandia Corporation. Under the terms of @@ -13,31 +15,81 @@ #ifndef WATERPROPSIAPWSPHI_H #define WATERPROPSIAPWSPHI_H -/* +/*! + * the WaterPropsIAPSWSphi class support low level calls for + * the real description of water. + * * Units Note: This class works with reduced units exclusively. */ - class WaterPropsIAPWSphi { public: + + //! Base constructor WaterPropsIAPWSphi(); - /* - * Calculate the base phi's, recalculating the internal polynomials + //! Calculate the Phi function, which is the base function + /*! + * The phi functino is basically the helmholtz free energy + * Eqn. (6.4) + * All internal polynomials are recalculated. + * + * @param tau Dimensionless temperature = T_c/T + * @param delta Dimensionless density = delta = rho / Rho_c */ double phi(double tau, double delta); - double phi_d(double tau, double delta); + + //! Delta derivative of phi + /*! + * @param tau Dimensionless temperature = T_c/T + * @param delta Dimensionless density = delta = rho / Rho_c + */ + double phi_d(double tau, double delta); + + //! 2nd derivative of phi wrt delta + /*! + * @param tau Dimensionless temperature = T_c/T + * @param delta Dimensionless density = delta = rho / Rho_c + */ double phi_dd(double tau, double delta); + + //! First derivative of phi wrt tau + /*! + * @param tau Dimensionless temperature = T_c/T + * @param delta Dimensionless density = delta = rho / Rho_c + */ double phi_t(double tau, double delta); + + //! Second derivative of phi wrt tau + /*! + * @param tau Dimensionless temperature = T_c/T + * @param delta Dimensionless density = delta = rho / Rho_c + */ double phi_tt(double tau, double delta); + + //! Second derivative of phi wrt tau, then delta + /*! + * @param tau Dimensionless temperature = T_c/T + * @param delta Dimensionless density = delta = rho / Rho_c + */ double phi_dt(double tau, double delta); + //! Internal check # 1 void check1(); + + //! Internal check # 2 void check2(); - /** - * Calculate the dimensionless pressure, pred: - * pred = pressure M / (rho RT) + + //! Calculate the dimensionless pressure at tau and delta; + /*! + * + * p/(rhoRT) = delta * phi_d() = 1.0 + delta phiR_d() + * + * @param tau Dimensionless temperature = T_c/T + * @param delta Dimensionless density = delta = rho / Rho_c + * + * note: this is done so much, we have a seperate routine. */ double pressure_rhoRT(double tau, double delta); @@ -46,6 +98,13 @@ public: * and the reduced temperature, tau. It takes an initial guess, deltaGuess. * DeltaGuess is important as this is a multivalued function below the * critical point. + * + * @param p_red Value of the dimensionless pressure + * @param tau Dimensionless temperature = T_c/T + * @param deltaGuess Initial guess for the dimensionless density + * + * @return + * Returns the dimensionless density. */ double dfind(double p_red, double tau, double deltaGuess); @@ -83,29 +142,65 @@ public: * 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 + * @param delta Dimensionless density = delta = rho / Rho_c */ void tdpolycalc(double tau, double delta); + //! Return the value of phiR(), res double phiR() const; + private: + + //! nau calculation double phi0() const; + //! calculation of d_phiR/d_d double phiR_d() const; + //! calculation of d_nau/d_d double phi0_d() const; + //! calculation of d2_res/d_dd double phiR_dd() const; + //! calculation of d2_nau/d_dd double phi0_dd() const; + //! calculation of d_nau/d_t double phi0_t() const; + //! calculation of d_res/d_t double phiR_t() const; + //! calculation of d2_res/d_tt double phiR_tt() const; + //! calculation of d2_nau/d_tt double phi0_tt() const; + //! calculation of d2_res/d_dt double phiR_dt() const; + //! calculation of d2_nau/d_dt double phi0_dt() const; + + /** + * intCheck() calculates all of the functions at a one point and + * prints out the result. It's used for conducting the internal + * check. + * + * @param tau Dimensionless temperature = T_c/T + * @param delta Dimensionless density = delta = rho / Rho_c + */ void intCheck(double tau, double delta); protected: + + //! Value of internally calculated polynomials of powers of TAU double TAUp[52]; + + //! Value of internally calculated polynomials of powers of delta double DELTAp[16]; + + //! Last tau that was used to calculate polynomials double TAUsave; + + //! sqrt of TAU double TAUsqrt; + + //! Last delta that was used to calculate polynomials double DELTAsave; }; #endif diff --git a/Cantera/src/units.h b/Cantera/src/units.h index a418c30e0..21ee851ec 100644 --- a/Cantera/src/units.h +++ b/Cantera/src/units.h @@ -129,9 +129,25 @@ namespace Cantera { /// pointer to the single instance of Unit static Unit* s_u; - std::map m_u; - std::map m_act_u; + //! Map between a string and a units double value + /*! + * This map maps the dimension string to the units value adjustment. Example + * - m_u["m"] = 1.0; + * - m_u["cm"] = 0.01; + */ + std::map m_u; + //! Map between a string and a units double value for activation energy units + /*! + * This map maps the dimension string to the units value adjustment. Example + * - m_act_u["K"] = GasConstant; + */ + std::map m_act_u; + + /*! + * Units class constructor, containing the default mappings between + * strings and units. + */ Unit(){ // length diff --git a/tools/doc/Cantera.cfg.in b/tools/doc/Cantera.cfg.in index 5b11be233..b10a2ac66 100755 --- a/tools/doc/Cantera.cfg.in +++ b/tools/doc/Cantera.cfg.in @@ -43,7 +43,7 @@ SUBGROUPING = YES # Build related configuration options #--------------------------------------------------------------------------- EXTRACT_ALL = NO -EXTRACT_PRIVATE = NO +EXTRACT_PRIVATE = YES EXTRACT_STATIC = NO EXTRACT_LOCAL_CLASSES = YES EXTRACT_LOCAL_METHODS = NO @@ -117,7 +117,8 @@ FILE_PATTERNS = Kinetics.h Kinetics.cpp \ IdealSolidSolnPhase.h IdealSolidSolnPhase.cpp \ StoichSubstanceSSTP.h StoichSubstanceSSTP.cpp \ DebyeHuckel.h DebyeHuckel.cpp \ - ConstDensityThermo.h ConstDensityThermo.cpp + ConstDensityThermo.h ConstDensityThermo.cpp \ + WaterPropsIAPWSphi.h WaterPropsIAPWSphi.cpp WaterPropsIAPWS.h WaterPropsIAPWS.cpp RECURSIVE = NO EXCLUDE = CVS examples converters zeroD EXCLUDE_SYMLINKS = NO