First commit of water properties routines. These are

under-the-hood routines for calculation of water electrolyte
thermochemistry.
This commit is contained in:
Harry Moffat 2006-07-04 00:01:53 +00:00
parent e7615eca8b
commit 29bd558cc2
14 changed files with 4542 additions and 3 deletions

View file

@ -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.
*/

View file

@ -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

391
Cantera/src/thermo/PDSS.cpp Normal file
View file

@ -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);
}
}

173
Cantera/src/thermo/PDSS.h Normal file
View file

@ -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

View file

@ -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;
}
}

View file

@ -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

View file

@ -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;
}
}

View file

@ -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

View file

@ -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 <math.h>
#include <stdio.h>
#include <stdlib.h>
/*
* 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);
}

View file

@ -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

File diff suppressed because it is too large Load diff

View file

@ -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

View file

@ -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;
}
}

View file

@ -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