diff --git a/Cantera/src/thermo/IdealMolalSoln.cpp b/Cantera/src/thermo/IdealMolalSoln.cpp index 15c3f3722..b7f7316bf 100644 --- a/Cantera/src/thermo/IdealMolalSoln.cpp +++ b/Cantera/src/thermo/IdealMolalSoln.cpp @@ -3,7 +3,7 @@ * @file IdealMolalSoln.cpp */ /* - * Copywrite (2005) Sandia Corporation. Under the terms of + * Copywrite (2006) Sandia Corporation. Under the terms of * Contract DE-AC04-94AL85000 with Sandia Corporation, the * U.S. Government retains certain rights in this software. */ diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index c986363f9..05ed803e2 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -24,10 +24,15 @@ CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) ifeq ($(do_electro),1) ELECTRO_OBJ = SingleSpeciesTP.o StoichSubstanceSSTP.o \ MolalityVPSSTP.o VPStandardStateTP.o \ - IdealSolidSolnPhase.o IdealMolalSoln.o + IdealSolidSolnPhase.o IdealMolalSoln.o \ + WaterPropsIAPWSphi.o WaterPropsIAPWS.o WaterProps.o \ + PDSS.o WaterPDSS.o WaterTP.o + ELECTRO_H = SingleSpeciesTP.h StoichSubstanceSSTP.h \ MolalityVPSSTP.h VPStandardStateTP.h \ - IdealSolidSolnPhase.h IdealMolalSoln.h + IdealSolidSolnPhase.h IdealMolalSoln.h \ + WaterPropsIAPWSphi.h WaterPropsIAPWS.h WaterProps.h \ + PDSS.h WaterPDSS.h WaterTP.h endif ifeq ($(do_issp),1) ISSP_OBJ = IdealSolidSolnPhase.o diff --git a/Cantera/src/thermo/PDSS.cpp b/Cantera/src/thermo/PDSS.cpp new file mode 100644 index 000000000..a504f42b6 --- /dev/null +++ b/Cantera/src/thermo/PDSS.cpp @@ -0,0 +1,391 @@ +/** + * @file PDSS.cpp + * + * Implementation of a pressure dependent standard state + * virtual function. + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ + +#include "ct_defs.h" +#include "xml.h" +#include "ctml.h" +#include "PDSS.h" +#include "importCTML.h" + + +#include "ThermoPhase.h" + +namespace Cantera { + /** + * Basic list of constructors and duplicators + */ + + PDSS::PDSS(ThermoPhase *tp, int spindex) : + m_temp(-1.0), + m_dens(-1.0), + m_tp(tp), + m_spindex(spindex), + m_mw(0.0) + { + constructPDSS(tp, spindex); + } + + + PDSS::PDSS(ThermoPhase *tp, int spindex, string inputFile, string id) : + m_temp(-1.0), + m_dens(-1.0), + m_tp(tp), + m_spindex(spindex), + m_mw(0.0) + { + constructPDSSFile(tp, spindex, inputFile, id); + } + + + PDSS::PDSS(ThermoPhase *tp, int spindex, XML_Node& phaseRoot, string id) : + m_temp(-1.0), + m_dens(-1.0), + m_tp(0), + m_spindex(0), + m_mw(0.0) + { + constructPDSSXML(tp, spindex, phaseRoot, id) ; + } + + + PDSS::PDSS(const PDSS &b) : + m_temp(-1.0), + m_dens(-1.0), + m_tp(0), + m_spindex(0), + m_mw(b.m_mw) + { + /* + * Use the assignment operator to do the brunt + * of the work for the copy construtor. + */ + *this = b; + } + + /** + * Assignment operator + */ + PDSS& PDSS::operator=(const PDSS&b) { + if (&b == this) return *this; + m_tp = b.m_tp; + m_spindex = b.m_spindex; + m_temp = b.m_temp; + m_dens = b.m_dens; + m_mw = b.m_mw; + return *this; + } + + PDSS::~PDSS() { + } + + void PDSS::constructPDSS(ThermoPhase *tp, int spindex) { + initThermo(); + } + + + /** + * constructPDSSXML: + * + * Initialization of a PDSS object using an + * xml file. + * + * This routine is a precursor to initThermo(XML_Node*) + * routine, which does most of the work. + * + * @param infile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void PDSS::constructPDSSXML(ThermoPhase *tp, int spindex, + XML_Node& phaseNode, string id) { + initThermo(); + } + + + /** + * constructPDSSFile(): + * + * Initialization of a PDSS object using an + * xml file. + * + * This routine is a precursor to initThermo(XML_Node*) + * routine, which does most of the work. + * + * @param infile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void PDSS::constructPDSSFile(ThermoPhase *tp, int spindex, + string inputFile, string id) { + + if (inputFile.size() == 0) { + throw CanteraError("PDSS::initThermo", + "input file is null"); + } + string path = findInputFile(inputFile); + ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("PDSS::initThermo","could not open " + +path+" for reading."); + } + /* + * The phase object automatically constructs an XML object. + * Use this object to store information. + */ + + XML_Node *fxml = new XML_Node(); + fxml->build(fin); + XML_Node *fxml_phase = findXMLPhase(fxml, id); + if (!fxml_phase) { + throw CanteraError("PDSS::initThermo", + "ERROR: Can not find phase named " + + id + " in file named " + inputFile); + } + constructPDSSXML(tp, spindex, *fxml_phase, id); + delete fxml; + } + + void PDSS:: + initThermoXML(XML_Node& phaseNode, string id) { + initThermo(); + } + + void PDSS::initThermo() { + } + + void PDSS:: + setParametersFromXML(const XML_Node& eosdata) { + } + + /** + * Return the molar enthalpy in units of J kmol-1 + */ + doublereal PDSS:: + enthalpy_mole() const { + throw CanteraError("PDSS::enthalpy_mole()", "unimplemented"); + return (0.0); + } + + /** + * Calculate the internal energy in mks units of + * J kmol-1 + */ + doublereal PDSS:: + intEnergy_mole() const { + throw CanteraError("PDSS::enthalpy_mole()", "unimplemented"); + return (0.0); + } + + /** + * Calculate the entropy in mks units of + * J kmol-1 K-1 + */ + doublereal PDSS:: + entropy_mole() const { + + throw CanteraError("PDSS::entropy_mole()", "unimplemented"); + return (0.0); + } + + /** + * Calculate the Gibbs free energy in mks units of + * J kmol-1 K-1. + */ + doublereal PDSS:: + gibbs_mole() const { + throw CanteraError("PDSS::gibbs_mole()", "unimplemented"); + return (0.0); + } + + /** + * Calculate the constant pressure heat capacity + * in mks units of J kmol-1 K-1 + */ + doublereal PDSS:: + cp_mole() const { + throw CanteraError("PDSS::cp_mole()", "unimplemented"); + return (0.0); + } + + /** + * Calculate the constant volume heat capacity + * in mks units of J kmol-1 K-1 + */ + doublereal PDSS:: + cv_mole() const { + throw CanteraError("PDSS::cv_mole()", "unimplemented"); + return (0.0); + } + + + /** + * Return the difference in enthalpy between current p + * and ref p0, in mks units of + * in units of J kmol-1 + */ + doublereal PDSS:: + enthalpyDelp_mole() const { + throw CanteraError("PDSS::enthalpy_mole()", "unimplemented"); + return (0.0); + } + + /** + * Calculate difference in the internal energy between current p + * and ref p0, in mks units of + * J kmol-1 + */ + doublereal PDSS:: + intEnergyDelp_mole() const { + throw CanteraError("PDSS::enthalpyDelp_mole()", "unimplemented"); + return (0.0); + } + + /** + * Return the difference in entropy between current p + * and ref p0, in mks units of + * J kmol-1 K-1 + */ + doublereal PDSS:: + entropyDelp_mole() const { + + throw CanteraError("PDSS::entropyDelp_mole()", "unimplemented"); + return (0.0); + } + + /** + * Calculate the difference in Gibbs free energy between current p and + * the ref p0, in mks units of + * J kmol-1 K-1. + */ + doublereal PDSS:: + gibbsDelp_mole() const { + throw CanteraError("PDSS::gibbsDelp_mole()", "unimplemented"); + return (0.0); + } + + /** + * Calculate the difference in the constant pressure heat capacity + * between the current p and the ref p0, + * in mks units of J kmol-1 K-1 + */ + doublereal PDSS:: + cpDelp_mole() const { + throw CanteraError("PDSS::cpDelp_mole()", "unimplemented"); + return (0.0); + } + + /** + * Calculate the difference in constant volume heat capacity + * between the current p and the ref p0 + * in mks units of J kmol-1 K-1 + */ + doublereal PDSS:: + cvDelp_mole() const { + throw CanteraError("PDSS::cvDelp_mole()", "unimplemented"); + return (0.0); + } + + /** + * Calculate the pressure (Pascals), given the temperature and density + * Temperature: kelvin + * rho: density in kg m-3 + */ + doublereal PDSS:: + pressure() const { + throw CanteraError("PDSS::pressure()", "unimplemented"); + return (0.0); + } + + void PDSS:: + setPressure(doublereal p) { + throw CanteraError("PDSS::pressure()", "unimplemented"); + } + + + /// critical temperature + doublereal PDSS::critTemperature() const { + throw CanteraError("PDSS::critTemperature()", "unimplemented"); + return (0.0); + } + + /// critical pressure + doublereal PDSS::critPressure() const { + throw CanteraError("PDSS::critPressure()", "unimplemented"); + return (0.0); + } + + /// critical density + doublereal PDSS::critDensity() const { + throw CanteraError("PDSS::critDensity()", "unimplemented"); + return (0.0); + } + + void PDSS::setDensity(double dens) { + m_dens = dens; + } + + /** + * Return the density of the standard state + * + * We assume that the storred density is current. + * Note, this is the density of the standard state, + * not of the mixture. + */ + double PDSS::density() const { + return m_dens; + } + + /** + * Return the temperature + * + * Obtain the temperature from the owning ThermoPhase object + * if you can. + */ + double PDSS::temperature() const { + if (m_tp) { + m_temp = m_tp->temperature(); + } + return m_temp; + } + + void PDSS::setTemperature(double temp) { + m_temp = temp; + } + + doublereal PDSS::molecularWeight() const { + return m_mw; + } + void PDSS::setMolecularWeight(double mw) { + m_mw = mw; + } + + void PDSS::setState_TP(double temp, double pres) { + throw CanteraError("PDSS::setState_TP()", "unimplemented"); + } + + /// saturation pressure + doublereal PDSS::satPressure(doublereal t){ + throw CanteraError("PDSS::satPressure()", "unimplemented"); + return (0.0); + } + + +} diff --git a/Cantera/src/thermo/PDSS.h b/Cantera/src/thermo/PDSS.h new file mode 100644 index 000000000..097f42b40 --- /dev/null +++ b/Cantera/src/thermo/PDSS.h @@ -0,0 +1,173 @@ +/** + * @file PDSS.h + * + * Declares class PDSS pressure dependent standard state + * for a single species + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ + +#ifndef CT_PDSS_H +#define CT_PDSS_H +#include "ct_defs.h" + +class XML_Node; +class ThermoPhase; + +class WaterPropsIAPWS; + +namespace Cantera { + + + /** + * Class for pressure dependent standard states. + * + * + */ + class PDSS { + + public: + + /** + * Basic list of constructors and duplicators + */ + PDSS(ThermoPhase *tp, int spindex); + PDSS(const PDSS &b); + PDSS& operator=(const PDSS&b); + PDSS(ThermoPhase *tp, int spindex, string inputFile, string id = ""); + PDSS(ThermoPhase *tp, int spindex, XML_Node& phaseRef, string id = ""); + virtual ~PDSS(); + + /** + * + * @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; + + /* + * Get the difference in the standard state thermodynamic properties + * between the reference pressure, po, and the current pressure. + */ + virtual doublereal enthalpyDelp_mole() const; + virtual doublereal intEnergyDelp_mole() const; + virtual doublereal entropyDelp_mole() const; + virtual doublereal gibbsDelp_mole() const; + virtual doublereal cpDelp_mole() const; + virtual doublereal cvDelp_mole() const; + + //@} + /// @name Mechanical Equation of State Properties --------------------- + //@{ + + virtual doublereal pressure() const; + virtual void setPressure(doublereal p); + + //@} + /// @name Partial Molar Properties of 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 pressure + virtual doublereal critPressure() const; + + /// critical density + virtual doublereal critDensity() const; + + /// saturation temperature + //virtual doublereal satTemperature(doublereal p) const; + + + + /// 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); + + doublereal molecularWeight() const; + void setMolecularWeight(double mw); + + virtual void constructPDSS(ThermoPhase *tp, int spindex); + virtual void constructPDSSFile(ThermoPhase *tp, int spindex, + string inputFile, string id); + virtual void constructPDSSXML(ThermoPhase *tp, int spindex, + XML_Node& phaseNode, string id); + virtual void initThermoXML(XML_Node& eosdata, string id); + virtual void initThermo(); + virtual void setParametersFromXML(const XML_Node& eosdata); + + protected: + + /** + * state of the system (temperature and density); + * This may redundant and may go away. Should be able to + * get this information from owning ThermoPhase object. + */ + mutable doublereal m_temp; + + /** + * state of the system (temperature and density); + * This may redundant and may go away. Should be able to + * get this information from owning ThermoPhase object. + */ + doublereal m_dens; + + /** + * Thermophase which this species belongs to. Note, in some + * applications (i.e., mostly testing applications, this may be a null + * value. Applications should test whether this is null before usage. + */ + ThermoPhase *m_tp; + + /** + * Species index in the thermophase corresponding to this species. + */ + int m_spindex; + + /** + * Molecular Weight of the species + */ + doublereal m_mw; + + }; + +} + +#endif + + + diff --git a/Cantera/src/thermo/WaterPDSS.cpp b/Cantera/src/thermo/WaterPDSS.cpp new file mode 100644 index 000000000..09fe2badb --- /dev/null +++ b/Cantera/src/thermo/WaterPDSS.cpp @@ -0,0 +1,405 @@ +/** + * @file WaterPDSS.cpp + * + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ +#include "ct_defs.h" +#include "xml.h" +#include "ctml.h" +#include "WaterPDSS.h" +#include "WaterPropsIAPWS.h" +#include "importCTML.h" + + +#include "ThermoPhase.h" + +namespace Cantera { + /** + * Basic list of constructors and duplicators + */ + + WaterPDSS::WaterPDSS(ThermoPhase *tp, int spindex) : + PDSS(tp, spindex), + m_sub(0), + m_iState(-1), + EW_Offset(0.0), + SW_Offset(0.0), + m_verbose(0), + m_allowGasPhase(false) + { + constructPDSS(tp, spindex); + } + + + WaterPDSS::WaterPDSS(ThermoPhase *tp, int spindex, + string inputFile, string id) : + PDSS(tp, spindex), + m_sub(0), + m_iState(-1), + m_mw(0.0), + EW_Offset(0.0), + SW_Offset(0.0), + m_verbose(0), + m_allowGasPhase(false) + { + constructPDSSFile(tp, spindex, inputFile, id); + } + + WaterPDSS::WaterPDSS(ThermoPhase *tp, int spindex, + XML_Node& phaseRoot, string id) : + PDSS(tp, spindex), + m_sub(0), + m_iState(-1), + m_mw(0.0), + EW_Offset(0.0), + SW_Offset(0.0), + m_verbose(0), + m_allowGasPhase(false) + { + constructPDSSXML(tp, spindex, phaseRoot, id) ; + } + + + + WaterPDSS::WaterPDSS(const WaterPDSS &b) : + PDSS(b), + m_sub(0), + m_iState(-1), + m_mw(b.m_mw), + EW_Offset(b.EW_Offset), + SW_Offset(b.SW_Offset), + m_verbose(b.m_verbose), + m_allowGasPhase(b.m_allowGasPhase) + { + m_sub = new WaterPropsIAPWS(*(b.m_sub)); + /* + * Use the assignment operator to do the brunt + * of the work for the copy construtor. + */ + *this = b; + } + + /** + * Assignment operator + */ + WaterPDSS& WaterPDSS::operator=(const WaterPDSS&b) { + if (&b == this) return *this; + m_sub->operator=(*(b.m_sub)); + PDSS::operator=(b); + m_verbose = b.m_verbose; + m_allowGasPhase = b.m_allowGasPhase; + return *this; + } + + WaterPDSS::~WaterPDSS() { + delete m_sub; + } + + void WaterPDSS::constructPDSS(ThermoPhase *tp, int spindex) { + initThermo(); + } + + /** + * constructPDSSXML: + * + * Initialization of a Debye-Huckel phase using an + * xml file. + * + * This routine is a precursor to initThermo(XML_Node*) + * routine, which does most of the work. + * + * @param infile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void WaterPDSS::constructPDSSXML(ThermoPhase *tp, int spindex, + XML_Node& phaseNode, string id) { + initThermo(); + } + + /** + * constructPDSSFile(): + * + * Initialization of a Debye-Huckel phase using an + * xml file. + * + * This routine is a precursor to initThermo(XML_Node*) + * routine, which does most of the work. + * + * @param infile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void WaterPDSS::constructPDSSFile(ThermoPhase *tp, int spindex, + string inputFile, string id) { + + if (inputFile.size() == 0) { + throw CanteraError("WaterTp::initThermo", + "input file is null"); + } + string path = findInputFile(inputFile); + ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("WaterPDSS::initThermo","could not open " + +path+" for reading."); + } + /* + * The phase object automatically constructs an XML object. + * Use this object to store information. + */ + + XML_Node *fxml = new XML_Node(); + fxml->build(fin); + XML_Node *fxml_phase = findXMLPhase(fxml, id); + if (!fxml_phase) { + throw CanteraError("WaterPDSS::initThermo", + "ERROR: Can not find phase named " + + id + " in file named " + inputFile); + } + constructPDSSXML(tp, spindex, *fxml_phase, id); + delete fxml; + } + + void WaterPDSS:: + initThermoXML(XML_Node& phaseNode, string id) { + initThermo(); + } + + void WaterPDSS::initThermo() { + if (m_sub) delete m_sub; + m_sub = new WaterPropsIAPWS(); + if (m_sub == 0) { + throw CanteraError("WaterPDSS::initThermo", + "could not create new substance object."); + } + /* + * Calculate the molecular weight. + * hard coded to Cantera's elements and Water. + */ + m_mw = 2 * 1.00794 + 15.9994; + + /* + * Set the baseline + */ + doublereal T = 298.15; + + doublereal presLow = 1.0E-2; + doublereal oneBar = 1.0E5; + doublereal dens = 1.0E-9; + doublereal dd = m_sub->density(T, presLow, WATER_GAS, dens); + setTemperature(T); + m_dens = dd; + SW_Offset = 0.0; + doublereal s = entropy_mole(); + s -= GasConstant * log(oneBar/presLow); + if (s != 188.835E3) { + SW_Offset = 188.835E3 - s; + } + s = entropy_mole(); + s -= GasConstant * log(oneBar/presLow); + //printf("s = %g\n", s); + + doublereal h = enthalpy_mole(); + if (h != -241.826E6) { + EW_Offset = -241.826E6 - h; + } + h = enthalpy_mole(); + + //printf("h = %g\n", h); + + + /* + * Set the initial state of the system to 298.15 K and + * 1 bar. + */ + setTemperature(298.15); + double rho0 = m_sub->density(298.15, OneAtm, WATER_LIQUID); + m_dens = rho0; + + + } + + void WaterPDSS:: + setParametersFromXML(const XML_Node& eosdata) { + + } + + /** + * Return the molar enthalpy in units of J kmol-1 + */ + doublereal WaterPDSS:: + enthalpy_mole() const { + double T = m_temp; + double dens = m_dens; + doublereal h = m_sub->enthalpy(T, dens); + return (h + EW_Offset); + } + + /** + * Calculate the internal energy in mks units of + * J kmol-1 + */ + doublereal WaterPDSS:: + intEnergy_mole() const { + double T = m_dens; + double dens = m_temp; + doublereal u = m_sub->intEnergy(T, dens); + return (u + EW_Offset); + } + + /** + * Calculate the entropy in mks units of + * J kmol-1 K-1 + */ + doublereal WaterPDSS:: + entropy_mole() const { + double T = m_temp; + double dens = m_dens; + doublereal s = m_sub->entropy(T, dens); + return (s + SW_Offset); + } + + /** + * Calculate the Gibbs free energy in mks units of + * J kmol-1 K-1. + */ + doublereal WaterPDSS:: + gibbs_mole() const { + double T = m_temp; + double dens = m_dens; + doublereal g = m_sub->Gibbs(T, dens); + return (g + EW_Offset - SW_Offset*T); + } + + /** + * Calculate the constant pressure heat capacity + * in mks units of J kmol-1 K-1 + */ + doublereal WaterPDSS:: + cp_mole() const { + double T = m_temp; + double dens = m_dens; + doublereal cp = m_sub->cp(T, dens); + return cp; + } + + /** + * Calculate the constant volume heat capacity + * in mks units of J kmol-1 K-1 + */ + doublereal WaterPDSS:: + cv_mole() const { + double T = m_temp; + double dens = m_dens; + doublereal cv = m_sub->cv(T, dens); + return cv; + } + + /** + * Calculate the pressure (Pascals), given the temperature and density + * Temperature: kelvin + * rho: density in kg m-3 + */ + doublereal WaterPDSS:: + pressure() const { + double T = m_temp; + double dens = m_dens; + doublereal p = m_sub->pressure(T, dens); + return p; + } + + void WaterPDSS:: + setTempPressure(doublereal t, doublereal p) { + m_temp = t; + setPressure(p); + } + + void WaterPDSS:: + setPressure(doublereal p) { + double T = m_temp; + double dens = m_dens; + int waterState = WATER_GAS; + double rc = m_sub->Rhocrit(); + if (dens > rc) { + waterState = WATER_LIQUID; + } +#ifdef DEBUG_HKM + //printf("waterPDSS: set pres = %g t = %g, waterState = %d\n", + // p, T, waterState); +#endif + doublereal dd = m_sub->density(T, p, waterState, dens); + if (dd <= 0.0) { + printf("throw an error\n"); + throw CanteraError("WaterPDSS:pressure", "Failed to set water state"); + } + m_dens = dd; + } + + + /// critical temperature + doublereal WaterPDSS::critTemperature() const { return m_sub->Tcrit(); } + + /// critical pressure + doublereal WaterPDSS::critPressure() const { return m_sub->Pcrit(); } + + /// critical density + doublereal WaterPDSS::critDensity() const { return m_sub->Rhocrit(); } + + void WaterPDSS::setDensity(double dens) { + m_dens = dens; + m_sub->setState(m_temp, m_dens); + } + + double WaterPDSS::density() const { + return m_dens; + } + + + double WaterPDSS::temperature() const { + return m_temp; + } + + + void WaterPDSS::setTemperature(double temp) { + m_temp = temp; + doublereal dd = m_dens; + m_sub->setState(temp, dd); + } + + doublereal WaterPDSS::molecularWeight() const { + return m_mw; + } + void WaterPDSS::setMolecularWeight(double mw) { + m_mw = mw; + } + + void WaterPDSS::setState_TP(double temp, double pres) { + m_temp = temp; + setPressure(pres); + } + + /// saturation pressure + doublereal WaterPDSS::satPressure(doublereal t){ + doublereal pp = m_sub->psat(t); + double dens = m_dens; + m_temp = t; + m_dens = dens; + return pp; + } + + + +} diff --git a/Cantera/src/thermo/WaterPDSS.h b/Cantera/src/thermo/WaterPDSS.h new file mode 100644 index 000000000..2561ffbcb --- /dev/null +++ b/Cantera/src/thermo/WaterPDSS.h @@ -0,0 +1,220 @@ +/** + * @file WaterPDSS.h + * + * Declares class PureFluid + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* $Author$ + * $Date$ + * $Revision$ + */ + +#ifndef CT_WATERPDSS_H +#define CT_WATERPDSS_H +#include "ct_defs.h" +#include "PDSS.h" +class XML_Node; +#include "ThermoPhase.h" + +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: + + /** + * Basic list of constructors and duplicators + */ + WaterPDSS(ThermoPhase *tp, int spindex); + WaterPDSS(const WaterPDSS &b); + WaterPDSS& operator=(const WaterPDSS&b); + WaterPDSS(ThermoPhase *tp, int spindex, + string inputFile, string id = ""); + WaterPDSS(ThermoPhase *tp, int spindex, + XML_Node& phaseRef, string id = ""); + virtual ~WaterPDSS(); + + /** + * + * @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 Mechanical Equation of State Properties --------------------- + //@{ + + virtual doublereal pressure() const; + virtual void setTempPressure(doublereal t, doublereal p); + virtual void setPressure(doublereal p); + + //@} + /// @name Partial Molar Properties of 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 pressure + virtual doublereal critPressure() const; + + /// critical density + virtual doublereal critDensity() const; + + /// saturation temperature + //virtual doublereal satTemperature(doublereal p) const; + + + + /// 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); + + doublereal molecularWeight() const; + void setMolecularWeight(double mw); + + virtual void constructPDSS(ThermoPhase *tp, int spindex); + virtual void constructPDSSFile(ThermoPhase *tp, int spindex, + string inputFile, string id); + virtual void constructPDSSXML(ThermoPhase *tp, int spindex, + XML_Node& phaseNode, string id); + virtual void initThermoXML(XML_Node& eosdata, string id); + virtual void initThermo(); + virtual void setParametersFromXML(const XML_Node& eosdata); + WaterPropsIAPWS *getWater() const { + return m_sub; + } +protected: + + +private: + mutable WaterPropsIAPWS *m_sub; + + + /** + * 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; + + /** + * Thermophase which this species belongs to + */ + ThermoPhase *m_tp; + + /** + * Species index in the thermophase corresponding to this species. + */ + int m_spindex; + + /* + * Molecular Weight + */ + doublereal m_mw; + + /** + * 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; + + 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; + }; + +} + +#endif + + + diff --git a/Cantera/src/thermo/WaterProps.cpp b/Cantera/src/thermo/WaterProps.cpp new file mode 100644 index 000000000..69198bbc0 --- /dev/null +++ b/Cantera/src/thermo/WaterProps.cpp @@ -0,0 +1,435 @@ +/** + * @file WaterProps.cpp + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ + +#ifndef MAX +#define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) +#endif + + +#include "WaterProps.h" +#include "ctml.h" +#include "WaterPDSS.h" +#include "WaterPropsIAPWS.h" + + +namespace Cantera { + + + /* + * default constructor -> object owns its own water evaluator + */ + WaterProps::WaterProps(): + m_waterIAPWS(0), + m_own_sub(false) + { + m_waterIAPWS = new WaterPropsIAPWS(); + m_own_sub = true; + } + + /* + * constructor -> object in slave mode, It doesn't own its + * own water evaluator. + */ + WaterProps::WaterProps(WaterPDSS *wptr) : + m_waterIAPWS(0), + m_own_sub(false) + { + m_waterIAPWS = wptr->getWater(); + m_own_sub = false; + } + + /** + * Copy constructor + */ + WaterProps::WaterProps(const WaterProps &b) + { + *this = b; + } + + /** + * Destructor + */ + WaterProps::~WaterProps() { + if (m_own_sub) { + delete m_waterIAPWS; + } + } + + /** + * Assignment operator + */ + WaterProps& WaterProps::operator=(const WaterProps&b) { + if (&b == this) return *this; + /* + * add content here. + */ + + return *this; + } + + /* + * Simple calculation of water density at atmospheric pressure. + * Valid up to boiling point. + * + * ifunc = 0 Returns the density in kg/m^3 + * ifunc = 1 returns the derivative of the density wrt T. + * ifunc = 3 returns the derivative of the density wrt P + * ifunc = 2 returns the 2nd derivative of the density wrt T + * + * Note -> needs augmenting with a T,P implementation. + * + * Verification: + * Agrees with the CRC values (6-10) for up to 4 sig digits. + * + * units = returns density in kg m-3. + * + * (static) + */ + double WaterProps::density_T(double T, double P, int ifunc) { + double Tc = T - 273.15; + const double U1 = 288.9414; + const double U2 = 508929.2; + const double U3 = 68.12963; + const double U4 = -3.9863; + + double tmp1 = Tc + U1; + double tmp4 = Tc + U4; + double t4t4 = tmp4 * tmp4; + double tmp3 = Tc + U3; + double rho = 1000. * (1.0 - tmp1*t4t4/(U2 * tmp3)); + + /* + * Impose an ideal gas lower bound on rho. We need this + * to ensure positivity of rho, even though it is + * grossly unrepresentative. + */ + double rhomin = P / (GasConstant * T); + if (rho < rhomin) { + rho = rhomin; + if (ifunc == 1) { + double drhodT = - rhomin / T; + return drhodT; + } else if (ifunc == 3) { + double drhodP = rhomin / P; + return drhodP; + } else if (ifunc == 2) { + double d2rhodT2 = 2.0 * rhomin / (T * T); + return d2rhodT2; + } + } + + if (ifunc == 1) { + double drhodT = 1000./U2 * ( + - tmp4 * tmp4 / (tmp3) + - tmp1 * 2 * tmp4 / (tmp3) + + tmp1 * t4t4 / (tmp3*tmp3) + ); + return drhodT; + } else if (ifunc == 3) { + return 0.0; + } else if (ifunc == 2) { + double t3t3 = tmp3 * tmp3; + double d2rhodT2 = 1000./U2 * + ((-4.0*tmp4-2.0*tmp1)/tmp3 + + (2.0*t4t4 + 4.0*tmp1*tmp4)/t3t3 + - 2.0*tmp1 * t4t4/(t3t3*tmp3)); + return d2rhodT2; + } + + return rho; + } + + + /** + * Dielectric constant for water: + * Bradley-Pitzer equation for the dielectric constant + * of water as a function of temperature and pressure. + * + * ifunc = 0 value + * ifunc = 1 Temperature deriviative + * ifunc = 2 second temperature derivative + * + * @param T temperature in Kelvin + * @param P Pressure in bar + * + * Range of validity 0 to 350C, 0 to 1 kbar pressure + * + * ifunc = 0 return value + * ifunc = 1 return temperature derivative + * ifunc = 2 return temperature second derivative + * ifunc = 3 return pressure first derivative + * + * Validation: + * Numerical experiments indicate that this function agrees with + * the Archer and Wang data in the CRC p. 6-10 to all 4 significant + * digits shown (0 to 100C). + * + * value at 25C, relEps = 78.38 + * + * (statically defined within the object) + */ + double WaterProps::relEpsilon(double T, double P_pascal, + int ifunc = 0) { + const double U1 = 3.4279E2; + const double U2 = -5.0866E-3; + const double U3 = 9.4690E-7; + const double U4 = -2.0525; + const double U5 = 3.1159E3; + const double U6 = -1.8289E2; + const double U7 = -8.0325E3; + const double U8 = 4.2142E6; + const double U9 = 2.1417; + double T2 = T * T; + + double eps1000 = U1 * exp(U2 * T + U3 * T2); + double C = U4 + U5/(U6 + T); + double B = U7 + U8/T + U9 * T; + + double Pbar = P_pascal * 1.0E-5; + double tmpBpar = B + Pbar; + double tmpB1000 = B + 1000.0; + double ltmp = log(tmpBpar/tmpB1000); + double epsRel = eps1000 + C * ltmp; + + if (ifunc == 1 || ifunc == 2) { + double tmpC = U6 + T; + double dCdT = - U5/(tmpC * tmpC); + + double dBdT = - U8/(T * T) + U9; + + double deps1000dT = eps1000 * (U2 + 2.0 * U3 * T); + + double dltmpdT = (dBdT/tmpBpar - dBdT/tmpB1000); + if (ifunc == 1) { + double depsReldT = deps1000dT + dCdT * ltmp + C * dltmpdT; + return depsReldT; + } + double T3 = T2 * T; + double d2CdT2 = - 2.0 * dCdT / tmpC; + double d2BdT2 = 2.0 * U8 / (T3); + + double d2ltmpdT2 = (d2BdT2*(1.0/tmpBpar - 1.0/tmpB1000) + + dBdT*dBdT*(1.0/(tmpB1000*tmpB1000) - 1.0/(tmpBpar*tmpBpar))); + + double d2eps1000dT2 = (deps1000dT * (U2 + 2.0 * U3 * T) + eps1000 * (2.0 * U3)); + + if (ifunc == 2) { + double d2epsReldT2 = (d2eps1000dT2 + d2CdT2 * ltmp + 2.0 * dCdT * dltmpdT + + C * d2ltmpdT2); + return d2epsReldT2; + } + } + if (ifunc == 3) { + double dltmpdP = 1.0E-5 / tmpBpar; + double depsReldP = C * dltmpdP; + return depsReldP; + } + + return epsRel; + } + + /** + * ADebye calculates the value of A_Debye as a function + * of temperature and pressure according to relations + * that take into account the temperature and pressure + * dependence of the water density and dieletric constant. + * + * A_Debye -> this expression appears on the top of the + * ln actCoeff term in the general Debye-Huckel + * expression + * It depends on temperature. And, therefore, + * most be recalculated whenever T or P changes. + * + * A_Debye = (1/(8 Pi)) sqrt(2 Na dw / 1000) + * (e e/(epsilon R T))^3/2 + * + * Units = sqrt(kg/gmol) ~ sqrt(1/I) + * + * Nominal value = 1.172576 sqrt(kg/gmol) + * based on: + * epsilon/epsilon_0 = 78.54 + * (water at 25C) + * epsilon_0 = 8.854187817E-12 C2 N-1 m-2 + * e = 1.60217653E-19 C + * F = 9.6485309E7 C kmol-1 + * R = 8.314472E3 kg m2 s-2 kmol-1 K-1 + * T = 298.15 K + * B_Debye = 3.28640E9 sqrt(kg/gmol)/m + * Na = 6.0221415E26 + * + * ifunc = 0 return value + * ifunc = 1 return temperature derivative + * ifunc = 2 return temperature second derivative + * ifunc = 3 return pressure first derivative + * + * Verification: + * With the epsRelWater value from the BP relation, + * and the water density from the WaterDens function, + * The A_Debye computed with this function agrees with + * the Pitzer table p. 99 to 4 significant digits at 25C. + * and 20C. (Aphi = ADebye/3) + * + * (statically defined within the object) + */ + double WaterProps::ADebye(double T, double P_input, int ifunc) { + const double e = 1.60217653E-19; + const double epsilon0 = 8.854187817E-12; + const double R = 8.314472E3; + double psat = satPressure(T); + double P; + if (psat > P_input) { + //printf("ADebye WARNING: p_input < psat: %g %g\n", + // P_input, psat); + P = psat; + } else { + P = P_input; + } + double epsRelWater = relEpsilon(T, P, 0); + //printf("releps calc = %g, compare to 78.38\n", epsRelWater); + //double B_Debye = 3.28640E9; + const double Na = 6.0221415E26; + + double epsilon = epsilon0 * epsRelWater; + double dw = density_IAPWS(T, P); + double tmp = sqrt( 2.0 * Na * dw / 1000.); + double tmp2 = e * e * Na / (epsilon * R * T); + double tmp3 = tmp2 * sqrt(tmp2); + double A_Debye = tmp * tmp3 / (8.0 * Pi); + + + /* + * dAdT = - 3/2 Ad/T + 1/2 Ad/dw d(dw)/dT - 3/2 Ad/eps d(eps)/dT + * + * dAdT = - 3/2 Ad/T - 1/2 Ad/Vw d(Vw)/dT - 3/2 Ad/eps d(eps)/dT + */ + if (ifunc == 1 || ifunc == 2) { + double dAdT = - 1.5 * A_Debye / T; + + double depsRelWaterdT = relEpsilon(T, P, 1); + dAdT -= A_Debye * (1.5 * depsRelWaterdT / epsRelWater); + + //int methodD = 1; + //double ddwdT = density_T_new(T, P, 1); + // double contrib1 = A_Debye * (0.5 * ddwdT / dw); + + /* + * calculate d(lnV)/dT _constantP, i.e., the cte + */ + double cte = coeffThermalExp_IAPWS(T, P); + double contrib2 = - A_Debye * (0.5 * cte); + + //dAdT += A_Debye * (0.5 * ddwdT / dw); + dAdT += contrib2; + +#ifdef DEBUG_HKM + //printf("dAdT = %g, contrib1 = %g, contrib2 = %g\n", + // dAdT, contrib1, contrib2); +#endif + + if (ifunc == 1) { + return dAdT; + } + + if (ifunc == 2) { + /* + * Get the second derivative of the dielectric constant wrt T + * -> we will take each of the terms in dAdT and differentiate + * it again. + */ + double d2AdT2 = 1.5 / T * (A_Debye/T - dAdT); + + double d2epsRelWaterdT2 = relEpsilon(T, P, 2); + + //double dT = -0.01; + //double TT = T + dT; + //double depsRelWaterdTdel = relEpsilon(TT, P, 1); + //double d2alt = (depsRelWaterdTdel- depsRelWaterdT ) / dT; + //printf("diff %g %g\n",d2epsRelWaterdT2, d2alt); + // HKM -> checks out, i.e., they are the same. + + d2AdT2 += 1.5 * (- dAdT * depsRelWaterdT / epsRelWater + - A_Debye / epsRelWater * + (d2epsRelWaterdT2 - depsRelWaterdT * depsRelWaterdT / epsRelWater)); + + double deltaT = -0.1; + double Tdel = T + deltaT; + double cte_del = coeffThermalExp_IAPWS(Tdel, P); + double dctedT = (cte_del - cte) / Tdel; + + + //double d2dwdT2 = density_T_new(T, P, 2); + + double contrib3 = 0.5 * ( -(dAdT * cte) -(A_Debye * dctedT)); + d2AdT2 += contrib3; + + return d2AdT2; + } + } + /* + * A_Debye = (1/(8 Pi)) sqrt(2 Na dw / 1000) + * (e e/(epsilon R T))^3/2 + * + * dAdP = + 1/2 Ad/dw d(dw)/dP - 3/2 Ad/eps d(eps)/dP + * + * dAdP = - 1/2 Ad/Vw d(Vw)/dP - 3/2 Ad/eps d(eps)/dP + * + * dAdP = + 1/2 Ad * kappa - 3/2 Ad/eps d(eps)/dP + * + * where kappa = - 1/Vw d(Vw)/dP_T (isothermal compressibility) + */ + if (ifunc == 3) { + + double dAdP = 0.0; + + double depsRelWaterdP = relEpsilon(T, P, 3); + dAdP -= A_Debye * (1.5 * depsRelWaterdP / epsRelWater); + + double kappa = isothermalCompressibility_IAPWS(T, P); + + //double ddwdP = density_T_new(T, P, 3); + dAdP += A_Debye * (0.5 * kappa); + + return dAdP; + } + + return A_Debye; + } + + double WaterProps::satPressure(double T) { + double pres = m_waterIAPWS->psat(T); + return pres; + } + + + double WaterProps::density_IAPWS(double temp, double press) { + + double dens; + dens = m_waterIAPWS->density(temp, press, WATER_LIQUID); + return dens; + } + + double WaterProps::coeffThermalExp_IAPWS(double temp, double press) { + + double cte; + cte = m_waterIAPWS->coeffThermExp(temp, press); + return cte; + } + + double WaterProps::isothermalCompressibility_IAPWS(double temp, double press) { + + double kappa; + kappa = m_waterIAPWS->isothermalCompressibility(temp, press); + return kappa; + } + + +} diff --git a/Cantera/src/thermo/WaterProps.h b/Cantera/src/thermo/WaterProps.h new file mode 100644 index 000000000..87353c05e --- /dev/null +++ b/Cantera/src/thermo/WaterProps.h @@ -0,0 +1,144 @@ +/** + * @file WaterProps.h + * + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ + +#ifndef CT_WATERPROPS_H +#define CT_WATERPROPS_H + + +#include "ct_defs.h" +class WaterPropsIAPWS; +namespace Cantera { + + class WaterPDSS; + + + /** + * Definition of the WaterProps class. This class is used to + * house several approximation routines for properties of water + * Most if not all of the member functions are static. + */ + class WaterProps { + + public: + + WaterProps(); + + WaterProps(WaterPDSS *wptr); + + WaterProps(const WaterProps &b); + + virtual ~WaterProps(); + + WaterProps& operator=(const WaterProps&b); + + /* + * Simple calculation of water density at atmospheric pressure. + * Valid up to boiling point. + * + * ifunc = 0 Returns the density in kg/m^3 + * ifunc = 1 returns the derivative of the density wrt T. + * ifunc = 2 returns the derivative of the density wrt P. + * ifunc = 3 returns the 2nd derivative of the density wrt T + * + * Note -> needs augmenting with a T,P implementation. + * + * Verification: + * Agrees with the CRC values (6-10) for up to 4 sig digits. + * + * units = returns density in kg m-3. + */ + static double density_T(double T, double P, int ifunc); + + /** + * Dielectric constant for water: + * Bradley-Pitzer equation for the dielectric constant + * of water as a function of temperature and pressure. + * + * ifunc = 0 value + * ifunc = 1 Temperature deriviative + * ifunc = 2 second temperature derivative + * + * @param T temperature in Kelvin + * @param P Pressure in bar + * + * Range of validity 0 to 350C, 0 to 1 kbar pressure + * + * ifunc = 0 return value + * ifunc = 1 return temperature derivative + * + * Validation: + * Numerical experiments indicate that this function agrees with + * the Archer and Wang data in the CRC p. 6-10 to all 4 significant + * digits shown (0 to 100C). + * + * value at 25C, relEps = 78.38 + */ + static double relEpsilon(double T, double P_pascal, int ifunc); + + /** + * ADebye calculates the value of A_Debye as a function + * of temperature and pressure according to relations + * that take into account the temperature and pressure + * dependence of the water density and dieletric constant. + * + * A_Debye -> this expression appears on the top of the + * ln actCoeff term in the general Debye-Huckel + * expression + * It depends on temperature. And, therefore, + * most be recalculated whenever T or P changes. + * + * A_Debye = (1/8Pi) sqrt(2Na dw/1000) + * (e e/(epsilon RT)^3/2 + * + * Units = sqrt(kg/gmol) + * + * Nominal value = 1.172576 sqrt(kg/gmol) + * based on: + * epsilon/epsilon_0 = 78.54 + * (water at 25C) + * epsilon_0 = 8.854187817E-12 C2 N-1 m-2 + * e = 1.60217653E-19 C + * F = 9.6485309E7 C kmol-1 + * R = 8.314472E3 kg m2 s-2 kmol-1 K-1 + * T = 298.15 K + * B_Debye = 3.28640E9 sqrt(kg/gmol)/m + * Na = 6.0221415E26 + * + * Verification: + * With the epsRelWater value from the BP relation, + * and the water density from the WaterDens function, + * The A_Debye computed with this function agrees with + * the Pitzer table p. 99 to 4 significant digits at 25C. + * and 20C. (Aphi = ADebye/3) + */ + double ADebye(double T, double P, int ifunc); + + double satPressure(double T); + + + double density_IAPWS(double T, double P); + double coeffThermalExp_IAPWS(double T, double P); + double isothermalCompressibility_IAPWS(double T, double P); + + protected: + + + WaterPropsIAPWS *m_waterIAPWS; + bool m_own_sub; + }; + + +} + + +#endif diff --git a/Cantera/src/thermo/WaterPropsIAPWS.cpp b/Cantera/src/thermo/WaterPropsIAPWS.cpp new file mode 100644 index 000000000..35c685eaf --- /dev/null +++ b/Cantera/src/thermo/WaterPropsIAPWS.cpp @@ -0,0 +1,587 @@ +/* + * @file WaterPropsIAPWS + * + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ + +#include "WaterPropsIAPWS.h" +#include +#include +#include +/* + * Critical Point values in mks units + */ +static const double T_c = 647.096; // Kelvin +static const double P_c = 22.064E6; // Pascals +static const double Rho_c = 322.; // kg m-3 +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 + */ +//static const double Rgas = 8.314472E3; // Joules kmol-1 K-1 +static const double Rgas = 8.314371E3; // Joules kmol-1 K-1 + + +WaterPropsIAPWS:: WaterPropsIAPWS() : + m_phi(0), + tau(-1.0), + delta(-1.0), + iState(-30000) +{ + m_phi = new WaterPropsIAPWSphi(); +} + +WaterPropsIAPWS::WaterPropsIAPWS(const WaterPropsIAPWS &b) : + m_phi(0), + tau(b.tau), + delta(b.delta), + iState(b.iState) +{ + m_phi = new WaterPropsIAPWSphi(); + m_phi->tdpolycalc(tau, delta); +} + +WaterPropsIAPWS & WaterPropsIAPWS::operator=(const WaterPropsIAPWS &b) { + if (this == &b) return *this; + tau = b.tau; + delta = b.delta; + iState = b.iState; + m_phi->tdpolycalc(tau, delta); + return *this; +} + +WaterPropsIAPWS::~WaterPropsIAPWS() { + delete (m_phi); + m_phi = 0; +} + + +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); +} + +/** + * Calculate the Helmholtz free energy in mks units of + * J kmol-1 K-1. + */ +double WaterPropsIAPWS::helmholtzFE(double temperature, double rho) { + setState(temperature, rho); + double retn = helmholtzFE_RT(); + double RT = Rgas * temperature; + return (retn * RT); +} +double WaterPropsIAPWS::helmholtzFE() const{ + double retn = helmholtzFE_RT(); + double temperature = T_c/tau; + double RT = Rgas * temperature; + return (retn * RT); +} + + +/** + * Calculate the pressure (Pascals), given the temperature and density + * Temperature: kelvin + * rho: density in kg m-3 + */ +double WaterPropsIAPWS::pressure(double temperature, double rho) { + calcDim(temperature, rho); + double retn = pressure_rhoRT(); + return (retn * rho * Rgas * temperature); +} +double WaterPropsIAPWS::pressure() const{ + double retn = pressure_rhoRT(); + double rho = delta * Rho_c; + double temperature = T_c / tau; + return (retn * rho * Rgas * temperature); +} + +/** + * Calculates the pressure in dimensionless form + * p/(rhoRT) at the currently stored tau and delta values + */ +double WaterPropsIAPWS::pressure_rhoRT() const { + double retn = m_phi->pressure_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 + * 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 + * + * If a problem is encountered, a negative 1 is returned. + */ +double WaterPropsIAPWS:: +density(double temperature, double pressure, int phase, double rhoguess) { + + double deltaGuess = 0.0; + if (rhoguess == -1.0) { + if (phase != -1) { + if (temperature > T_c) { + rhoguess = pressure * M_water / (Rgas * temperature); + } else { + if (phase != WATER_LIQUID) { + rhoguess = pressure * M_water / (Rgas * temperature); + } else { + /* + * Provide a guess about the liquid density + */ + rhoguess = 1000.; + } + } + } else { + /* + * Assume the Gas phase initial guess, if nothing is + * specified to the routine + */ + rhoguess = pressure * M_water / (Rgas * temperature); + } + + } + double p_red = pressure * M_water / (Rgas * temperature * Rho_c); + deltaGuess = rhoguess / Rho_c; + calcDim(temperature, rhoguess); + double delta_retn = m_phi->dfind(p_red, tau, deltaGuess); + double density_retn; + if (delta_retn >0.0) { + delta = delta_retn; + + /* + * Dimensionalize the density before returning + */ + density_retn = delta_retn * Rho_c; + /* + * Determine the internal state + */ + if (temperature > T_c) { + iState = WATER_SUPERCRIT; + } else { + if (delta_retn < 1.0) { + iState = WATER_GAS; + } else { + iState = WATER_LIQUID; + } + } + + } else { + density_retn = -1.0; + } + return density_retn; +} + +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. + * + * Input + * temperature (kelvin) + * + * return: + * psat (Pascals) + */ +double WaterPropsIAPWS::psat_est(double temperature) { + + static const double A[8] = { + -7.8889166E0, + 2.5514255E0, + -6.716169E0, + 33.2239495E0, + -105.38479E0, + 174.35319E0, + -148.39348E0, + 48.631602E0 + }; + double ps; + if (temperature < 314.) { + double pl = 6.3573118E0 - 8858.843E0 / temperature + + 607.56335E0 * pow(temperature, -0.6); + ps = 0.1 * exp(pl); + } else { + double v = temperature / 647.25; + double w = fabs(1.0-v); + double b = 0.0; + for (int i = 0; i < 8; i++) { + double z = i + 1; + b += A[i] * pow(w, ((z+1.0)/2.0)); + } + double q = b / v; + ps = 22.093*exp(q); + } + /* + * Original correlation was in cgs. Convert to mks + */ + ps *= 1.0E6; + return ps; +} + +/** + * Returns the coefficient of thermal expansion as a function + * of temperature and pressure. + * alpha = d (ln V) / dT at constant P. + * + * Currently this function is calculated using a differencing scheme. + */ +double WaterPropsIAPWS::coeffThermExp(double temperature, double pressure) { + + double deltaT = 0.01; + double psat_at=0.0; + double rhoguess = -1; + int phase = -1; + if (temperature > T_c) { + rhoguess = pressure * M_water / (Rgas * temperature); + } else { + psat_at = psat(temperature); + if (pressure >= psat_at) { + phase = WATER_LIQUID; + deltaT = -0.01; + } else + phase = WATER_GAS; + } + + double dens_base = density(temperature, pressure, phase, rhoguess); + + if (dens_base == -1.0) { + printf("problems\n"); + exit(-1); + } + double temp_del = temperature + deltaT; + + double dens_del = density(temp_del, pressure, phase, dens_base); + + double Vavg = 0.5 * (1./dens_del + 1./dens_base); + double retn = 1.0 / Vavg * (1./dens_del - 1.0/dens_base)/deltaT; + return retn; +} + +/** + * Returns the coefficient of isothermal compressibility + * of temperature and pressure. + * kappa = - d (ln V) / dP at constant T. + * + * Currently this function is calculated using an inaccurate + * one-sided differencing scheme. + */ +double WaterPropsIAPWS:: +isothermalCompressibility(double temperature, double pressure) { + /* + * Difference amount is large, because we are solving for + * density underneath + */ + double deltaP = -0.001 * pressure; + double psat_at=0.0; + double rhoguess = -1; + int phase = -1; + if (temperature > T_c) { + rhoguess = pressure * M_water / (Rgas * temperature); + deltaP = +0.0001 * pressure; + phase = WATER_SUPERCRIT; + } else { + psat_at = psat(temperature); + if (pressure >= psat_at) { + phase = WATER_LIQUID; + deltaP = +0.0001 * pressure; + } else + phase = WATER_GAS; + deltaP = -0.0001 * pressure; + } + double dens_base = density(temperature, pressure, phase, rhoguess); + if (dens_base == -1.0) { + printf("problems\n"); + exit(-1); + } + double pres_del = pressure + deltaP; + double dens_del = density(temperature, pres_del, phase, dens_base); + double Vavg = 0.5 * (1./dens_del + 1./dens_base); + double retn = -1.0 / Vavg * (1./dens_del - 1.0/dens_base)/deltaP; + 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(); + return (gRT * Rgas * temperature); +} +double WaterPropsIAPWS:: +Gibbs() const { + double gRT = 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. + */ +void WaterPropsIAPWS:: +corr(double temperature, double pressure, double &densLiq, + double &densGas, double &delGRT) { + + densLiq = density(temperature, pressure, WATER_LIQUID, densLiq); + if (densLiq <= 0.0) { + printf("error liq\n"); + exit(-1); + } + setState(temperature, densLiq); + double gibbsLiqRT = Gibbs_RT(); + + densGas = density(temperature, pressure, WATER_GAS, densGas); + if (densGas <= 0.0) { + printf("error gas\n"); + exit(-1); + } + setState(temperature, densGas); + double gibbsGasRT = Gibbs_RT(); + + delGRT = gibbsLiqRT - gibbsGasRT; +} + +void WaterPropsIAPWS:: +corr1(double temperature, double pressure, double &densLiq, + double &densGas, double &pcorr) { + + densLiq = density(temperature, pressure, WATER_LIQUID, densLiq); + setState(temperature, densLiq); + double prL = m_phi->phiR(); + + densGas = density(temperature, pressure, WATER_GAS, densGas); + setState(temperature, densGas); + double prG = m_phi->phiR(); + + double rhs = (prL - prG) + log(densLiq/densGas); + rhs /= (1.0/densGas - 1.0/densLiq); + + pcorr = rhs * Rgas * temperature / M_water; +} + +/** + * Calculate the saturation pressure given the temperature. + * p : Pascals : Newtons/m**2 + */ +static int method = 1; +double WaterPropsIAPWS:: +psat(double temperature) { + double densLiq = -1.0, densGas = -1.0, delGRT = 0.0; + double dp, pcorr; + double p = psat_est(temperature); + bool conv = false; + for (int i = 0; i < 30; i++) { + if (method == 1) { + corr(temperature, p, densLiq, densGas, delGRT); + double delV = M_water * (1.0/densLiq - 1.0/densGas); + dp = - delGRT * Rgas * temperature / delV; + } else { + corr1(temperature, p, densLiq, densGas, pcorr); + dp = pcorr - p; + } + p += dp; + + if ((method == 1) && delGRT < 1.0E-8) { + conv = true; + break; + } else { + if (fabs(dp/p) < 1.0E-9) { + conv = true; + break; + } + } + } + return p; +} + +/** + * Sets the internal state of the object to the + * specified temperature and density. + */ +void WaterPropsIAPWS:: +setState(double temperature, double rho) { + calcDim(temperature, 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(); + return (hRT * Rgas * temperature); +} +double WaterPropsIAPWS:: +enthalpy() const { + double temperature = T_c/tau; + double hRT = 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(); + return (uRT * Rgas * temperature); +} +double WaterPropsIAPWS:: +intEnergy() const{ + double temperature = T_c / tau; + double uRT = 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(); + return (sR * Rgas); +} + +/** + * Calculate the enthalpy in mks units of + * J kmol-1 K-1. + */ +double WaterPropsIAPWS:: +entropy() const { + double sR = 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(); + 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(); + return (cpR * Rgas); +} +double WaterPropsIAPWS:: +cp() const { + double cpR = cp_R(); + return (cpR * Rgas); +} + +double WaterPropsIAPWS:: +molarVolume(double temperature, double rho) { + setState(temperature, rho); + return (M_water / rho); +} + +double WaterPropsIAPWS:: +molarVolume() const { + double rho = delta * Rho_c; + return (M_water / rho); +} diff --git a/Cantera/src/thermo/WaterPropsIAPWS.h b/Cantera/src/thermo/WaterPropsIAPWS.h new file mode 100644 index 000000000..fe63aff77 --- /dev/null +++ b/Cantera/src/thermo/WaterPropsIAPWS.h @@ -0,0 +1,203 @@ +/** + * @file WaterPropsIAPWS.h + * + */ +/* + * Copywrite (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ + +#ifndef WATERPROPSIAPWS_H +#define WATERPROPSIAPWS_H + +#include "WaterPropsIAPWSphi.h" + +/* + * These constants are defined and used in the interphase + * to describe desired phases. + */ +#define WATER_GAS 0 +#define WATER_LIQUID 1 +#define WATER_SUPERCRIT 2 + +/** + * Class for calculating the properties of water. + * + * + * Note, the base thermodynamic state for this class is the one + * 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 + * + */ +class WaterPropsIAPWS { +public: + WaterPropsIAPWS(); + WaterPropsIAPWS(const WaterPropsIAPWS &b); + WaterPropsIAPWS & operator=(const WaterPropsIAPWS &b); + ~WaterPropsIAPWS(); + + + void setState(double temperature, double rho); + + /** + * Calculate the Helmholtz free energy in mks units of + * J kmol-1 K-1. + */ + double helmholtzFE(double temperature, double rho); + double helmholtzFE() const; + + /** + * Calculate the Gibbs free energy in mks units of + * J kmol-1 K-1. + */ + double Gibbs(double temperature, double rho); + double Gibbs() const; + + /** + * Calculate the enthalpy in mks units of + * J kmol-1 + */ + double enthalpy(double temperature, double rho); + double enthalpy() const; + + /** + * Calculate the internal energy in mks units of + * J kmol-1 + */ + double intEnergy(double temperature, double rho); + double intEnergy() const; + + /** + * Calculate the entropy in mks units of + * J kmol-1 K-1 + */ + double entropy(double temperature, double rho); + double entropy() const; + + /** + * Calculate the constant volume heat capacity + * in mks units of J kmol-1 K-1 + */ + double cv(double temperature, double rho); + double cv() const; + + /** + * Calculate the constant pressure heat capacity + * in mks units of J kmol-1 K-1 + */ + double cp(double temperature, double rho); + double cp() const; + + 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 + */ + double pressure(double temperature, double rho); + 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. + * + * 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 + */ + double density(double temperature, double pressure, + int phase = -1, double rhoguess = -1.0); + 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) + */ + 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. + * + * + */ + double coeffThermExp(double temperature, double 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 + */ + double isothermalCompressibility(double temperature, double pressure); + + /** + * Utility routine in the calculation of the saturation pressure + */ + void corr(double temperature, double pressure, double &densLiq, + double &densGas, double &delGRT); + 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 + */ + double psat(double temperature); + + double Tcrit() { return 647.096;} + double Pcrit() { return 22.064E6;} + + double Rhocrit() { return 322.;} + + +private: + /** + * Calculate the dimensionless temp and rho and store internally. + */ + 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. + */ + double helmholtzFE_RT() const; + double Gibbs_RT() const; + double enthalpy_RT() const; + double intEnergy_RT() const; + double entropy_R() const; + double cv_R() const; + double cp_R() const; + double pressure_rhoRT() const; + +protected: + WaterPropsIAPWSphi *m_phi; + + double tau; + double delta; + int iState; +}; +#endif + diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.cpp b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp new file mode 100644 index 000000000..ccf1e1bba --- /dev/null +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp @@ -0,0 +1,1244 @@ +/** + * @file WaterPropsIAPWSphi.cpp + * + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ + +#include "WaterPropsIAPWSphi.h" + +#include +#include + +/* + * Critical Point values in mks units: Note, these aren't used in this + * routine, except for internal checks. All calculations here are done + * in dimensionless units. + */ +static const double T_c = 647.096; // Kelvin +static const double P_c = 22.064E6; // Pascals +static const double Rho_c = 322.; // kg m-3 +static const double M_water = 18.015268; // kg kmol-1 + +/* + * The added constants were calculated so that u = s = 0 + * for liquid at the triple point. These where determined + * by the program testPress. I'm not quite satisfied with + * the result, but will let it stand for the moment. + * H didn't turn out to be .611872 J/kg, but .611782 J/kg. + * There may be a slight error here somehow. + */ +static const double ni0[9] = { + 0.0, + -8.32044648201 - 0.000000001739715, + 6.6832105268 + 0.000000000793232, + 3.00632, + 0.012436, + 0.97315, + 1.27950, + 0.96956, + 0.24873 +}; + +static const double gammi0[9] = { + 0.0, + 0.0, + 0.0, + 0.0, + 1.28728967, + 3.53734222, + 7.74073708, + 9.24437796, + 27.5075105 +}; + +static const int ciR[56] = { + 0, // 0 + 0, // 1 + 0, + 0, + 0, + 0, // 5 + 0, + 0, + 1, + 1, + 1, // 10 + 1, + 1, + 1, + 1, + 1, // 15 + 1, + 1, + 1, + 1, + 1, // 20 + 1, + 1, + 2, + 2, + 2, // 25 + 2, + 2, + 2, + 2, + 2, // 30 + 2, + 2, + 2, + 2, + 2, // 35 + 2, + 2, + 2, + 2, + 2, // 40 + 2, + 2, + 3, + 3, + 3, // 45 + 3, + 4, + 6, + 6, + 6, // 50 + 6, + 0, + 0, + 0, + 0 // 55 +}; + +static const int diR[55] = { + 0, // 0 + 1, // 1 + 1, + 1, + 2, + 2, // 5 + 3, + 4, + 1, + 1, + 1, // 10 + 2, + 2, + 3, + 4, + 4, // 15 + 5, + 7, + 9, + 10, + 11, // 20 + 13, + 15, + 1, + 2, + 2, // 25 + 2, + 3, + 4, + 4, + 4, // 30 + 5, + 6, + 6, + 7, + 9, // 35 + 9, + 9, + 9, + 9, + 10, // 40 + 10, + 12, + 3, + 4, + 4, // 45 + 5, + 14, + 3, + 6, + 6, // 50 + 6, + 3, + 3, + 3 // 54 +}; + +static const int tiR[55] = { + 0, // 0 + 0, // 1 + 0, + 0, + 0, + 0, // 5 + 0, + 0, + 4, // 8 + 6, + 12, // 10 + 1, + 5, + 4, + 2, + 13, // 15 + 9, + 3, + 4, + 11, + 4, // 20 + 13, + 1, + 7, + 1, + 9, // 25 + 10, + 10, + 3, + 7, + 10, // 30 + 10, + 6, + 10, + 10, + 1, // 35 + 2, + 3, + 4, + 8, + 6, // 40 + 9, + 8, + 16, + 22, + 23, // 45 + 23, + 10, + 50, + 44, + 46, // 50 + 50, + 0, + 1, + 4 // 54 +}; + +static const double ni[57] = { + +0.0, + +0.12533547935523E-1, // 1 + +0.78957634722828E1, // 2 + -0.87803203303561E1, // 3 + +0.31802509345418E0, // 4 + -0.26145533859358E0, // 5 + -0.78199751687981E-2, // 6 + +0.88089493102134E-2, // 7 + -0.66856572307965E0, // 8 + +0.20433810950965, // 9 + -0.66212605039687E-4, // 10 + -0.19232721156002E0, // 11 + -0.25709043003438E0, // 12 + +0.16074868486251E0, // 13 + -0.40092828925807E-1, // 14 + +0.39343422603254E-6, // 15 + -0.75941377088144E-5, // 16 + +0.56250979351888E-3, // 17 + -0.15608652257135E-4, // 18 + +0.11537996422951E-8, // 19 + +0.36582165144204E-6, // 20 + -0.13251180074668E-11,// 21 + -0.62639586912454E-9, // 22 + -0.10793600908932E0, // 23 + +0.17611491008752E-1, // 24 + +0.22132295167546E0, // 25 + -0.40247669763528E0, // 26 + +0.58083399985759E0, // 27 + +0.49969146990806E-2, // 28 + -0.31358700712549E-1, // 29 + -0.74315929710341E0, // 30 + +0.47807329915480E0, // 31 + +0.20527940895948E-1, // 32 + -0.13636435110343E0, // 33 + +0.14180634400617E-1, // 34 + +0.83326504880713E-2, // 35 + -0.29052336009585E-1, // 36 + +0.38615085574206E-1, // 37 + -0.20393486513704E-1, // 38 + -0.16554050063734E-2, // 39 + +0.19955571979541E-2, // 40 + +0.15870308324157E-3, // 41 + -0.16388568342530E-4, // 42 + +0.43613615723811E-1, // 43 + +0.34994005463765E-1, // 44 + -0.76788197844621E-1, // 45 + +0.22446277332006E-1, // 46 + -0.62689710414685E-4, // 47 + -0.55711118565645E-9, // 48 + -0.19905718354408E0, // 49 + +0.31777497330738E0, // 50 + -0.11841182425981E0, // 51 + -0.31306260323435E2, // 52 + +0.31546140237781E2, // 53 + -0.25213154341695E4, // 54 + -0.14874640856724E0, // 55 + +0.31806110878444E0 // 56 +}; + + +static const double alphai[3] = { + +20., + +20., + +20. +}; + +static const double betai[3] = { + +150., + +150., + +250. +}; + +static const double gammai[3] = { + +1.21, + +1.21, + +1.25 +}; + +static const double epsi[3] = { + +1.0, + +1.0, + +1.0 +}; + +static const double ai[2] = { + +3.5, + +3.5 +}; + +static const double bi[2] = { + +0.85, + +0.95 +}; + +static const double Bi[2] = { + +0.2, + +0.2 +}; + +static const double Ci[2] = { + +28.0, + +32.0 +}; + +static const double Di[2] = { + +700., + +800. +}; + +static const double Ai[2] = { + +0.32, + +0.32 +}; + +static const double Bbetai[2] = { + +0.3, + +0.3 +}; + +/** + * Constructor for the object. + */ +WaterPropsIAPWSphi::WaterPropsIAPWSphi() : + TAUsave(-1.0), + TAUsqrt(-1.0), + DELTAsave(-1.0) +{ +} + +/** + * intCheck() calculates all of the functions at a one point and + * prints out the result. It's used for conducting the internal + * check. + */ +void WaterPropsIAPWSphi::intCheck(double tau, double delta) { + tdpolycalc(tau, delta); + double nau = phi0(); + double res = phiR(); + double res_d = phiR_d(); + double nau_d = phi0_d(); + double res_dd = phiR_dd(); + double nau_dd = phi0_dd(); + double res_t = phiR_t(); + double nau_t = phi0_t(); + double res_tt = phiR_tt(); + double nau_tt = phi0_tt(); + double res_dt = phiR_dt(); + double nau_dt = phi0_dt(); + + printf("nau = %20.13e\t\tres = %20.13e\n", nau, res); + printf("nau_d = %20.13e\t\tres_d = %20.13e\n", nau_d, res_d); + printf("nau_dd = %20.13e\t\tres_dd = %20.13e\n", nau_dd, res_dd); + printf("nau_t = %20.13e\t\tres_t = %20.13e\n", nau_t, res_t); + printf("nau_tt = %20.13e\t\tres_tt = %20.13e\n", nau_tt, res_tt); + printf("nau_dt = %20.13e\t\tres_dt = %20.13e\n", nau_dt, res_dt); +} + +void WaterPropsIAPWSphi::check1() { + double T = 500.; + double rho = 838.025; + double tau = T_c/T; + double delta = rho / Rho_c; + printf(" T = 500 K, rho = 838.025 kg m-3\n"); + intCheck(tau, delta); +} + +void WaterPropsIAPWSphi::check2() { + double T = 647; + double rho = 358.0; + double tau = T_c/T; + double delta = rho / Rho_c; + printf(" T = 647 K, rho = 358.0 kg m-3\n"); + intCheck(tau, delta); +} + +/* + * Calculate the polynomials in tau and delta, and store them in static + * storage. + */ +void WaterPropsIAPWSphi::tdpolycalc(double tau, double delta) { + if ((tau != TAUsave) || 1) { + TAUsave = tau; + TAUsqrt = sqrt(tau); + TAUp[0] = 1.0; + for (int i = 1; i < 51; i++) { + TAUp[i] = TAUp[i-1] * tau; + } + } + if ((delta != DELTAsave) || 1) { + DELTAsave = delta; + DELTAp[0] = 1.0; + for (int i = 1; i <= 15; i++) { + DELTAp[i] = DELTAp[i-1] * delta; + } + } +} + +/* + * Calculate Eqn. 6.5 for phi0, the ideal gas part of the + * dimensionless Helmholtz free energy. + */ +double WaterPropsIAPWSphi::phi0() const { + double tau = TAUsave; + double delta = DELTAsave; + double retn = log(delta) + ni0[1] + ni0[2]*tau + ni0[3]*log(tau); + + retn += ni0[4] * log(1.0 - exp(-gammi0[4]*tau)); + retn += ni0[5] * log(1.0 - exp(-gammi0[5]*tau)); + retn += ni0[6] * log(1.0 - exp(-gammi0[6]*tau)); + retn += ni0[7] * log(1.0 - exp(-gammi0[7]*tau)); + retn += ni0[8] * log(1.0 - exp(-gammi0[8]*tau)); + return retn; +} + +/* + * Calculate Eqn. 6.6 for phiR, the residual part of the + * dimensionless Helmholtz free energy. + * + * tau = dimensionless temperature + * delta = dimensionless pressure + */ +double WaterPropsIAPWSphi::phiR() const { + double tau = TAUsave; + double delta = DELTAsave; + int i, j; + + /* + * Write out the first seven polynomials in the expression + */ + double T375 = pow(tau, 0.375); + double val = (ni[1] * delta / TAUsqrt + + ni[2] * delta * TAUsqrt * T375 + + ni[3] * delta * tau + + ni[4] * DELTAp[2] * TAUsqrt + + ni[5] * DELTAp[2] * T375 * T375 + + ni[6] * DELTAp[3] * T375 + + ni[7] * DELTAp[4] * tau); + /* + * Next, do polynomial contributions 8 to 51 + */ + for (i = 8; i <= 51; i++) { + val += (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * exp(-DELTAp[ciR[i]])); + } + + /* + * Next do contributions 52 to 54 + */ + for (j = 0; j < 3; j++) { + i = 52 + j; + double dtmp = delta - epsi[j]; + double ttmp = tau - gammai[j]; + val += (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * + exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp)); + } + + /* + * Next do contributions 55 and 56 + */ + for (j = 0; j < 2; j++) { + i = 55 + j; + double deltam1 = delta - 1.0; + double dtmp2 = deltam1 * deltam1; + double atmp = 0.5 / Bbetai[j]; + double theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp); + double triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]); + double ttmp = tau - 1.0; + + double triagtmp = pow(triag, bi[j]); + + double phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp); + val += (ni[i] * triagtmp * delta * phi); + } + + return val; +} + +/* + * Calculate the Phi function, which is basically the helmholtz free energy + * Eqn. (6.4) + */ +double WaterPropsIAPWSphi::phi(double tau, double delta) { + tdpolycalc(tau, delta); + double nau = phi0(); + double res = phiR(); + double retn = nau + res; + return retn; +} + + +/* + * Calculate d_phiR_d(delta), the first derivative of phiR + * wrt delta + * + * tau = dimensionless temperature + * delta = dimensionless pressure + */ +double WaterPropsIAPWSphi::phiR_d() const { + double tau = TAUsave; + double delta = DELTAsave; + int i, j; + + /* + * Write out the first seven polynomials in the expression + */ + double T375 = pow(tau, 0.375); + double val = (ni[1] / TAUsqrt + + ni[2] * TAUsqrt * T375 + + ni[3] * tau + + ni[4] * 2.0 * delta * TAUsqrt + + ni[5] * 2.0 * delta * T375 * T375 + + ni[6] * 3.0 * DELTAp[2] * T375 + + ni[7] * 4.0 * DELTAp[3] * tau); + /* + * Next, do polynomial contributions 8 to 51 + */ + for (i = 8; i <= 51; i++) { + val += ( (ni[i] * exp(-DELTAp[ciR[i]]) * DELTAp[diR[i] - 1] * + TAUp[tiR[i]]) * (diR[i] - ciR[i]* DELTAp[ciR[i]])); + } + + /* + * Next do contributions 52 to 54 + */ + for (j = 0; j < 3; j++) { + i = 52 + j; + double dtmp = delta - epsi[j]; + double ttmp = tau - gammai[j]; + double tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * + exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp)); + val += tmp * (diR[i]/delta - 2.0 * alphai[j] * dtmp); + } + + /* + * Next do contributions 55 and 56 + */ + for (j = 0; j < 2; j++) { + i = 55 + j; + double deltam1 = delta - 1.0; + double dtmp2 = deltam1 * deltam1; + double atmp = 0.5 / Bbetai[j]; + double theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp); + double triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]); + double ttmp = tau - 1.0; + + double triagtmp = pow(triag, bi[j]); + double triagtmpm1 = pow(triag, bi[j]-1.0); + double atmpM1 = atmp - 1.0; + double ptmp = pow(dtmp2,atmpM1); + double p2tmp = pow(dtmp2, ai[j]-1.0); + double dtriagddelta = + deltam1 *(Ai[j] * theta * 2.0 / Bbetai[j] * ptmp + + 2.0*Bi[j]*ai[j]*p2tmp); + + double phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp); + double dphiddelta = -2.0*Ci[j]*deltam1*phi; + double dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta; + + double tmp = ni[i] * (triagtmp * (phi + delta*dphiddelta) + + dtriagtmpddelta * delta * phi); + val += tmp; + } + + return val; +} + +/* + * Calculate d_phi0_d(delta), the first derivative of phi0 + * wrt delta + * + * tau = dimensionless temperature + * delta = dimensionless pressure + */ +double WaterPropsIAPWSphi::phi0_d() const { + double delta = DELTAsave; + return (1.0/delta); +} + +/* + * Calculate the dPhidDelta function, which is basically the derivative + * of helmholtz free energy wrt delta + * Eqn. (6.4) + */ +double WaterPropsIAPWSphi::phi_d(double tau, double delta) { + tdpolycalc(tau, delta); + double nau = phi0_d(); + double res = phiR_d(); + double retn = nau + res; + return retn; +} + +/** + * Calculate the dimensionless pressure at tau and delta; + * + * p/(rhoRT) = delta * phi_d() + * + * note: this is done so much, we have a seperate routine. + */ +double WaterPropsIAPWSphi::pressure_rhoRT(double tau, double delta) { + tdpolycalc(tau, delta); + double res = phiR_d(); + double retn = 1.0 + delta * res; + return retn; +} + +/* + * Calculate d_phiR_d(delta), the second derivative of phiR + * wrt delta + * + * tau = dimensionless temperature + * delta = dimensionless pressure + */ +double WaterPropsIAPWSphi::phiR_dd() const { + double tau = TAUsave; + double delta = DELTAsave; + int i, j; + double atmp; + + /* + * Write out the first seven polynomials in the expression + */ + double T375 = pow(tau, 0.375); + double val = (ni[4] * 2.0 * TAUsqrt + + ni[5] * 2.0 * T375 * T375 + + ni[6] * 6.0 * delta * T375 + + ni[7] * 12.0 * DELTAp[2] * tau); + /* + * Next, do polynomial contributions 8 to 51 + */ + for (i = 8; i <= 51; i++) { + double dtmp = DELTAp[ciR[i]]; + double tmp = ni[i] * exp(-dtmp) * TAUp[tiR[i]]; + if (diR[i] == 1) { + atmp = 1.0/delta; + } else { + atmp = DELTAp[diR[i] - 2]; + } + tmp *= atmp *((diR[i] - ciR[i]*dtmp)*(diR[i]-1.0-ciR[i]*dtmp) - + ciR[i]*ciR[i]*dtmp); + val += tmp; + } + + /* + * Next do contributions 52 to 54 + */ + for (j = 0; j < 3; j++) { + i = 52 + j; + double dtmp = delta - epsi[j]; + double ttmp = tau - gammai[j]; + double tmp = (ni[i] * TAUp[tiR[i]] * + exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp)); + double deltmp = DELTAp[diR[i]]; + double deltmpM1 = deltmp/delta; + double deltmpM2 = deltmpM1 / delta; + double d2tmp = dtmp * dtmp; + + val += tmp * (-2.0*alphai[j]*deltmp + + 4.0 * alphai[j] * alphai[j] * deltmp * d2tmp - + 4.0 * diR[i] * alphai[j] * deltmpM1 * dtmp + + diR[i] * (diR[i] - 1.0) * deltmpM2); + } + + /* + * Next do contributions 55 and 56 + */ + for (j = 0; j < 2; j++) { + i = 55 + j; + double deltam1 = delta - 1.0; + double dtmp2 = deltam1 * deltam1; + atmp = 0.5 / Bbetai[j]; + double theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp); + double triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]); + double ttmp = tau - 1.0; + + double triagtmp = pow(triag, bi[j]); + double triagtmpm1 = pow(triag, bi[j]-1.0); + double atmpM1 = atmp - 1.0; + double ptmp = pow(dtmp2,atmpM1); + double p2tmp = pow(dtmp2, ai[j]-1.0); + double dtriagddelta = + deltam1 *(Ai[j] * theta * 2.0 / Bbetai[j] * ptmp + + 2.0*Bi[j]*ai[j]*p2tmp); + + double phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp); + double dphiddelta = -2.0*Ci[j]*deltam1*phi; + double dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta; + + + double d2phiddelta2 = 2.0 * Ci[j] * phi * (2.0*Ci[j]*dtmp2 - 1.0); + + double pptmp = ptmp / dtmp2; + double d2triagddelta2 = dtriagddelta / deltam1; + d2triagddelta2 += + dtmp2 *(4.0*Bi[j]*ai[j]*(ai[j]-1.0)*pow(dtmp2,ai[j]-2.0) + + 2.0*Ai[j]*Ai[j]/(Bbetai[j]*Bbetai[j])*ptmp*ptmp + + Ai[j]*theta*4.0/Bbetai[j]*(atmp-1.0)*pptmp); + + double d2triagtmpd2delta = + bi[j] * (triagtmpm1 * d2triagddelta2 + + (bi[j]-1.0)*triagtmpm1/triag*dtriagddelta*dtriagddelta); + + double ctmp = (triagtmp * (2.0*dphiddelta + delta*d2phiddelta2) + + 2.0*dtriagtmpddelta*(phi + delta * dphiddelta) + + d2triagtmpd2delta * delta * phi); + + val += ni[i] * ctmp; + } + + return val; +} + +/* + * Calculate d_phi0_d(delta), the first derivative of phi0 + * wrt delta + * + * tau = dimensionless temperature + * delta = dimensionless pressure + */ +double WaterPropsIAPWSphi::phi0_dd() const { + double delta = DELTAsave; + return (-1.0/(delta*delta)); +} + +/* + * Calculate the dPhidDelta function, which is basically the derivative + * of helmholtz free energy wrt delta + * Eqn. (6.4) + */ +double WaterPropsIAPWSphi::phi_dd(double tau, double delta) { + tdpolycalc(tau, delta); + double nau = phi0_dd(); + double res = phiR_dd(); + double retn = nau + res; + return retn; +} + +/* + * Calculate d_phi0/d(tau) + */ +double WaterPropsIAPWSphi::phi0_t() const { + double tau = TAUsave; + double retn = ni0[2] + ni0[3]/tau;; + retn += (ni0[4] * gammi0[4] * (1.0/(1.0 - exp(-gammi0[4]*tau)) - 1.0)); + retn += (ni0[5] * gammi0[5] * (1.0/(1.0 - exp(-gammi0[5]*tau)) - 1.0)); + retn += (ni0[6] * gammi0[6] * (1.0/(1.0 - exp(-gammi0[6]*tau)) - 1.0)); + retn += (ni0[7] * gammi0[7] * (1.0/(1.0 - exp(-gammi0[7]*tau)) - 1.0)); + retn += (ni0[8] * gammi0[8] * (1.0/(1.0 - exp(-gammi0[8]*tau)) - 1.0)); + return retn; +} + +/* + * Calculate Eqn. 6.6 for dphiRdtau, the derivative residual part of the + * dimensionless Helmholtz free energy wrt temperature + * + * tau = dimensionless temperature + * delta = dimensionless pressure + */ +double WaterPropsIAPWSphi::phiR_t() const { + double tau = TAUsave; + double delta = DELTAsave; + int i, j; + double atmp, tmp; + + /* + * Write out the first seven polynomials in the expression + */ + double T375 = pow(tau, 0.375); + double val = ((-0.5) *ni[1] * delta / TAUsqrt / tau + + ni[2] * delta * 0.875 / TAUsqrt * T375 + + ni[3] * delta + + ni[4] * DELTAp[2] * 0.5 / TAUsqrt + + ni[5] * DELTAp[2] * 0.75 * T375 * T375 / tau + + ni[6] * DELTAp[3] * 0.375 * T375 / tau + + ni[7] * DELTAp[4]); + /* + * Next, do polynomial contributions 8 to 51 + */ + for (i = 8; i <= 51; i++) { + tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]-1] * exp(-DELTAp[ciR[i]])); + val += tiR[i] * tmp; + } + + /* + * Next do contributions 52 to 54 + */ + for (j = 0; j < 3; j++) { + i = 52 + j; + double dtmp = delta - epsi[j]; + double ttmp = tau - gammai[j]; + tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * + exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp)); + val += tmp *(tiR[i]/tau - 2.0 * betai[j]*ttmp); + } + + /* + * Next do contributions 55 and 56 + */ + for (j = 0; j < 2; j++) { + i = 55 + j; + double deltam1 = delta - 1.0; + double dtmp2 = deltam1 * deltam1; + atmp = 0.5 / Bbetai[j]; + double theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp); + double triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]); + double ttmp = tau - 1.0; + + double triagtmp = pow(triag, bi[j]); + + double phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp); + + + double dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag; + + double dphidtau = - 2.0 * Di[j] * ttmp * phi; + + val += ni[i] * delta * (dtriagtmpdtau * phi + triagtmp * dphidtau); + } + + return val; +} + +/** + * Calculate the dPhidtau function, which is basically the derivative + * of helmholtz free energy wrt tau + * Eqn. (6.4) + */ +double WaterPropsIAPWSphi::phi_t(double tau, double delta) { + tdpolycalc(tau, delta); + double nau = phi0_t(); + double res = phiR_t(); + double retn = nau + res; + return retn; +} + +/** + * Calculate d2_phi0/dtau2 + */ +double WaterPropsIAPWSphi::phi0_tt() const { + double tau = TAUsave; + double tmp, itmp; + double retn = - ni0[3]/(tau * tau); + for (int i = 4; i <= 8; i++) { + tmp = exp(-gammi0[i]*tau); + itmp = 1.0 - tmp; + retn -= (ni0[i] * gammi0[i] * gammi0[i] * tmp / (itmp * itmp)); + } + return retn; +} + +/** + * Calculate Eqn. 6.6 for dphiRdtau, the second derivative residual part of the + * dimensionless Helmholtz free energy wrt temperature + * + * tau = dimensionless temperature + * delta = dimensionless pressure + */ +double WaterPropsIAPWSphi::phiR_tt() const { + double tau = TAUsave; + double delta = DELTAsave; + int i, j; + double atmp, tmp; + + /* + * Write out the first seven polynomials in the expression + */ + double T375 = pow(tau, 0.375); + double val = ((-0.5) * (-1.5) * ni[1] * delta / (TAUsqrt * tau * tau) + + ni[2] * delta * 0.875 * (-0.125) * T375 / (TAUsqrt * tau) + + ni[4] * DELTAp[2] * 0.5 * (-0.5)/ (TAUsqrt * tau) + + ni[5] * DELTAp[2] * 0.75 *(-0.25) * T375 * T375 / (tau * tau) + + ni[6] * DELTAp[3] * 0.375 *(-0.625) * T375 / (tau * tau)); + /* + * Next, do polynomial contributions 8 to 51 + */ + for (i = 8; i <= 51; i++) { + if (tiR[i] > 1) { + tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]-2] * exp(-DELTAp[ciR[i]])); + val += tiR[i] * (tiR[i] - 1.0) * tmp; + } + } + + /* + * Next do contributions 52 to 54 + */ + for (j = 0; j < 3; j++) { + i = 52 + j; + double dtmp = delta - epsi[j]; + double ttmp = tau - gammai[j]; + tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * + exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp)); + atmp = tiR[i]/tau - 2.0 * betai[j]*ttmp; + val += tmp *(atmp * atmp - tiR[i]/(tau*tau) - 2.0*betai[j]); + } + + /* + * Next do contributions 55 and 56 + */ + for (j = 0; j < 2; j++) { + i = 55 + j; + double deltam1 = delta - 1.0; + double dtmp2 = deltam1 * deltam1; + atmp = 0.5 / Bbetai[j]; + double theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp); + double triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]); + double ttmp = tau - 1.0; + + double triagtmp = pow(triag, bi[j]); + double triagtmpM1 = triagtmp / triag; + + double phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp); + + + double dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag; + + double dphidtau = - 2.0 * Di[j] * ttmp * phi; + + double d2triagtmpdtau2 = + (2 * bi[j] * triagtmpM1 + + 4 * theta * theta * bi[j] * (bi[j]-1.0) * triagtmpM1 / triag); + + double d2phidtau2 = 2.0*Di[j]*phi *(2.0*Di[j]*ttmp*ttmp - 1.0); + + tmp = (d2triagtmpdtau2 * phi + + 2 * dtriagtmpdtau * dphidtau + + triagtmp * d2phidtau2); + val += ni[i] * delta * tmp; + } + + return val; +} + +/* + * Calculate the d2Phidtau2 function, which is basically the second derivative + * of helmholtz free energy wrt tau + * Eqn. (6.4) + */ +double WaterPropsIAPWSphi::phi_tt(double tau, double delta) { + tdpolycalc(tau, delta); + double nau = phi0_tt(); + double res = phiR_tt(); + double retn = nau + res; + return retn; +} + +/** + * Calculate d2_phi0/dtauddelta + */ +double WaterPropsIAPWSphi::phi0_dt() const { + return 0.0; +} + +/* + * Calculate d2_phiR_d(delta)d(tau), the mixed derivative of phi + * wrt delta and tau. + * + * tau = dimensionless temperature + * delta = dimensionless pressure + */ +double WaterPropsIAPWSphi::phiR_dt() const { + double tau = TAUsave; + double delta = DELTAsave; + int i, j; + double tmp; + /* + * Write out the first seven polynomials in the expression + */ + double T375 = pow(tau, 0.375); + double val = (ni[1] * (-0.5) / (TAUsqrt * tau) + + ni[2] * (0.875) * T375 / TAUsqrt + + ni[3] + + ni[4] * 2.0 * delta * (0.5) / TAUsqrt + + ni[5] * 2.0 * delta * (0.75) * T375 * T375 / tau + + ni[6] * 3.0 * DELTAp[2] * 0.375 * T375 / tau + + ni[7] * 4.0 * DELTAp[3]); + /* + * Next, do polynomial contributions 8 to 51 + */ + for (i = 8; i <= 51; i++) { + tmp = (ni[i] * tiR[i] * exp(-DELTAp[ciR[i]]) * DELTAp[diR[i] - 1] * + TAUp[tiR[i] - 1]); + val += tmp * (diR[i] - ciR[i] * DELTAp[ciR[i]]); + } + + /* + * Next do contributions 52 to 54 + */ + for (j = 0; j < 3; j++) { + i = 52 + j; + double dtmp = delta - epsi[j]; + double ttmp = tau - gammai[j]; + tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * + exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp)); + val += tmp * ((diR[i]/delta - 2.0 * alphai[j] * dtmp) * + (tiR[i]/tau - 2.0 * betai[j] * ttmp)); + } + + /* + * Next do contributions 55 and 56 + */ + for (j = 0; j < 2; j++) { + i = 55 + j; + double deltam1 = delta - 1.0; + double dtmp2 = deltam1 * deltam1; + double atmp = 0.5 / Bbetai[j]; + double theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp); + double triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]); + double ttmp = tau - 1.0; + + double triagtmp = pow(triag, bi[j]); + double triagtmpm1 = pow(triag, bi[j]-1.0); + double atmpM1 = atmp - 1.0; + double ptmp = pow(dtmp2,atmpM1); + double p2tmp = pow(dtmp2, ai[j]-1.0); + double dtriagddelta = + deltam1 *(Ai[j] * theta * 2.0 / Bbetai[j] * ptmp + + 2.0*Bi[j]*ai[j]*p2tmp); + + double phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp); + double dphiddelta = -2.0*Ci[j]*deltam1*phi; + double dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta; + + + double dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag; + + double dphidtau = - 2.0 * Di[j] * ttmp * phi; + + double d2phiddeltadtau = 4.0 * Ci[j] * Di[j] * deltam1 * ttmp * phi; + + double d2triagtmpddeltadtau = + ( -Ai[j] * bi[j] * 2.0 / Bbetai[j] * triagtmpm1 * deltam1 * ptmp + -2.0 * theta * bi[j] * (bi[j] - 1.0) * triagtmpm1 / triag * dtriagddelta); + + + double tmp = ni[i] * (triagtmp * (dphidtau + delta*d2phiddeltadtau) + + delta * dtriagtmpddelta * dphidtau + + dtriagtmpdtau * (phi + delta * dphiddelta) + + d2triagtmpddeltadtau * delta * phi); + val += tmp; + } + + return val; +} + +/** + * This program computes the reduced density, given the reduced pressure + * and the reduced temperature, tau. It takes an initial guess, deltaGuess. + * DeltaGuess is important as this is a multivalued function below the + * critical point. + * + */ +double WaterPropsIAPWSphi::dfind(double p_red, double tau, double deltaGuess) { + double dd = deltaGuess; + bool conv = false; + double deldd = dd; + double pcheck = 1.0E-30 + 1.0E-8 * p_red; + for (int n = 0; n < 100; n++) { + /* + * Calculate the internal polynomials, and then calculate the + * phi deriv functions needed by this routine. + */ + tdpolycalc(tau, dd); + double q1 = phiR_d(); + double q2 = phiR_dd(); + + /* + * Calculate the predicted reduced pressure, pred0, based on the + * current tau and dd. + */ + double pred0 = dd + dd * dd * q1; + /* + * Calculate the derivative of the predicted reduced pressure + * wrt the reduced density, dd, This is dpddelta + */ + double dpddelta = 1.0 + 2.0 * dd * q1 + dd * dd * q2; + /* + * If dpddelta is negative, then we are in the middle of the + * 2 phase region, beyond the stability curve. We need to adjust + * the initial guess outwards and start a new iteration. + */ + if (dpddelta <= 0.0) { + if (deltaGuess > 1.0) dd = dd * 1.05; + if (deltaGuess < 1.0) dd = dd * 0.95; + continue; + } + /* + * Check for convergence + */ + if (fabs(pred0-p_red) < pcheck) { + conv = true; + break; + } + + /* + * Dampen and crop the update + */ + double dpdx = dpddelta; + if (n < 10) { + dpdx = dpddelta * 1.1; + } + if (dpdx < 0.1) dpdx = 0.1; + + /* + * Formulate the update to reduced density using + * Newton's method. Then, crop it to a max value + * of 0.02 + */ + deldd = - (pred0 - p_red) / dpdx; + if (fabs(deldd) > 0.05) { + deldd = deldd * 0.05 / fabs(deldd); + } + /* + * updated the reduced density value + */ + dd = dd + deldd; + if (fabs(deldd/dd) < 1.0E-14) { + conv = true; + break; + } + /* + * Check for negative densities + */ + if (dd <= 0.0) { + dd = 1.0E-24; + } + } + /* + * Check for convergence, and return 0.0 if it wasn't achieved. + */ + if (! conv) { + dd = 0.0; + } + return dd; +} + +/** + * Calculate the dimensionless gibbs free energy g/RT. + */ +double WaterPropsIAPWSphi::gibbs_RT() const { + double delta = DELTAsave; + double rd = phiR_d(); + double g = 1.0 + phi0() + phiR() + delta * rd; + return g; +} + +/** + * Calculate the dimensionless enthalpy h/RT. + */ +double WaterPropsIAPWSphi::enthalpy_RT() const { + double delta = DELTAsave; + double tau = TAUsave; + double rd = phiR_d(); + double nt = phi0_t(); + double rt = phiR_t(); + double hRT = 1.0 + tau * (nt + rt) + delta * rd; + return hRT; +} + +/** + * Calculate the dimensionless entropy s/R. + */ +double WaterPropsIAPWSphi::entropy_R() const { + double tau = TAUsave; + double nt = phi0_t(); + double rt = phiR_t(); + double p0 = phi0(); + double pR = phiR(); + double sR = tau * (nt + rt) - p0 - pR; + return sR; +} + +/** + * Calculate the dimensionless internal energy, u/RT. + */ +double WaterPropsIAPWSphi::intEnergy_RT() const { + double tau = TAUsave; + double nt = phi0_t(); + double rt = phiR_t(); + double uR = tau * (nt + rt); + return uR; +} + +/** + * Calculate the dimensionless constant volume Heat Capacity, Cv/R + */ +double WaterPropsIAPWSphi::cv_R() const { + double tau = TAUsave; + double ntt = phi0_tt(); + double rtt = phiR_tt(); + double cvR = - tau * tau * (ntt + rtt); + return cvR; +} + +/** + * Calculate the dimensionless constant pressure Heat Capacity, Cp/R + */ +double WaterPropsIAPWSphi::cp_R() const { + double tau = TAUsave; + double delta = DELTAsave; + double cvR = cv_R(); + //double nd = phi0_d(); + double rd = phiR_d(); + double rdd = phiR_dd(); + double rdt = phiR_dt(); + double num = (1.0 + delta * rd - delta * tau * rdt); + double cpR = cvR + (num * num / + (1.0 + 2.0 * delta * rd + delta * delta * rdd)); + return cpR; +} + diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.h b/Cantera/src/thermo/WaterPropsIAPWSphi.h new file mode 100644 index 000000000..9a9eeb463 --- /dev/null +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.h @@ -0,0 +1,111 @@ +/** + * @file WaterPropsIAPWSphi.h + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ + +#ifndef WATERPROPSIAPWSPHI_H +#define WATERPROPSIAPWSPHI_H + +/* + * Units Note: This class works with reduced units exclusively. + */ + +class WaterPropsIAPWSphi { + +public: + WaterPropsIAPWSphi(); + + /* + * Calculate the base phi's, recalculating the internal polynomials + */ + double phi(double tau, double delta); + double phi_d(double tau, double delta); + double phi_dd(double tau, double delta); + double phi_t(double tau, double delta); + double phi_tt(double tau, double delta); + double phi_dt(double tau, double delta); + + void check1(); + void check2(); + + /** + * Calculate the dimensionless pressure, pred: + * pred = pressure M / (rho RT) + */ + double pressure_rhoRT(double tau, double delta); + + /** + * This program computes the reduced density, given the reduced pressure + * and the reduced temperature, tau. It takes an initial guess, deltaGuess. + * DeltaGuess is important as this is a multivalued function below the + * critical point. + */ + double dfind(double p_red, double tau, double deltaGuess); + + /** + * Calculate the dimensionless gibbs free energy + */ + double gibbs_RT() const; + + /** + * Calculate the dimensionless enthalpy, h/RT + */ + double enthalpy_RT() const; + + /** + * Calculate the dimensionless entropy, s/R + */ + double entropy_R() const; + + /** + * Calculate the dimensionless internal energy, u/RT + */ + double intEnergy_RT() const; + + /** + * Calculate the dimensionless constant volume heat capacity, Cv/R + */ + double cv_R() const; + + /** + * Calculate the dimensionless constant pressure heat capacity, Cv/R + */ + double cp_R() const; + + /** + * 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. + */ + void tdpolycalc(double tau, double delta); + + double phiR() const; +private: + double phi0() const; + double phiR_d() const; + double phi0_d() const; + double phiR_dd() const; + double phi0_dd() const; + double phi0_t() const; + double phiR_t() const; + double phiR_tt() const; + double phi0_tt() const; + double phiR_dt() const; + double phi0_dt() const; + void intCheck(double tau, double delta); + +protected: + double TAUp[52]; + double DELTAp[16]; + double TAUsave; + double TAUsqrt; + double DELTAsave; +}; +#endif diff --git a/Cantera/src/thermo/WaterTP.cpp b/Cantera/src/thermo/WaterTP.cpp new file mode 100644 index 000000000..18ae092f1 --- /dev/null +++ b/Cantera/src/thermo/WaterTP.cpp @@ -0,0 +1,425 @@ +/** + * @file WaterTP.cpp + * + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ + +#include "xml.h" +#include "WaterTP.h" +#include "WaterPropsIAPWS.h" +#include "importCTML.h" + +namespace Cantera { + /** + * Basic list of constructors and duplicators + */ + + WaterTP::WaterTP() : + ThermoPhase(), + m_sub(0), + m_subflag(0), + m_mw(0.0), + EW_Offset(0.0), + SW_Offset(0.0), + m_verbose(0), + m_allowGasPhase(false) + { + constructPhase(); + } + + + WaterTP::WaterTP(string inputFile, string id) : + ThermoPhase(), + m_sub(0), + m_subflag(0), + m_mw(0.0), + EW_Offset(0.0), + SW_Offset(0.0), + m_verbose(0), + m_allowGasPhase(false) + { + constructPhaseFile(inputFile, id); + } + + + WaterTP::WaterTP(XML_Node& phaseRoot, string id) : + ThermoPhase(), + m_sub(0), + m_subflag(0), + m_mw(0.0), + EW_Offset(0.0), + SW_Offset(0.0), + m_verbose(0), + m_allowGasPhase(false) + { + constructPhaseXML(phaseRoot, id) ; + } + + + + WaterTP::WaterTP(const WaterTP &b) : + ThermoPhase(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_allowGasPhase(b.m_allowGasPhase) + { + m_sub = new WaterPropsIAPWS(*(b.m_sub)); + /* + * Use the assignment operator to do the brunt + * of the work for the copy construtor. + */ + *this = b; + } + + /** + * Assignment operator + */ + WaterTP& WaterTP::operator=(const WaterTP&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_allowGasPhase = b.m_allowGasPhase; + return *this; + } + + + ThermoPhase *WaterTP::duplMyselfAsThermoPhase() { + WaterTP* wtp = new WaterTP(*this); + return (ThermoPhase *) wtp; + } + + WaterTP::~WaterTP() { + delete m_sub; + } + + + + void WaterTP::constructPhase() { + throw CanteraError("constructPhaseXML", "unimplemented"); + + } + + + /** + * constructPhase: + * + * Initialization of a Debye-Huckel phase using an + * xml file. + * + * This routine is a precursor to initThermo(XML_Node*) + * routine, which does most of the work. + * + * @param infile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void WaterTP::constructPhaseXML(XML_Node& phaseNode, string id) { + + /* + * Call the Cantera importPhase() function. This will import + * all of the species into the phase. This will also handle + * all of the solvent and solute standard states. + */ + bool m_ok = importPhase(phaseNode, this); + if (!m_ok) { + throw CanteraError("initThermo","importPhase failed "); + } + + } + + + + + /** + * initThermo(): + * + * Initialization of a Debye-Huckel phase using an + * xml file. + * + * This routine is a precursor to initThermo(XML_Node*) + * routine, which does most of the work. + * + * @param infile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void WaterTP::constructPhaseFile(string inputFile, string id) { + + if (inputFile.size() == 0) { + throw CanteraError("WaterTp::initThermo", + "input file is null"); + } + string path = findInputFile(inputFile); + ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("WaterTP::initThermo","could not open " + +path+" for reading."); + } + /* + * The phase object automatically constructs an XML object. + * Use this object to store information. + */ + XML_Node &phaseNode_XML = xml(); + XML_Node *fxml = new XML_Node(); + fxml->build(fin); + XML_Node *fxml_phase = findXMLPhase(fxml, id); + if (!fxml_phase) { + throw CanteraError("WaterTP::initThermo", + "ERROR: Can not find phase named " + + id + " in file named " + inputFile); + } + fxml_phase->copy(&phaseNode_XML); + constructPhaseXML(*fxml_phase, id); + delete fxml; + } + + + + void WaterTP::initThermo() { + + + + + } + + void WaterTP:: + initThermoXML(XML_Node& phaseNode, string id) { + if (m_sub) delete m_sub; + m_sub = new WaterPropsIAPWS(); + if (m_sub == 0) { + throw CanteraError("WaterTP::initThermo", + "could not create new substance object."); + } + /* + * Calculate the molecular weight. Note while there may + * be a very good calculated weight in the steam table + * class, using this weight may lead to codes exhibiting + * mass loss issues. We need to grab the elemental + * atomic weights used in the Element class and calculate + * a consistent H2O molecular weight based on that. + */ + int nH = elementIndex("H"); + if (nH < 0) { + throw CanteraError("WaterTP::initThermo", + "H not an element"); + } + double mw_H = atomicWeight(nH); + int nO = elementIndex("O"); + if (nO < 0) { + throw CanteraError("WaterTP::initThermo", + "O not an element"); + } + double mw_O = atomicWeight(nO); + m_mw = 2.0 * mw_H + mw_O; + m_weight[0] = m_mw; + setMolecularWeight(0,m_mw); + double one = 1.0; + setMoleFractions(&one); + + /* + * Set the baseline + */ + doublereal T = 298.15; + + doublereal presLow = 1.0E-2; + doublereal oneBar = 1.0E5; + doublereal dens = density(); + doublereal dd = m_sub->density(T, presLow, WATER_GAS, dens); + setTemperature(T); + setDensity(dd); + SW_Offset = 0.0; + doublereal s = entropy_mole(); + s -= GasConstant * log(oneBar/presLow); + if (s != 188.835E3) { + SW_Offset = 188.835E3 - s; + } + s = entropy_mole(); + s -= GasConstant * log(oneBar/presLow); + printf("s = %g\n", s); + + doublereal h = enthalpy_mole(); + if (h != -241.826E6) { + EW_Offset = -241.826E6 - h; + } + h = enthalpy_mole(); + + printf("h = %g\n", h); + + + /* + * Set the initial state of the system to 298.15 K and + * 1 bar. + */ + setTemperature(298.15); + double rho0 = m_sub->density(298.15, OneAtm, WATER_LIQUID); + setDensity(rho0); + + /* + * We have to do something with the thermo function here. + */ + if (m_spthermo) { + delete m_spthermo; + m_spthermo = 0; + } + } + + void WaterTP:: + setParametersFromXML(const XML_Node& eosdata) { + eosdata._require("model","PureFluid"); + m_subflag = atoi(eosdata["fluid_type"].c_str()); + if (m_subflag < 0) + throw CanteraError("WaterTP::setParametersFromXML", + "missing or negative substance flag"); + } + + /** + * Return the molar enthalpy in units of J kmol-1 + */ + doublereal WaterTP:: + enthalpy_mole() const { + double T = temperature(); + double dens = density(); + doublereal h = m_sub->enthalpy(T, dens); + return (h + EW_Offset); + } + + /** + * Calculate the internal energy in mks units of + * J kmol-1 + */ + doublereal WaterTP:: + intEnergy_mole() const { + double T = temperature(); + double dens = density(); + doublereal u = m_sub->intEnergy(T, dens); + return (u + EW_Offset); + } + + /** + * Calculate the entropy in mks units of + * J kmol-1 K-1 + */ + doublereal WaterTP:: + entropy_mole() const { + double T = temperature(); + double dens = density(); + doublereal s = m_sub->entropy(T, dens); + return (s + SW_Offset); + } + + /** + * Calculate the Gibbs free energy in mks units of + * J kmol-1 K-1. + */ + doublereal WaterTP:: + gibbs_mole() const { + double T = temperature(); + double dens = density(); + doublereal g = m_sub->Gibbs(T, dens); + return (g + EW_Offset - SW_Offset*T); + } + + /** + * Calculate the constant pressure heat capacity + * in mks units of J kmol-1 K-1 + */ + doublereal WaterTP:: + cp_mole() const { + double T = temperature(); + double dens = density(); + doublereal cp = m_sub->cp(T, dens); + return cp; + } + + /** + * Calculate the constant volume heat capacity + * in mks units of J kmol-1 K-1 + */ + doublereal WaterTP:: + cv_mole() const { + double T = temperature(); + double dens = density(); + doublereal cv = m_sub->cv(T, dens); + return cv; + } + + /** + * Calculate the pressure (Pascals), given the temperature and density + * Temperature: kelvin + * rho: density in kg m-3 + */ + doublereal WaterTP:: + pressure() const { + double T = temperature(); + double dens = density(); + doublereal p = m_sub->pressure(T, dens); + return p; + } + + void WaterTP:: + setPressure(doublereal p) { + double T = temperature(); + double dens = density(); + int waterState = WATER_GAS; + double rc = m_sub->Rhocrit(); + if (dens > rc) { + waterState = WATER_LIQUID; + } + doublereal dd = m_sub->density(T, p, waterState, dens); + if (dd <= 0.0) { + throw CanteraError("setPressure", "error"); + } + setDensity(dd); + } + + + /// critical temperature + doublereal WaterTP::critTemperature() const { return m_sub->Tcrit(); } + + /// critical pressure + doublereal WaterTP::critPressure() const { return m_sub->Pcrit(); } + + /// critical density + doublereal WaterTP::critDensity() const { return m_sub->Rhocrit(); } + + + + void WaterTP::setTemperature(double temp) { + State::setTemperature(temp); + doublereal dd = density(); + m_sub->setState(temp, dd); + } + + + + /// saturation pressure + doublereal WaterTP::satPressure(doublereal t){ + doublereal pp = m_sub->psat(t); + double dens = density(); + setTemperature(t); + setDensity(dens); + return pp; + } + + + +} diff --git a/Cantera/src/thermo/WaterTP.h b/Cantera/src/thermo/WaterTP.h new file mode 100644 index 000000000..8d643e84d --- /dev/null +++ b/Cantera/src/thermo/WaterTP.h @@ -0,0 +1,196 @@ +/** + * @file WaterTP.h + * + * Declares a ThermoPhase class consisting of + * pure water. + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ + +#ifndef CT_WATERTP_H +#define CT_WATERTP_H + +#include "ThermoPhase.h" + +class WaterPropsIAPWS; + +namespace Cantera { + + + /** + *Class for single-component water. This is designed to cover just the + * liquid part of water. + * + * + * 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 WaterTP : public ThermoPhase { + + public: + + /** + * Basic list of constructors and duplicators + */ + WaterTP(); + WaterTP(const WaterTP &b); + WaterTP& operator=(const WaterTP&b); + WaterTP(string inputFile, string id = ""); + WaterTP(XML_Node& phaseRef, string id = ""); + virtual ~WaterTP(); + ThermoPhase *duplMyselfAsThermoPhase(); + + /** + * + * @name Utilities + * @{ + */ + virtual int eosType() 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 Mechanical Equation of State Properties --------------------- + //@{ + + virtual doublereal pressure() const; + virtual void setPressure(doublereal p); + + /** + * @} + * @name Potential Energy + * @{ + */ + + /** + * @} + * @name Activities, Standard States, and Activity Concentrations + * @{ + */ + + //@} + /// @name Partial Molar Properties of 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 pressure + virtual doublereal critPressure() const; + + /// critical density + virtual doublereal critDensity() const; + + /// saturation temperature + //virtual doublereal satTemperature(doublereal p) const; + + + + /// saturation pressure + virtual doublereal satPressure(doublereal t); + + + virtual void setTemperature(double temp); + + virtual void constructPhase(); + virtual void constructPhaseFile(string inputFile, string id); + virtual void constructPhaseXML(XML_Node& phaseNode, string id); + + + virtual void initThermoXML(XML_Node& eosdata, string id); + virtual void initThermo(); + 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; + +private: + mutable WaterPropsIAPWS *m_sub; + int m_subflag; + doublereal m_mw; + + /** + * 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; + + 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; + }; + +} + +#endif + + +