research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

WaterSSTP update

Browse source 
Worked on doxygen documentation
  Changed the XML definition
  Worked on making it fall in line with other ThermoPhase objects.
This commit is contained in:
Harry Moffat 2007-03-16 17:22:38 +00:00
parent bec5229624
commit 871ab9f64c
11 changed files with 332 additions and 227 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
Cantera/src/IdealGasPhase.h
View file

@ -278,7 +278,7 @@ namespace Cantera {
* XML_Node * const xs = xc->findNameID("phase", "silane");
* IdealGasPhase *silaneGas = new IdealGasPhase(*xs);
* @endcode
*
* <HR>
* <H2> XML Example </H2>
* <HR>

4
Cantera/src/ThermoPhase.h
View file

@ -1459,7 +1459,3 @@ namespace Cantera {
#endif

4
Cantera/src/thermo/DebyeHuckel.h
View file

@ -523,7 +523,7 @@ namespace Cantera {
* or
*
* @code
* char iFile[80];
* char iFile[80], file_ID[80];
* strcpy(iFile, "DH_NaCl.xml");
* sprintf(file_ID,"%s#NaCl_electrolyte", iFile);
* XML_Node *xm = get_XML_NameID("phase", file_ID, 0);
@ -533,7 +533,7 @@ namespace Cantera {
* or by the following call to importPhase():
*
* @code
* char iFile[80];
* char iFile[80], file_ID[80];
* strcpy(iFile, "DH_NaCl.xml");
* sprintf(file_ID,"%s#NaCl_electrolyte", iFile);
* XML_Node *xm = get_XML_NameID("phase", file_ID, 0);

157
Cantera/src/thermo/WaterPropsIAPWS.cpp
View file

@ -25,6 +25,8 @@ static const double M_water = 18.015268; // kg kmol-1
/*
* Note, this is the Rgas value quoted in the paper. For consistency
* we have to use that value and not the updated value
*
* The Ratio of R/M = 0.46151805 kJ kg-1 K-1 , which is Eqn. (6.3) in the paper.
*/
//static const double Rgas = 8.314472E3; // Joules kmol-1 K-1
static const double Rgas = 8.314371E3; // Joules kmol-1 K-1
@ -67,15 +69,18 @@ WaterPropsIAPWS::~WaterPropsIAPWS() {
void WaterPropsIAPWS::calcDim(double temperature, double rho) {
tau = T_c / temperature;
delta = rho / Rho_c;
}
/**
* Calculate the Helmholtz Free energy in dimensionless units
*
*/
double WaterPropsIAPWS::helmholtzFE_RT() const{
double retn = m_phi->phi(tau, delta);
return (retn);
/*
* Determine the internal state
*/
if (temperature > T_c) {
iState = WATER_SUPERCRIT;
} else {
if (delta < 1.0) {
iState = WATER_GAS;
} else {
iState = WATER_LIQUID;
}
}
}
/*
@ -84,18 +89,17 @@ double WaterPropsIAPWS::helmholtzFE_RT() const{
*/
double WaterPropsIAPWS::helmholtzFE(double temperature, double rho) {
setState(temperature, rho);
double retn = helmholtzFE_RT();
double retn = m_phi->phi(tau, delta);
double RT = Rgas * temperature;
return (retn * RT);
}
double WaterPropsIAPWS::helmholtzFE() const{
double retn = helmholtzFE_RT();
double retn = m_phi->phi(tau, delta);
double temperature = T_c/tau;
double RT = Rgas * temperature;
return (retn * RT);
}
/*
* Calculate the pressure (Pascals), given the temperature and density
* Temperature: kelvin
@ -103,25 +107,16 @@ double WaterPropsIAPWS::helmholtzFE() const{
*/
double WaterPropsIAPWS::pressure(double temperature, double rho) {
calcDim(temperature, rho);
double retn = pressureM_rhoRT();
double retn = m_phi->pressureM_rhoRT(tau, delta);
return (retn * rho * Rgas * temperature/M_water);
}
double WaterPropsIAPWS::pressure() const{
double retn = pressureM_rhoRT();
double retn = m_phi->pressureM_rhoRT(tau, delta);
double rho = delta * Rho_c;
double temperature = T_c / tau;
return (retn * rho * Rgas * temperature/M_water);
}
/*
* Calculates the pressure in dimensionless form
* pM/(rhoRT) at the currently stored tau and delta values
*/
double WaterPropsIAPWS::pressureM_rhoRT() const {
double retn = m_phi->pressureM_rhoRT(tau, delta);
return retn;
}
/*
* Calculates the density given the temperature and the pressure,
* and a guess at the density. Note, below T_c, this is a
@ -150,7 +145,8 @@ density(double temperature, double pressure, int phase, double rhoguess) {
rhoguess = pressure * M_water / (Rgas * temperature);
} else {
/*
* Provide a guess about the liquid density
* Provide a guess about the liquid density that is
* relatively high -> convergnce from above seems robust.
*/
rhoguess = 1000.;
}
@ -166,7 +162,7 @@ density(double temperature, double pressure, int phase, double rhoguess) {
}
double p_red = pressure * M_water / (Rgas * temperature * Rho_c);
deltaGuess = rhoguess / Rho_c;
calcDim(temperature, rhoguess);
setState(temperature, rhoguess);
double delta_retn = m_phi->dfind(p_red, tau, deltaGuess);
double density_retn;
if (delta_retn >0.0) {
@ -177,17 +173,11 @@ density(double temperature, double pressure, int phase, double rhoguess) {
*/
density_retn = delta_retn * Rho_c;
/*
* Determine the internal state
* Set the internal state -> this may be
* a duplication. However, let's just be sure.
*/
if (temperature > T_c) {
iState = WATER_SUPERCRIT;
} else {
if (delta_retn < 1.0) {
iState = WATER_GAS;
} else {
iState = WATER_LIQUID;
}
}
setState(temperature, density_retn);
} else {
density_retn = -1.0;
@ -327,34 +317,25 @@ isothermalCompressibility(double temperature, double pressure) {
return retn;
}
/**
* Calculate the Gibbs Free energy in dimensionless units
*
*/
double WaterPropsIAPWS::
Gibbs_RT() const{
double gRT = m_phi->gibbs_RT();
return gRT;
}
/**
/*
* Calculate the Gibbs free energy in mks units of
* J kmol-1 K-1.
*/
double WaterPropsIAPWS::
Gibbs(double temperature, double rho) {
setState(temperature, rho);
double gRT = Gibbs_RT();
double gRT = m_phi->gibbs_RT();
return (gRT * Rgas * temperature);
}
double WaterPropsIAPWS::
Gibbs() const {
double gRT = Gibbs_RT();
double gRT = m_phi->gibbs_RT();
double temperature = T_c/tau;
return (gRT * Rgas * temperature);
}
/**
/*
* Calculate the Gibbs free energy in mks units of
* J kmol-1 K-1.
*/
@ -368,7 +349,7 @@ corr(double temperature, double pressure, double &densLiq,
exit(-1);
}
setState(temperature, densLiq);
double gibbsLiqRT = Gibbs_RT();
double gibbsLiqRT = m_phi->gibbs_RT();
densGas = density(temperature, pressure, WATER_GAS, densGas);
if (densGas <= 0.0) {
@ -376,7 +357,7 @@ corr(double temperature, double pressure, double &densLiq,
exit(-1);
}
setState(temperature, densGas);
double gibbsGasRT = Gibbs_RT();
double gibbsGasRT = m_phi->gibbs_RT();
delGRT = gibbsLiqRT - gibbsGasRT;
}
@ -443,133 +424,91 @@ setState(double temperature, double rho) {
m_phi->tdpolycalc(tau, delta);
}
/**
* Calculate the enthalpy in dimensionless units
*
*/
double WaterPropsIAPWS::
enthalpy_RT() const{
double hRT = m_phi->enthalpy_RT();
return hRT;
}
/**
/*
* Calculate the enthalpy in mks units of
* J kmol-1 K-1.
*/
double WaterPropsIAPWS::
enthalpy(double temperature, double rho) {
setState(temperature, rho);
double hRT = enthalpy_RT();
double hRT = m_phi->enthalpy_RT();
return (hRT * Rgas * temperature);
}
double WaterPropsIAPWS::
enthalpy() const {
double temperature = T_c/tau;
double hRT = enthalpy_RT();
double hRT = m_phi->enthalpy_RT();
return (hRT * Rgas * temperature);
}
/**
* Calculate the internal Energy in dimensionless units
*
*/
double WaterPropsIAPWS::
intEnergy_RT() const {
double uRT = m_phi->intEnergy_RT();
return uRT;
}
/**
/*
* Calculate the internal Energy in mks units of
* J kmol-1 K-1.
*/
double WaterPropsIAPWS::
intEnergy(double temperature, double rho) {
setState(temperature, rho);
double uRT = intEnergy_RT();
double uRT = m_phi->intEnergy_RT();
return (uRT * Rgas * temperature);
}
double WaterPropsIAPWS::
intEnergy() const{
double temperature = T_c / tau;
double uRT = intEnergy_RT();
double uRT = m_phi->intEnergy_RT();
return (uRT * Rgas * temperature);
}
/**
* Calculate the enthalpy in dimensionless units
*
*/
double WaterPropsIAPWS::
entropy_R() const {
double sR = m_phi->entropy_R();
return sR;
}
/**
/*
* Calculate the enthalpy in mks units of
* J kmol-1 K-1.
*/
double WaterPropsIAPWS::
entropy(double temperature, double rho) {
setState(temperature, rho);
double sR = entropy_R();
double sR = m_phi->entropy_R();
return (sR * Rgas);
}
/**
/*
* Calculate the enthalpy in mks units of
* J kmol-1 K-1.
*/
double WaterPropsIAPWS::
entropy() const {
double sR = entropy_R();
double sR = m_phi->entropy_R();
return (sR * Rgas);
}
/**
* Calculate the dimensionless Heat capacity at constant volume
*/
double WaterPropsIAPWS::
cv_R() const {
double cvR = m_phi->cv_R();
return cvR;
}
/**
/*
* Calculate heat capacity at constant volume
* J kmol-1 K-1.
*/
double WaterPropsIAPWS::
cv(double temperature, double rho) {
setState(temperature, rho);
double cvR = cv_R();
double cvR = m_phi->cv_R();
return (cvR * Rgas);
}
/**
* Calculate the dimensionless Heat capacity at constant pressure
*/
double WaterPropsIAPWS::
cp_R() const {
double cpR = m_phi->cp_R();
return cpR;
}
/**
/*
* Calculate heat capacity at constant pressure
* J kmol-1 K-1.
*/
double WaterPropsIAPWS::
cp(double temperature, double rho) {
setState(temperature, rho);
double cpR = cp_R();
double cpR = m_phi->cp_R();
return (cpR * Rgas);
}
double WaterPropsIAPWS::
cp() const {
double cpR = cp_R();
double cpR = m_phi->cp_R();
return (cpR * Rgas);
}

198
Cantera/src/thermo/WaterPropsIAPWS.h
View file

@ -1,6 +1,6 @@
/**
* @file WaterPropsIAPWS.h
*
* Definitions for a class for calculating the equation of state of water.
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
@ -28,9 +28,41 @@
#define WATER_SUPERCRIT 2
//@}
//! Class for calculating the equation of state of water.
/*!
* Class for calculating the properties of water.
*
* The reference is W. Wagner, A. Prub, "The IAPWS Formulation 1995 for the Themodynamic
* Properties of Ordinary Water Substance for General and Scientific Use,"
* J. Phys. Chem. Ref. Dat, 31, 387, 2002.
*
* This class provides a very complicated polynomial for the specific helmholtz free
* energy of water, as a function of temperature and density.
*
* \f[
* \frac{M\hat{f}(\rho,T)}{R T} = \phi(\delta, \tau) =
* \phi^o(\delta, \tau) + \phi^r(\delta, \tau)
* \f]
*
* where
*
* \f[
* \delta = \rho / \rho_c \mbox{\qquad and \qquad} \tau = T_c / T
* \f]
*
* The following constants are assumed
*
* \f[
* T_c = 647.096\mbox{\ K}
* \f]
* \f[
* \rho_c = 322 \mbox{\ kg\ m$^{-3}$}
* \f]
* \f[
* R/M = 0.46151805 \mbox{\ kJ\ kg$^{-1}$\ K$^{-1}$}
* \f]
*
* The free energy is a unique single-valued function of the temperature and density
* over its entire range.
*
* Note, the base thermodynamic state for this class is the one
* used in the steam tables, i.e., the liquid at the triple point
@ -41,6 +73,59 @@
* - psat(273.16) = 611.655 Pascal
* - rho(273.16, psat) = 999.793 kg m-3
*
* Therefore, to use this class within %Cantera, offsets to u() and s() must be used
* to put the water class onto the same basis as other thermodynamic quantities.
* For example, in the WaterSSTP class, these offsets are calculated in the following way.
* The thermodynamic base state for water is set to the NIST basis here
* by specifying constants EW_Offset and SW_Offset. These offsets are
* calculated on the fly so that the following properties hold:
*
* - Delta_Hfo_idealGas(298.15, 1bar) = -241.826 kJ/gmol
* - So_idealGas(298.15, 1bar) = 188.835 J/gmolK
*
* The offsets are calculated by actually computing the above quantities and then
* calculating the correction factor.
*
* This class provides an interface to the #WaterPropsIAPWSphi class, which actually
* calculates the \f$ \phi^o(\delta, \tau) \f$ and the \f$ \phi^r(\delta, \tau) \f$
* polynomials in dimensionless form.
*
* All thermodynamic results from this class are returned in dimensional form. This
* is because the gas constant (and molecular weight) used within this class is allowed to be potentially
* different than that used elsewhere in %Cantera. Therefore, everything has to be
* in dimensional units. Note, however, the thermodynamic basis is set to that used
* in the steam tables. (u = s = 0 for liquid water at the triple point).
*
* This class is not a %ThermoPhase. However, it does maintain an internal state of
* the object that is dependent on temperature and density. The internal state
* is characterized by an internally storred \f$ \tau\f$ and a \f$ \delta \f$ value,
* and an iState value, which indicates whether the point is a liquid, a gas,
* or a supercritical fluid.
* Along with that the \f$ \tau\f$ and a \f$ \delta \f$ values are polynomials of
* \f$ \tau\f$ and a \f$ \delta \f$ that are kept by the #WaterPropsIAPWSphi class.
* Therefore, whenever \f$ \tau\f$ or \f$ \delta \f$ is changed, the function setState()
* must be called in order for the internal state to be kept up to date.
*
* The class is pretty straightfoward. However, one function deserves mention.
* the #density() function calculates the density that is consistent with
* a particular value of the temperature and pressure. It may therefore be
* multivalued or potentially there may be no answer from this function. It therefore
* takes a phase guess and a density guess as optional parameters. If no guesses are
* supplied to density(), a gas phase guess is assumed. This may or may not be what
* is wanted. Therefore, density() should usually at leat be supplied with a phase
* guess so that it may manufacture an appropriate density guess.
* #density() manufactures the initial density guess, nondimensionalizes everything,
* and then calls #WaterPropsIAPWSphi::dfind(), which does the iterative calculation
* to find the density condition that matches the desired input pressure.
*
* The phase guess defines are located in the .h file. they are
*
* - WATER_GAS
* - WATER_LIQUID
* - WATER_SUPERCRIT
*
* @ingroup thermoprops
*
*/
class WaterPropsIAPWS {
public:
@ -174,10 +259,20 @@ public:
double pressure() const;
//! Calculates the density given the temperature and the pressure,
//! and a guess at the density.
//! and a guess at the density. Sets the internal state.
/*!
* Note, below T_c, this is a multivalued function.
*
* The #density() function calculates the density that is consistent with
* a particular value of the temperature and pressure. It may therefore be
* multivalued or potentially there may be no answer from this function. It therefore
* takes a phase guess and a density guess as optional parameters. If no guesses are
* supplied to density(), a gas phase guess is assumed. This may or may not be what
* is wanted. Therefore, density() should usually at leat be supplied with a phase
* guess so that it may manufacture an appropriate density guess.
* #density() manufactures the initial density guess, nondimensionalizes everything,
* and then calls #WaterPropsIAPWSphi::dfind(), which does the iterative calculation
* to find the density condition that matches the desired input pressure.
*
* @param temperature: Kelvin
* @param pressure : Pressure in Pascals (Newton/m**2)
@ -186,7 +281,8 @@ public:
* @param rhoguess : guessed density of the water
* : -1.0 no guessed density
* @return
* Returns the density
* Returns the density. If an error is encountered in the calculation
* the value of -1.0 is returned.
*/
double density(double temperature, double pressure,
int phase = -1, double rhoguess = -1.0);
@ -197,19 +293,6 @@ public:
*/
double density() const;
//! This function returns an estimated value for the saturation pressure.
/*!
* It does this via a polynomial fit of the vapor pressure curve.
* units = (Pascals)
*
* @param temperature Input temperature (Kelvin)
*
* @return
* Returns the estimated saturation pressure
*/
double psat_est(double temperature);
//! Returns the coefficient of thermal expansion as a function of temperature and pressure.
/*!
* alpha = d (ln V) / dT at constant P.
@ -235,31 +318,19 @@ public:
* returns the isothermal compressibility
*/
double isothermalCompressibility(double temperature, double pressure);
//! Utility routine in the calculation of the saturation pressure
//! This function returns an estimated value for the saturation pressure.
/*!
* @param temperature temperature (kelvin)
* @param pressure pressure (Pascal)
* @param densLiq Output density of liquid
* @param densGas output Density of gas
* @param delGRT output delGRT
*/
void corr(double temperature, double pressure, double &densLiq,
double &densGas, double &delGRT);
* It does this via a polynomial fit of the vapor pressure curve.
* units = (Pascals)
*
* @param temperature Input temperature (Kelvin)
*
* @return
* Returns the estimated saturation pressure
*/
double psat_est(double temperature);
//! Utility routine in the calculation of the saturation pressure
/*!
* @param temperature temperature (kelvin)
* @param pressure pressure (Pascal)
* @param densLiq Output density of liquid
* @param densGas output Density of gas
* @param pcorr output corrected pressure
*/
void corr1(double temperature, double pressure, double &densLiq,
double &densGas, double &pcorr);
//! This function returns the saturation pressure given the
//! temperature as an input parameter.
/*!
@ -297,35 +368,29 @@ private:
*/
void calcDim(double temperature, double rho);
/*
* Dimensionless versions of thermo functions. Note these are
* private, because R value is specific to the class. We only
* show the dimensional functions in the interface.
//! Utility routine in the calculation of the saturation pressure
/*!
* @param temperature temperature (kelvin)
* @param pressure pressure (Pascal)
* @param densLiq Output density of liquid
* @param densGas output Density of gas
* @param delGRT output delGRT
*/
double helmholtzFE_RT() const;
void corr(double temperature, double pressure, double &densLiq,
double &densGas, double &delGRT);
//! Returns the dimensionless gibbs free energy
double Gibbs_RT() const;
//! Utility routine in the calculation of the saturation pressure
/*!
* @param temperature temperature (kelvin)
* @param pressure pressure (Pascal)
* @param densLiq Output density of liquid
* @param densGas output Density of gas
* @param pcorr output corrected pressure
*/
void corr1(double temperature, double pressure, double &densLiq,
double &densGas, double &pcorr);
//! Returns the dimensionless enthalpy
double enthalpy_RT() const;
//! Returns the dimensionless internal energy
double intEnergy_RT() const;
//! Returns the dimensionless entropy
double entropy_R() const;
//! Returns the dimensionless heat capacity at constant volume
double cv_R() const;
//! Returns the dimensionless heat capacity at constant pressure
double cp_R() const;
//! Return the current dimensionless pressure
double pressureM_rhoRT() const;
protected:
private:
//! pointer to the underlying object that does the calculations.
WaterPropsIAPWSphi *m_phi;
@ -346,4 +411,3 @@ protected:
int iState;
};
#endif

16
Cantera/src/thermo/WaterPropsIAPWSphi.h
View file

@ -1,7 +1,7 @@
/**
* @file WaterPropsIAPWSphi.h
*
* Lowest level of the classes which support a real water model.
* This class calculates dimensionless quantitites.
*/
/*
* Copywrite (2006) Sandia Corporation. Under the terms of
@ -18,6 +18,10 @@
/*!
* the WaterPropsIAPSWSphi class support low level calls for
* the real description of water.
*
* The reference is W. Wagner, A. Prub, "The IAPWS Formulation 1995 for the Themodynamic
* Properties of Ordinary Water Substance for General and Scientific Use,"
* J. Phys. Chem. Ref. Dat, 31, 387, 2002.
*
* Units Note: This class works with reduced units exclusively.
*/
@ -80,7 +84,6 @@ public:
//! Internal check # 2
void check2();
//! Calculate the dimensionless pressure at tau and delta;
/*!
*
@ -138,9 +141,10 @@ public:
*/
double cp_R() const;
/**
* Calculates internal polynomials in tau and delta. This
* routine is used to store the internal state of tau and delta
//! Calculates internal polynomials in tau and delta.
/*!
* This routine is used to store the internal state of tau and delta
* for later use by the other routines in the class.
*
* @param tau Dimensionless temperature = T_c/T
@ -186,7 +190,7 @@ private:
*/
void intCheck(double tau, double delta);
protected:
private:
//! Value of internally calculated polynomials of powers of TAU
double TAUp[52];

47
Cantera/src/thermo/WaterSSTP.cpp
View file

@ -1,7 +1,7 @@
/**
* @file WaterSSTP.cpp
* Definitions for the Object WaterSSTP, which creates a
* single species ThermoPhase object for real liquid water.
* Declarations for the object WaterSSTP, which creates a
* single species %ThermoPhase object for real liquid water.
*/
/*
* Copywrite (2006) Sandia Corporation. Under the terms of
@ -25,25 +25,25 @@ namespace Cantera {
WaterSSTP::WaterSSTP() :
SingleSpeciesTP(),
m_sub(0),
m_subflag(0),
m_mw(0.0),
EW_Offset(0.0),
SW_Offset(0.0),
m_verbose(0),
m_ready(false),
m_allowGasPhase(false)
{
constructPhase();
//constructPhase();
}
WaterSSTP::WaterSSTP(std::string inputFile, std::string id) :
SingleSpeciesTP(),
m_sub(0),
m_subflag(0),
m_mw(0.0),
EW_Offset(0.0),
SW_Offset(0.0),
m_verbose(0),
m_ready(false),
m_allowGasPhase(false)
{
constructPhaseFile(inputFile, id);
@ -53,11 +53,11 @@ namespace Cantera {
WaterSSTP::WaterSSTP(XML_Node& phaseRoot, std::string id) :
SingleSpeciesTP(),
m_sub(0),
m_subflag(0),
m_mw(0.0),
EW_Offset(0.0),
SW_Offset(0.0),
m_verbose(0),
m_ready(false),
m_allowGasPhase(false)
{
constructPhaseXML(phaseRoot, id) ;
@ -68,14 +68,14 @@ namespace Cantera {
WaterSSTP::WaterSSTP(const WaterSSTP &b) :
SingleSpeciesTP(b),
m_sub(0),
m_subflag(b.m_subflag),
m_mw(b.m_mw),
EW_Offset(b.EW_Offset),
SW_Offset(b.SW_Offset),
m_verbose(b.m_verbose),
m_ready(false),
m_allowGasPhase(b.m_allowGasPhase)
{
m_sub = new WaterPropsIAPWS(*(b.m_sub));
m_sub = new WaterPropsIAPWS(*(b.m_sub));
/*
* Use the assignment operator to do the brunt
* of the work for the copy construtor.
@ -89,9 +89,9 @@ namespace Cantera {
WaterSSTP& WaterSSTP::operator=(const WaterSSTP&b) {
if (&b == this) return *this;
m_sub->operator=(*(b.m_sub));
m_subflag = b.m_subflag;
m_mw = b.m_mw;
m_verbose = b.m_verbose;
m_ready = b.m_ready;
m_allowGasPhase = b.m_allowGasPhase;
return *this;
}
@ -109,7 +109,7 @@ namespace Cantera {
void WaterSSTP::constructPhase() {
throw CanteraError("constructPhaseXML", "unimplemented");
throw CanteraError("WaterSSTP::constructPhase()", "unimplemented");
}
@ -183,10 +183,16 @@ namespace Cantera {
void WaterSSTP::initThermo() {
SingleSpeciesTP::initThermo();
}
void WaterSSTP::
initThermoXML(XML_Node& phaseNode, std::string id) {
/*
* Do initializations that don't depend on knowing the XML file
*/
initThermo();
if (m_sub) delete m_sub;
m_sub = new WaterPropsIAPWS();
if (m_sub == 0) {
@ -238,7 +244,7 @@ namespace Cantera {
}
s = entropy_mole();
s -= GasConstant * log(oneBar/presLow);
printf("s = %g\n", s);
//printf("s = %g\n", s);
doublereal h = enthalpy_mole();
if (h != -241.826E6) {
@ -246,7 +252,7 @@ namespace Cantera {
}
h = enthalpy_mole();
printf("h = %g\n", h);
//printf("h = %g\n", h);
/*
@ -264,15 +270,16 @@ namespace Cantera {
delete m_spthermo;
m_spthermo = 0;
}
/*
* Set the flag to say we are ready to calculate stuff
*/
m_ready = true;
}
void WaterSSTP::
setParametersFromXML(const XML_Node& eosdata) {
eosdata._require("model","PureFluid");
m_subflag = atoi(eosdata["fluid_type"].c_str());
if (m_subflag < 0)
throw CanteraError("WaterSSTP::setParametersFromXML",
"missing or negative substance flag");
eosdata._require("model","PureLiquidWater");
}
/*
@ -315,6 +322,9 @@ namespace Cantera {
double dens = density();
doublereal g = m_sub->Gibbs(T, dens);
*grt = (g + EW_Offset - SW_Offset*T) / (GasConstant * T);
if (!m_ready) {
throw CanteraError("waterSSTP::", "Phase not ready");
}
}
/*
@ -326,6 +336,9 @@ namespace Cantera {
double dens = density();
doublereal g = m_sub->Gibbs(T, dens);
*gss = (g + EW_Offset - SW_Offset*T);
if (!m_ready) {
throw CanteraError("waterSSTP::", "Phase not ready");
}
}
void WaterSSTP::getCp_R(doublereal* cpr) const {

112
Cantera/src/thermo/WaterSSTP.h
View file

@ -24,19 +24,37 @@ namespace Cantera {
//! Class for single-component water. This is designed to cover just the
//! liquid part of water.
/*!
* The reference is W. Wagner, A. Prub, "The IAPWS Formulation 1995 for the Themodynamic
* Properties of Ordinary Water Substance for General and Scientific Use,"
* J. Phys. Chem. Ref. Dat, 31, 387, 2002.
*
* <HR>
* <H2> Specification of Species Standard %State Properties </H2>
* <HR>
*
* The offsets used in the steam tables are different than NIST's.
* They assume u_liq(TP) = 0.0, s_liq(TP) = 0.0, where TP is the
* triple point conditions:
*
* - u(273.16, rho) = 0.0
* - s(273.16, rho) = 0.0
* - psat(273.16) = 611.655 Pascal
* - rho(273.16, psat) = 999.793 kg m-3
*
* These "steam table" assumptions are used by the WaterPropsIAPWS class.
* Therefore, offsets must be calculated to make the thermodynamic
* properties calculated within this class to be consistent with
* thermo properties within Cantera.
*
* Notes:
* Base state for thermodynamic properties:
*
* The thermodynamic base state for water is set to the NIST basis here
* by specifying constants EW_Offset and SW_Offset. These offsets are
* by specifying constants, #EW_Offset and #SW_Offset, one for energy
* quantities and one for entropy quantities. The offsets are
* specified so that the following properties hold:
*
* Delta_Hfo_gas(298.15) = -241.826 kJ/gmol
* So_gas(298.15, 1bar) = 188.835 J/gmolK
* - Delta_Hfo_idealgas(298.15) = -241.826 kJ/gmol
* - So_idealgas(298.15, 1bar) = 188.835 J/gmolK
*
* (http://webbook.nist.gov)
* ref -> (http://webbook.nist.gov)
*
* The "o" here refers to a hypothetical ideal gas state. The way
* we achieve this in practice is to evaluate at a very low pressure
@ -47,10 +65,69 @@ namespace Cantera {
*
* So(1bar) = S(P0) + RT ln(1bar/P0)
*
* The offsets used in the steam tables are different than NIST's.
* They assume u_liq(TP) = 0.0, s_liq(TP) = 0.0, where TP is the
* triple point conditions.
* <HR>
* <H2> %Application within %Kinetics Managers </H2>
* <HR>
*
* This is unimplemented.
*
* <HR>
* <H2> Instantiation of the Class </H2>
* <HR>
*
* The constructor for this phase is NOT located in the default ThermoFactory
* for %Cantera. However, a new %WaterSSTP object may be created by
* the following code snippets, combined with an XML file given in the
* XML example section.
*
* @code
* WaterSSTP *w = new WaterSSTP("waterSSTPphase.xml","");
* @endcode
*
* or
*
* @code
* char iFile[80], file_ID[80];
* strcpy(iFile, "waterSSTPphase.xml");
* sprintf(file_ID,"%s#water", iFile);
* XML_Node *xm = get_XML_NameID("phase", file_ID, 0);
* WaterSSTP *w = new WaterSSTP(*xm);
* @endcode
*
* or by the following call to importPhase():
*
* @code
* char iFile[80], file_ID[80];
* strcpy(iFile, "waterSSTPphase.xml");
* sprintf(file_ID,"%s#water", iFile);
* XML_Node *xm = get_XML_NameID("phase", file_ID, 0);
* WaterSSTP water;
* importPhase(*xm, &water);
* @endcode
*
* <HR>
* <H2> XML Example </H2>
* <HR>
*
* An example of an XML Element named phase setting up a WaterSSTP object with
* id "water" is given below.
*
* @verbatim
<!-- phase water -->
<phase dim="3" id="water">
<elementArray datasrc="elements.xml">O H </elementArray>
<speciesArray datasrc="#species_data">H2O</speciesArray>
<state>
<temperature units="K">300.0</temperature>
<pressure units="Pa">101325.0</pressure>
</state>
<thermo model="PureLiquidWater"/>
<kinetics model="none"/>
</phase>
@endverbatim
*
* Note the model "PureLiquidWater" indicates the usage of the WaterSSTP object.
*
* @ingroup thermoprops
*
*/
@ -406,14 +483,8 @@ namespace Cantera {
*/
virtual void setParametersFromXML(const XML_Node& eosdata);
protected:
void Set(int n, double x, double y) const;
void setTPXState() const;
void check(doublereal v = 0.0) const;
void reportTPXError() const;
protected:
/**
* @internal
* This internal routine must be overwritten because
@ -422,8 +493,11 @@ namespace Cantera {
void _updateThermo() const;
private:
//! Pointer to the WaterPropsIAPWS that calculates the real properties
//! of water.
mutable WaterPropsIAPWS *m_sub;
int m_subflag;
//! Molecular weight of Water -> Cantera assumption
doublereal m_mw;
/**
@ -442,6 +516,8 @@ namespace Cantera {
bool m_verbose;
bool m_ready;
/**
* Since this phase represents a liquid phase, it's an error to
* return a gas-phase answer. However, if the below is true, then

2
test_problems/cathermo/testWaterTP/output_blessed.txt
View file

@ -1,5 +1,3 @@
s = 188835
h = -2.41826e+08
psat(273.16) = 611.655
Comparisons to NIST: (see http://webbook.nist.gov):

15
test_problems/cathermo/testWaterTP/testWaterSSTP.cpp
View file

@ -4,6 +4,7 @@
#include "stdio.h"
#include "math.h"
#include "WaterSSTP.h"
#include "importCTML.h"
#include <new>
using namespace std;
using namespace Cantera;
@ -23,7 +24,21 @@ int main () {
double pres;
try {
WaterSSTP *w = new WaterSSTP("waterTPphase.xml","");
delete w;
char iFile[80], file_ID[80];
strcpy(iFile, "waterTPphase.xml");
sprintf(file_ID,"%s#water", iFile);
XML_Node *xm = get_XML_NameID("phase", file_ID, 0);
w = new WaterSSTP(*xm);
delete w;
strcpy(iFile, "waterTPphase.xml");
sprintf(file_ID,"%s#water", iFile);
xm = get_XML_NameID("phase", file_ID, 0);
w = new WaterSSTP();
importPhase(*xm, w);
/*
* Print out the triple point conditions

2
test_problems/cathermo/testWaterTP/waterTPphase.xml
View file

@ -10,7 +10,7 @@
<temperature units="K">300.0</temperature>
<pressure units="Pa">101325.0</pressure>
</state>
<thermo model="PureFluid" fluid_type="0"/>
<thermo model="PureLiquidWater"/>
<kinetics model="none"/>
</phase>