From 96930d57d6fcc72125de61a33ad622e8cf0b103e Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 21 Aug 2008 01:46:11 +0000 Subject: [PATCH] Incremental update on HKFT implementation Added thermalExpansionCoeff() and its temp derivative to objects. --- Cantera/src/thermo/HKFT_PDSS.cpp | 77 ++++++-- Cantera/src/thermo/HKFT_PDSS.h | 1 + Cantera/src/thermo/PDSS.cpp | 13 +- Cantera/src/thermo/PDSS.h | 9 + Cantera/src/thermo/WaterPDSS.cpp | 22 +++ Cantera/src/thermo/WaterPDSS.h | 304 ++++++++++++++++--------------- Cantera/src/thermo/WaterSSTP.cpp | 15 +- Cantera/src/thermo/WaterSSTP.h | 10 + 8 files changed, 293 insertions(+), 158 deletions(-) diff --git a/Cantera/src/thermo/HKFT_PDSS.cpp b/Cantera/src/thermo/HKFT_PDSS.cpp index 31f08870b..6e4f79482 100644 --- a/Cantera/src/thermo/HKFT_PDSS.cpp +++ b/Cantera/src/thermo/HKFT_PDSS.cpp @@ -356,11 +356,13 @@ namespace Cantera { double t2 = temp * temp; double val = ag_coeff[0] + ag_coeff[1] * temp + ag_coeff[2] * t2; return val; + } else if (ifunc == 1) { + return ag_coeff[1] + ag_coeff[2] * 2.0 * temp; } - if (ifunc != 1) { - throw CanteraError("HKFT_PDSS::ag", "unimplemented"); + if (ifunc != 2) { + return 0.0; } - return ag_coeff[1] + ag_coeff[2] * 2.0 * temp; + return ag_coeff[2] * 2.0;; } //! Internal formula for the calculation of b_g() @@ -373,11 +375,13 @@ namespace Cantera { double t2 = temp * temp; double val = bg_coeff[0] + bg_coeff[1] * temp + bg_coeff[2] * t2; return val; + } else if (ifunc == 1) { + return bg_coeff[1] + bg_coeff[2] * 2.0 * temp; } - if (ifunc != 1) { - throw CanteraError("HKFT_PDSS::bg", "unimplemented"); + if (ifunc != 2) { + return 0.0; } - return bg_coeff[1] + bg_coeff[2] * 2.0 * temp; + return bg_coeff[2] * 2.0; } double HKFT_PDSS::f(const double temp, const double pres, const int ifunc) { @@ -402,7 +406,14 @@ namespace Cantera { fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0); return fac1 * fac2; } else if (ifunc == 1) { - fac1 = (4.8 * pow(T1,3.8) + 16.0 * af_coeff[0] * pow(T1, 16.0)) / 300.; + fac1 = (4.8 * pow(T1,3.8) + 16.0 * af_coeff[0] * pow(T1, 15.0)) / 300.; + return fac1 * fac2; + } else if (ifunc == 2) { + fac1 = (4.8 * 3.8 * pow(T1,2.8) + 16.0 * 15.0 * af_coeff[0] * pow(T1, 14.0)) / (300. * 300.); + return fac1 * fac2; + } else if (ifunc == 3) { + fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0); + fac2 = (3.0 * af_coeff[1] * p2 + 4.0 * af_coeff[2] * p3 )/ 1.0E5; return fac1 * fac2; } else { throw CanteraError("HKFT_PDSS::gg", "unimplemented"); @@ -417,18 +428,58 @@ namespace Cantera { m_densWaterSS = m_waterSS->density(); // density in gm cm-3 double dens = m_densWaterSS * 1.0E-3; - if (ifunc == 0) { - if (dens >= 1.0) { - return 0.0; + double gval = afunc * pow((1.0-dens), bfunc); + if (dens >= 1.0) { + return 0.0; + } + if (ifunc == 0) { + return gval; + + } else if (ifunc == 1 || ifunc == 2) { + double afuncdT = ag(temp, 1); + double bfuncdT = bg(temp, 1); + double alpha = m_waterSS->thermalExpansionCoeff(); + + double fac1 = afuncdT * gval / afunc; + double fac2 = bfuncdT * gval * log(1.0 - dens); + double fac3 = gval * alpha * bfunc * dens / (1.0 - dens); + + double dgdt = fac1 + fac2 + fac3; + if (ifunc == 1) { + return dgdt; } + + double afuncdT2 = ag(temp, 2); + double bfuncdT2 = bg(temp, 2); + + double dfac1dT = dgdt * afuncdT / afunc + afuncdT2 * gval / afunc + - afuncdT * afuncdT * gval / (afunc * afunc); + + double ddensdT = - alpha * dens; + double dfac2dT = bfuncdT2 * gval * log(1.0 - dens) + + bfuncdT * dgdt * log(1.0 - dens) + - bfuncdT * gval /(1.0 - dens) * ddensdT; + + double dalphadT = m_waterSS->dthermalExpansionCoeffdT(); - double gval = afunc * pow((1.0-dens), bfunc); - double fval = f(temp, pres, ifunc); - return gval - fval; + double dfac3dT = dgdt * alpha * bfunc * dens / (1.0 - dens) + + gval * dalphadT * bfunc * dens / (1.0 - dens) + + gval * alpha * bfuncdT * dens / (1.0 - dens) + + gval * alpha * bfunc * ddensdT / (1.0 - dens) + - gval * alpha * bfunc * dens / ((1.0 - dens) * (1.0 - dens)) * ddensdT; + + return dfac1dT + dfac2dT + dfac3dT; + } else { throw CanteraError("HKFT_PDSS::gg", "unimplemented"); } return 0.0; } + double HKFT_PDSS::gstar(const double temp, const double pres, const int ifunc) { + double gval = g(temp, pres, ifunc); + double fval = f(temp, pres, ifunc); + return gval - fval; + } + } diff --git a/Cantera/src/thermo/HKFT_PDSS.h b/Cantera/src/thermo/HKFT_PDSS.h index f9348f43e..4a30d73db 100644 --- a/Cantera/src/thermo/HKFT_PDSS.h +++ b/Cantera/src/thermo/HKFT_PDSS.h @@ -144,6 +144,7 @@ namespace Cantera { double bg(const double temp, const int ifunc = 0) const; double g(const double temp, const double pres, const int ifunc = 0); double f(const double temp, const double pres, const int ifunc = 0); + double gstar(const double temp, const double pres, const int ifunc = 0); protected: diff --git a/Cantera/src/thermo/PDSS.cpp b/Cantera/src/thermo/PDSS.cpp index ad80d58e0..e16b55e39 100644 --- a/Cantera/src/thermo/PDSS.cpp +++ b/Cantera/src/thermo/PDSS.cpp @@ -397,7 +397,18 @@ namespace Cantera { throw CanteraError("PDSS::pressure()", "unimplemented"); } - + // Return the volumetric thermal expansion coefficient. Units: 1/K. + /* + * The thermal expansion coefficient is defined as + * \f[ + * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P + * \f] + */ + doublereal PDSS::thermalExpansionCoeff() const { + throw CanteraError("PDSS::thermalExpansionCoeff()", "unimplemented"); + return (0.0); + } + /// critical temperature doublereal PDSS::critTemperature() const { throw CanteraError("PDSS::critTemperature()", "unimplemented"); diff --git a/Cantera/src/thermo/PDSS.h b/Cantera/src/thermo/PDSS.h index 487f75bf8..0ec40fdbb 100644 --- a/Cantera/src/thermo/PDSS.h +++ b/Cantera/src/thermo/PDSS.h @@ -213,6 +213,15 @@ namespace Cantera { virtual doublereal pressure() const; virtual void setPressure(doublereal p); + //! Return the volumetric thermal expansion coefficient. Units: 1/K. + /*! + * The thermal expansion coefficient is defined as + * \f[ + * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P + * \f] + */ + virtual doublereal thermalExpansionCoeff() const; + //@} /// @name Partial Molar Properties of the Solution ----------------- //@{ diff --git a/Cantera/src/thermo/WaterPDSS.cpp b/Cantera/src/thermo/WaterPDSS.cpp index 1dff633f1..81479a8f7 100644 --- a/Cantera/src/thermo/WaterPDSS.cpp +++ b/Cantera/src/thermo/WaterPDSS.cpp @@ -387,6 +387,28 @@ namespace Cantera { } + // Return the volumetric thermal expansion coefficient. Units: 1/K. + /* + * The thermal expansion coefficient is defined as + * \f[ + * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P + * \f] + */ + doublereal WaterPDSS::thermalExpansionCoeff() const { + doublereal pres = pressure(); + doublereal val = m_sub->coeffThermExp(m_temp, pres); + return val; + } + + doublereal WaterPDSS::dthermalExpansionCoeffdT() const { + doublereal pres = pressure(); + double tt = m_temp - 0.04; + doublereal vald = m_sub->coeffThermExp(tt, pres); + doublereal val2 = m_sub->coeffThermExp(m_temp, pres); + doublereal val = (val2 - vald) / 0.04; + return val; + } + /// critical temperature doublereal WaterPDSS::critTemperature() const { return m_sub->Tcrit(); } diff --git a/Cantera/src/thermo/WaterPDSS.h b/Cantera/src/thermo/WaterPDSS.h index 5b3cfdfef..40673013a 100644 --- a/Cantera/src/thermo/WaterPDSS.h +++ b/Cantera/src/thermo/WaterPDSS.h @@ -25,178 +25,196 @@ class WaterPropsIAPWS; namespace Cantera { + /** + * Class for the liquid water pressure dependent + * standard state + * + * + * 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 + * 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 + * + * (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 + * and then use the theoretical ideal gas results to scale up to + * higher pressures: + * + * Ho(1bar) = H(P0) + * + * 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. + * + * + */ + class WaterPDSS : public PDSS { + + public: + /** - * Class for the liquid water pressure dependent - * standard state - * - * - * 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 - * 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 - * - * (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 - * and then use the theoretical ideal gas results to scale up to - * higher pressures: - * - * Ho(1bar) = H(P0) - * - * 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. - * - * + * Basic list of constructors and duplicators */ - class WaterPDSS : public PDSS { - - public: - - /** - * Basic list of constructors and duplicators - */ - WaterPDSS(); - WaterPDSS(ThermoPhase *tp, int spindex); - WaterPDSS(const WaterPDSS &b); - WaterPDSS& operator=(const WaterPDSS&b); - WaterPDSS(ThermoPhase *tp, int spindex, - std::string inputFile, std::string id = ""); - WaterPDSS(ThermoPhase *tp, int spindex, - XML_Node& phaseRef, std::string id = ""); - virtual ~WaterPDSS(); + WaterPDSS(); + WaterPDSS(ThermoPhase *tp, int spindex); + WaterPDSS(const WaterPDSS &b); + WaterPDSS& operator=(const WaterPDSS&b); + WaterPDSS(ThermoPhase *tp, int spindex, + std::string inputFile, std::string id = ""); + WaterPDSS(ThermoPhase *tp, int spindex, + XML_Node& phaseRef, std::string id = ""); + virtual ~WaterPDSS(); - /** - * - * @name Utilities - * @{ - */ - virtual int pdssType() const { return -1; } + /** + * + * @name Utilities + * @{ + */ + virtual int pdssType() const { return -1; } - /** - * @} - * @name Molar Thermodynamic Properties of the Solution -------------- - * @{ - */ - virtual doublereal enthalpy_mole() const; - virtual doublereal intEnergy_mole() const; - virtual doublereal entropy_mole() const; - virtual doublereal gibbs_mole() const; - virtual doublereal cp_mole() const; - virtual doublereal cv_mole() const; + /** + * @} + * @name Molar Thermodynamic Properties of the Solution -------------- + * @{ + */ + virtual doublereal enthalpy_mole() const; + virtual doublereal intEnergy_mole() const; + virtual doublereal entropy_mole() const; + virtual doublereal gibbs_mole() const; + virtual doublereal cp_mole() const; + virtual doublereal cv_mole() const; - //@} - /// @name Mechanical Equation of State Properties --------------------- - //@{ + //@} + /// @name Mechanical Equation of State Properties --------------------- + //@{ - virtual doublereal pressure() const; - virtual void setTempPressure(doublereal t, doublereal p); - virtual void setPressure(doublereal p); + virtual doublereal pressure() const; + virtual void setTempPressure(doublereal t, doublereal p); + virtual void setPressure(doublereal p); - //@} - /// @name Partial Molar Properties of the Solution ----------------- - //@{ + //! Return the volumetric thermal expansion coefficient. Units: 1/K. + /*! + * The thermal expansion coefficient is defined as + * \f[ + * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P + * \f] + */ + virtual doublereal thermalExpansionCoeff() const; + + //! Return the derivative of the volumetric thermal expansion coefficient. Units: 1/K2. + /*! + * The thermal expansion coefficient is defined as + * \f[ + * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P + * \f] + */ + virtual doublereal dthermalExpansionCoeffdT() const; - virtual void getChemPotentials(doublereal* mu) const { - mu[0] = gibbs_mole(); - } + //@} + /// @name Partial Molar Properties of the Solution ----------------- + //@{ - //@} - /// @name Properties of the Standard State of the Species - // in the Solution -- - //@{ + virtual void getChemPotentials(doublereal* mu) const { + mu[0] = gibbs_mole(); + } + + //@} + /// @name Properties of the Standard State of the Species + // in the Solution -- + //@{ - /// critical temperature - virtual doublereal critTemperature() const; + /// critical temperature + virtual doublereal critTemperature() const; - /// critical pressure - virtual doublereal critPressure() const; + /// critical pressure + virtual doublereal critPressure() const; - /// critical density - virtual doublereal critDensity() const; + /// critical density + virtual doublereal critDensity() const; - /// saturation temperature - //virtual doublereal satTemperature(doublereal p) const; + /// saturation temperature + //virtual doublereal satTemperature(doublereal p) const; - /// saturation pressure - virtual doublereal satPressure(doublereal t); + /// saturation pressure + virtual doublereal satPressure(doublereal t); - virtual void setDensity(double dens); - double density() const; - virtual void setTemperature(double temp); - double temperature() const; - virtual void setState_TP(double temp, double pres); + virtual void setDensity(double dens); + double density() const; + virtual void setTemperature(double temp); + double temperature() const; + virtual void setState_TP(double temp, double pres); - doublereal molecularWeight() const; - void setMolecularWeight(double mw); + doublereal molecularWeight() const; + void setMolecularWeight(double mw); - virtual void constructPDSS(ThermoPhase *tp, int spindex); - virtual void constructPDSSFile(ThermoPhase *tp, int spindex, - std::string inputFile, std::string id); - virtual void constructPDSSXML(ThermoPhase *tp, int spindex, - XML_Node& phaseNode, std::string id); - virtual void initThermoXML(XML_Node& eosdata, std::string id); - virtual void initThermo(); - virtual void setParametersFromXML(const XML_Node& eosdata); - WaterPropsIAPWS *getWater() const { - return m_sub; - } -protected: + virtual void constructPDSS(ThermoPhase *tp, int spindex); + virtual void constructPDSSFile(ThermoPhase *tp, int spindex, + std::string inputFile, std::string id); + virtual void constructPDSSXML(ThermoPhase *tp, int spindex, + XML_Node& phaseNode, std::string id); + virtual void initThermoXML(XML_Node& eosdata, std::string id); + virtual void initThermo(); + virtual void setParametersFromXML(const XML_Node& eosdata); + WaterPropsIAPWS *getWater() const { + return m_sub; + } + protected: -private: - mutable WaterPropsIAPWS *m_sub; + private: + mutable WaterPropsIAPWS *m_sub; - /** - * state of the system (temperature and density); - */ - doublereal m_temp; - doublereal m_dens; + /** + * state of the system (temperature and density); + */ + doublereal m_temp; + doublereal m_dens; - /* - * state of the fluid - * 0 gas - * 1 liquid - * 2 supercrit - */ - int m_iState; + /* + * state of the fluid + * 0 gas + * 1 liquid + * 2 supercrit + */ + int m_iState; - /** - * Offset constants used to obtain consistency with the NIST database. - * This is added to all internal energy and enthalpy results. - * units = J kmol-1. - */ - double EW_Offset; + /** + * Offset constants used to obtain consistency with the NIST database. + * This is added to all internal energy and enthalpy results. + * units = J kmol-1. + */ + double EW_Offset; - /* - * Offset constant used to obtain consistency with NIST convention. - * This is added to all internal entropy results. - * units = J kmol-1 K-1. - */ - double SW_Offset; + /* + * Offset constant used to obtain consistency with NIST convention. + * This is added to all internal entropy results. + * units = J kmol-1 K-1. + */ + double SW_Offset; - bool m_verbose; + bool m_verbose; - /** - * Since this phase represents a liquid phase, it's an error to - * return a gas-phase answer. However, if the below is true, then - * a gas-phase answer is allowed. This is used to check the thermodynamic - * consistency with ideal-gas thermo functions for example. - */ - bool m_allowGasPhase; - }; + /** + * Since this phase represents a liquid phase, it's an error to + * return a gas-phase answer. However, if the below is true, then + * a gas-phase answer is allowed. This is used to check the thermodynamic + * consistency with ideal-gas thermo functions for example. + */ + bool m_allowGasPhase; + }; } diff --git a/Cantera/src/thermo/WaterSSTP.cpp b/Cantera/src/thermo/WaterSSTP.cpp index bf79f4e4d..72db00445 100644 --- a/Cantera/src/thermo/WaterSSTP.cpp +++ b/Cantera/src/thermo/WaterSSTP.cpp @@ -499,8 +499,21 @@ namespace Cantera { } setDensity(dd); } + + // Return the volumetric thermal expansion coefficient. Units: 1/K. + /* + * The thermal expansion coefficient is defined as + * \f[ + * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P + * \f] + */ + doublereal WaterSSTP::thermalExpansionCoeff() const { + doublereal pres = pressure(); + double T = temperature(); + doublereal val = m_sub->coeffThermExp(T, pres); + return val; + } - // critical temperature doublereal WaterSSTP::critTemperature() const { return m_sub->Tcrit(); } diff --git a/Cantera/src/thermo/WaterSSTP.h b/Cantera/src/thermo/WaterSSTP.h index 3566760ce..c2fcbb597 100644 --- a/Cantera/src/thermo/WaterSSTP.h +++ b/Cantera/src/thermo/WaterSSTP.h @@ -186,6 +186,16 @@ namespace Cantera { virtual doublereal pressure() const; virtual void setPressure(doublereal p); + //! Return the volumetric thermal expansion coefficient. Units: 1/K. + /*! + * The thermal expansion coefficient is defined as + * \f[ + * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P + * \f] + */ + virtual doublereal thermalExpansionCoeff() const; + + /** * @} * @name Potential Energy