631 lines
15 KiB
C++
631 lines
15 KiB
C++
/**
|
|
* @file WaterSSTP.cpp
|
|
* Definitions for a %ThermoPhase class consisting of pure water (see \ref thermoprops
|
|
* and class \link Cantera::WaterSSTP WaterSSTP\endlink).
|
|
*/
|
|
/*
|
|
* Copyright (2006) Sandia Corporation. Under the terms of
|
|
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
|
|
* U.S. Government retains certain rights in this software.
|
|
*/
|
|
|
|
#include "cantera/thermo/WaterSSTP.h"
|
|
#include "cantera/thermo/WaterPropsIAPWS.h"
|
|
#include "cantera/thermo/WaterProps.h"
|
|
#include "cantera/thermo/ThermoFactory.h"
|
|
#include "cantera/base/xml.h"
|
|
#include "cantera/base/stringUtils.h"
|
|
|
|
#include <cmath>
|
|
|
|
using namespace std;
|
|
|
|
namespace Cantera
|
|
{
|
|
/**
|
|
* Basic list of constructors and duplicators
|
|
*/
|
|
|
|
WaterSSTP::WaterSSTP() :
|
|
SingleSpeciesTP(),
|
|
m_sub(0),
|
|
m_waterProps(0),
|
|
m_mw(0.0),
|
|
EW_Offset(0.0),
|
|
SW_Offset(0.0),
|
|
m_ready(false),
|
|
m_allowGasPhase(false)
|
|
{
|
|
//constructPhase();
|
|
}
|
|
|
|
|
|
WaterSSTP::WaterSSTP(std::string inputFile, std::string id) :
|
|
SingleSpeciesTP(),
|
|
m_sub(0),
|
|
m_waterProps(0),
|
|
m_mw(0.0),
|
|
EW_Offset(0.0),
|
|
SW_Offset(0.0),
|
|
m_ready(false),
|
|
m_allowGasPhase(false)
|
|
{
|
|
constructPhaseFile(inputFile, id);
|
|
}
|
|
|
|
|
|
WaterSSTP::WaterSSTP(XML_Node& phaseRoot, std::string id) :
|
|
SingleSpeciesTP(),
|
|
m_sub(0),
|
|
m_waterProps(0),
|
|
m_mw(0.0),
|
|
EW_Offset(0.0),
|
|
SW_Offset(0.0),
|
|
m_ready(false),
|
|
m_allowGasPhase(false)
|
|
{
|
|
constructPhaseXML(phaseRoot, id) ;
|
|
}
|
|
|
|
|
|
|
|
WaterSSTP::WaterSSTP(const WaterSSTP& b) :
|
|
SingleSpeciesTP(b),
|
|
m_sub(0),
|
|
m_waterProps(0),
|
|
m_mw(b.m_mw),
|
|
EW_Offset(b.EW_Offset),
|
|
SW_Offset(b.SW_Offset),
|
|
m_ready(false),
|
|
m_allowGasPhase(b.m_allowGasPhase)
|
|
{
|
|
m_sub = new WaterPropsIAPWS(*(b.m_sub));
|
|
m_waterProps = new WaterProps(m_sub);
|
|
|
|
/*
|
|
* Use the assignment operator to do the brunt
|
|
* of the work for the copy constructor.
|
|
*/
|
|
*this = b;
|
|
}
|
|
|
|
/*
|
|
* Assignment operator
|
|
*/
|
|
WaterSSTP& WaterSSTP::operator=(const WaterSSTP& b)
|
|
{
|
|
if (&b == this) {
|
|
return *this;
|
|
}
|
|
m_sub->operator=(*(b.m_sub));
|
|
|
|
if (!m_waterProps) {
|
|
m_waterProps = new WaterProps(m_sub);
|
|
}
|
|
m_waterProps->operator=(*(b.m_waterProps));
|
|
|
|
|
|
m_mw = b.m_mw;
|
|
m_ready = b.m_ready;
|
|
m_allowGasPhase = b.m_allowGasPhase;
|
|
return *this;
|
|
}
|
|
|
|
|
|
ThermoPhase* WaterSSTP::duplMyselfAsThermoPhase() const
|
|
{
|
|
WaterSSTP* wtp = new WaterSSTP(*this);
|
|
return (ThermoPhase*) wtp;
|
|
}
|
|
|
|
WaterSSTP::~WaterSSTP()
|
|
{
|
|
delete m_sub;
|
|
delete m_waterProps;
|
|
}
|
|
|
|
|
|
|
|
|
|
/*
|
|
* @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 WaterSSTP::constructPhaseXML(XML_Node& phaseNode, std::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 ");
|
|
}
|
|
|
|
}
|
|
|
|
/*
|
|
* constructPhaseFile
|
|
*
|
|
*
|
|
* This routine is a precursor to constructPhaseXML(XML_Node*)
|
|
* routine, which does most of the work.
|
|
*
|
|
* @param inputFile 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 WaterSSTP::constructPhaseFile(std::string inputFile, std::string id)
|
|
{
|
|
|
|
if (inputFile.size() == 0) {
|
|
throw CanteraError("WaterSSTP::constructPhaseFile",
|
|
"input file is null");
|
|
}
|
|
std::string path = findInputFile(inputFile);
|
|
std::ifstream fin(path.c_str());
|
|
if (!fin) {
|
|
throw CanteraError("WaterSSTP::constructPhaseFile","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("WaterSSTP::constructPhaseFile",
|
|
"ERROR: Can not find phase named " +
|
|
id + " in file named " + inputFile);
|
|
}
|
|
fxml_phase->copy(&phaseNode_XML);
|
|
constructPhaseXML(*fxml_phase, id);
|
|
delete fxml;
|
|
}
|
|
|
|
|
|
|
|
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) {
|
|
throw CanteraError("WaterSSTP::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.
|
|
*/
|
|
size_t nH = elementIndex("H");
|
|
if (nH == npos) {
|
|
throw CanteraError("WaterSSTP::initThermo",
|
|
"H not an element");
|
|
}
|
|
double mw_H = atomicWeight(nH);
|
|
size_t nO = elementIndex("O");
|
|
if (nO == npos) {
|
|
throw CanteraError("WaterSSTP::initThermo",
|
|
"O not an element");
|
|
}
|
|
double mw_O = atomicWeight(nO);
|
|
m_mw = 2.0 * mw_H + mw_O;
|
|
setMolecularWeight(0,m_mw);
|
|
double one = 1.0;
|
|
setMoleFractions(&one);
|
|
|
|
/*
|
|
* Set the baseline
|
|
*/
|
|
doublereal T = 298.15;
|
|
Phase::setDensity(7.0E-8);
|
|
Phase::setTemperature(T);
|
|
|
|
doublereal presLow = 1.0E-2;
|
|
doublereal oneBar = 1.0E5;
|
|
doublereal dd = m_sub->density(T, presLow, WATER_GAS, 7.0E-8);
|
|
setDensity(dd);
|
|
setTemperature(T);
|
|
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);
|
|
|
|
m_waterProps = new WaterProps(m_sub);
|
|
|
|
|
|
/*
|
|
* We have to do something with the thermo function here.
|
|
*/
|
|
if (m_spthermo) {
|
|
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","PureLiquidWater");
|
|
}
|
|
|
|
/*
|
|
* Return the molar dimensionless enthalpy
|
|
*/
|
|
void WaterSSTP::getEnthalpy_RT(doublereal* hrt) const
|
|
{
|
|
double T = temperature();
|
|
doublereal h = m_sub->enthalpy();
|
|
*hrt = (h + EW_Offset)/(GasConstant*T);
|
|
}
|
|
|
|
/*
|
|
* Calculate the internal energy in mks units of
|
|
* J kmol-1
|
|
*/
|
|
void WaterSSTP::getIntEnergy_RT(doublereal* ubar) const
|
|
{
|
|
doublereal u = m_sub->intEnergy();
|
|
*ubar = (u + EW_Offset)/GasConstant;
|
|
}
|
|
|
|
/*
|
|
* Calculate the dimensionless entropy
|
|
*/
|
|
void WaterSSTP::getEntropy_R(doublereal* sr) const
|
|
{
|
|
doublereal s = m_sub->entropy();
|
|
sr[0] = (s + SW_Offset) / GasConstant;
|
|
}
|
|
|
|
/*
|
|
* Calculate the Gibbs free energy in mks units of
|
|
* J kmol-1 K-1.
|
|
*/
|
|
void WaterSSTP::getGibbs_RT(doublereal* grt) const
|
|
{
|
|
double T = temperature();
|
|
doublereal g = m_sub->Gibbs();
|
|
*grt = (g + EW_Offset - SW_Offset*T) / (GasConstant * T);
|
|
if (!m_ready) {
|
|
throw CanteraError("waterSSTP::", "Phase not ready");
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Calculate the Gibbs free energy in mks units of
|
|
* J kmol-1 K-1.
|
|
*/
|
|
void WaterSSTP::getStandardChemPotentials(doublereal* gss) const
|
|
{
|
|
double T = temperature();
|
|
doublereal g = m_sub->Gibbs();
|
|
*gss = (g + EW_Offset - SW_Offset*T);
|
|
if (!m_ready) {
|
|
throw CanteraError("waterSSTP::", "Phase not ready");
|
|
}
|
|
}
|
|
|
|
void WaterSSTP::getCp_R(doublereal* cpr) const
|
|
{
|
|
doublereal cp = m_sub->cp();
|
|
cpr[0] = cp / GasConstant;
|
|
}
|
|
|
|
/*
|
|
* Calculate the constant volume heat capacity
|
|
* in mks units of J kmol-1 K-1
|
|
*/
|
|
doublereal WaterSSTP::cv_mole() const
|
|
{
|
|
doublereal cv = m_sub->cv();
|
|
return cv;
|
|
}
|
|
|
|
// @name Thermodynamic Values for the Species Reference State
|
|
|
|
|
|
void WaterSSTP::getEnthalpy_RT_ref(doublereal* hrt) const
|
|
{
|
|
doublereal p = pressure();
|
|
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, OneAtm, waterState, dens);
|
|
if (dd <= 0.0) {
|
|
throw CanteraError("setPressure", "error");
|
|
}
|
|
doublereal h = m_sub->enthalpy();
|
|
*hrt = (h + EW_Offset) / (GasConstant * T);
|
|
dd = m_sub->density(T, p, waterState, dens);
|
|
}
|
|
|
|
void WaterSSTP::getGibbs_RT_ref(doublereal* grt) const
|
|
{
|
|
doublereal p = pressure();
|
|
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, OneAtm, waterState, dens);
|
|
if (dd <= 0.0) {
|
|
throw CanteraError("setPressure", "error");
|
|
}
|
|
m_sub->setState_TR(T, dd);
|
|
doublereal g = m_sub->Gibbs();
|
|
*grt = (g + EW_Offset - SW_Offset*T)/ (GasConstant * T);
|
|
dd = m_sub->density(T, p, waterState, dens);
|
|
|
|
}
|
|
|
|
void WaterSSTP::getGibbs_ref(doublereal* g) const
|
|
{
|
|
getGibbs_RT_ref(g);
|
|
doublereal rt = _RT();
|
|
for (size_t k = 0; k < m_kk; k++) {
|
|
g[k] *= rt;
|
|
}
|
|
}
|
|
|
|
void WaterSSTP::getEntropy_R_ref(doublereal* sr) const
|
|
{
|
|
doublereal p = pressure();
|
|
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, OneAtm, waterState, dens);
|
|
|
|
if (dd <= 0.0) {
|
|
throw CanteraError("setPressure", "error");
|
|
}
|
|
m_sub->setState_TR(T, dd);
|
|
|
|
doublereal s = m_sub->entropy();
|
|
*sr = (s + SW_Offset)/ (GasConstant);
|
|
dd = m_sub->density(T, p, waterState, dens);
|
|
|
|
}
|
|
|
|
void WaterSSTP::getCp_R_ref(doublereal* cpr) const
|
|
{
|
|
doublereal p = pressure();
|
|
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, OneAtm, waterState, dens);
|
|
m_sub->setState_TR(T, dd);
|
|
if (dd <= 0.0) {
|
|
throw CanteraError("setPressure", "error");
|
|
}
|
|
doublereal cp = m_sub->cp();
|
|
*cpr = cp / (GasConstant);
|
|
dd = m_sub->density(T, p, waterState, dens);
|
|
}
|
|
|
|
void WaterSSTP::getStandardVolumes_ref(doublereal* vol) const
|
|
{
|
|
doublereal p = pressure();
|
|
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, OneAtm, waterState, dens);
|
|
if (dd <= 0.0) {
|
|
throw CanteraError("setPressure", "error");
|
|
}
|
|
*vol = meanMolecularWeight() /dd;
|
|
dd = m_sub->density(T, p, waterState, dens);
|
|
}
|
|
|
|
/*
|
|
* Calculate the pressure (Pascals), given the temperature and density
|
|
* Temperature: kelvin
|
|
* rho: density in kg m-3
|
|
*/
|
|
doublereal WaterSSTP::pressure() const
|
|
{
|
|
doublereal p = m_sub->pressure();
|
|
return p;
|
|
}
|
|
|
|
void WaterSSTP::
|
|
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);
|
|
}
|
|
|
|
// Returns the isothermal compressibility. Units: 1/Pa.
|
|
/*
|
|
* The isothermal compressibility is defined as
|
|
* \f[
|
|
* \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T
|
|
* \f]
|
|
* or
|
|
* \f[
|
|
* \kappa_T = \frac{1}{\rho}\left(\frac{\partial \rho}{\partial P}\right)_T
|
|
* \f]
|
|
*/
|
|
doublereal WaterSSTP::isothermalCompressibility() const
|
|
{
|
|
doublereal val = m_sub->isothermalCompressibility();
|
|
return val;
|
|
}
|
|
|
|
// Return the volumetric thermal expansion coefficient. Units: 1/K.
|
|
/*
|
|
* The thermal expansion coefficient is defined as
|
|
* \f[
|
|
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
|
|
* \f]
|
|
*/
|
|
doublereal WaterSSTP::thermalExpansionCoeff() const
|
|
{
|
|
doublereal val = m_sub->coeffThermExp();
|
|
return val;
|
|
}
|
|
|
|
doublereal WaterSSTP::dthermalExpansionCoeffdT() const
|
|
{
|
|
doublereal pres = pressure();
|
|
doublereal dens_save = density();
|
|
double T = temperature();
|
|
double tt = T - 0.04;
|
|
doublereal dd = m_sub->density(tt, pres, WATER_LIQUID, dens_save);
|
|
if (dd < 0.0) {
|
|
throw CanteraError("WaterSSTP::dthermalExpansionCoeffdT",
|
|
"Unable to solve for the density at T = " + fp2str(tt) + ", P = " + fp2str(pres));
|
|
}
|
|
doublereal vald = m_sub->coeffThermExp();
|
|
m_sub->setState_TR(T, dens_save);
|
|
doublereal val2 = m_sub->coeffThermExp();
|
|
doublereal val = (val2 - vald) / 0.04;
|
|
return val;
|
|
}
|
|
|
|
|
|
// critical temperature
|
|
doublereal WaterSSTP::critTemperature() const
|
|
{
|
|
return m_sub->Tcrit();
|
|
}
|
|
|
|
// critical pressure
|
|
doublereal WaterSSTP::critPressure() const
|
|
{
|
|
return m_sub->Pcrit();
|
|
}
|
|
|
|
// critical density
|
|
doublereal WaterSSTP::critDensity() const
|
|
{
|
|
return m_sub->Rhocrit();
|
|
}
|
|
|
|
|
|
void WaterSSTP::setTemperature(const doublereal temp)
|
|
{
|
|
Phase::setTemperature(temp);
|
|
doublereal dd = density();
|
|
m_sub->setState_TR(temp, dd);
|
|
}
|
|
|
|
void WaterSSTP::setDensity(const doublereal dens)
|
|
{
|
|
Phase::setDensity(dens);
|
|
doublereal temp = temperature();
|
|
m_sub->setState_TR(temp, dens);
|
|
}
|
|
|
|
// saturation pressure
|
|
doublereal WaterSSTP::satPressure(doublereal t) const
|
|
{
|
|
doublereal tsave = temperature();
|
|
doublereal dsave = density();
|
|
doublereal pp = m_sub->psat(t);
|
|
m_sub->setState_TR(tsave, dsave);
|
|
return pp;
|
|
}
|
|
|
|
// Return the fraction of vapor at the current conditions
|
|
doublereal WaterSSTP::vaporFraction() const
|
|
{
|
|
if (temperature() >= m_sub->Tcrit()) {
|
|
double dens = density();
|
|
if (dens >= m_sub->Rhocrit()) {
|
|
return 0.0;
|
|
}
|
|
return 1.0;
|
|
}
|
|
/*
|
|
* If below tcrit we always return 0 from this class
|
|
*/
|
|
return 0.0;
|
|
}
|
|
|
|
|
|
}
|