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

Interrum update for standing up HKFT standard state

Browse source 
Gibbs free energies are working and checked against excel spread sheet!
This commit is contained in:
Harry Moffat 2008-08-24 01:56:24 +00:00
parent 2d9abce4a1
commit 34bae00f15
13 changed files with 1084 additions and 1997 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

646
Cantera/src/thermo/HKFT_PDSS.cpp
View file

@ -1,646 +0,0 @@
/*
* $Id$
*/
#include "ct_defs.h"
#include "xml.h"
#include "ctml.h"
#include "HKFT_PDSS.h"
#include "WaterPDSS.h"
#include "WaterProps.h"
#include "ThermoPhase.h"
using namespace std;
namespace Cantera {
/**
* Basic list of constructors and duplicators
*/
HKFT_PDSS::HKFT_PDSS(ThermoPhase *tp, int spindex) :
PDSS(tp, spindex)
{
}
HKFT_PDSS::HKFT_PDSS(ThermoPhase *tp, int spindex, std::string inputFile, std::string id) :
PDSS(tp, spindex, inputFile, id)
{
}
HKFT_PDSS::HKFT_PDSS(ThermoPhase *tp, int spindex, XML_Node& phaseRoot, std::string id) :
PDSS(tp, spindex, phaseRoot, id)
{
}
HKFT_PDSS::HKFT_PDSS(const HKFT_PDSS &b) :
PDSS(b)
{
/*
* Use the assignment operator to do the brunt
* of the work for the copy construtor.
*/
*this = b;
}
/**
* Assignment operator
*/
HKFT_PDSS& HKFT_PDSS::operator=(const HKFT_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;
}
/**
* Destructor for the HKFT_PDSS class
*/
HKFT_PDSS::~HKFT_PDSS() {
}
void HKFT_PDSS::constructHKFT_PDSS(ThermoPhase *tp, int spindex) {
initThermo();
}
/**
* constructHKFT_PDSSXML:
*
* 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 HKFT_PDSS::constructHKFT_PDSSXML(ThermoPhase *tp, int spindex,
XML_Node& phaseNode, std::string id) {
initThermo();
}
/**
* constructHKFT_PDSSFile():
*
* 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 HKFT_PDSS::constructHKFT_PDSSFile(ThermoPhase *tp, int spindex,
std::string inputFile, std::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("HKFT_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("HKFT_PDSS::initThermo",
"ERROR: Can not find phase named " +
id + " in file named " + inputFile);
}
constructHKFT_PDSSXML(tp, spindex, *fxml_phase, id);
delete fxml;
}
void HKFT_PDSS::
initThermoXML(XML_Node& phaseNode, std::string id) {
initThermo();
}
void HKFT_PDSS::initThermo() {
/*
* Section to initialize m_Z_pr_tr and m_Y_pr_tr
*/
double temp = 273.15 + 25.;
double pres = OneAtm;
double relepsilon = m_wprops->relEpsilon(temp, pres, 0);
m_Z_pr_tr = -1.0 / relepsilon;
//double m_Z_pr_tr = -0.0127803;
//printf("m_Z_pr_tr = %20.10g\n", m_Z_pr_tr );
double drelepsilondT = m_wprops->relEpsilon(temp, pres, 1);
//double m_Y_pr_tr = -5.799E-5;
m_Y_pr_tr = drelepsilondT / (relepsilon * relepsilon);
//printf("m_Y_pr_tr = %20.10g\n", m_Y_pr_tr );
}
void HKFT_PDSS::
setParametersFromXML(const XML_Node& eosdata) {
}
/**
* Return the molar enthalpy in units of J kmol-1
*/
doublereal HKFT_PDSS::
enthalpy_mole() const {
throw CanteraError("HKFT_PDSS::enthalpy_mole()", "unimplemented");
return (0.0);
}
/**
* Calculate the internal energy in mks units of
* J kmol-1
*/
doublereal HKFT_PDSS::
intEnergy_mole() const {
throw CanteraError("HKFT_PDSS::enthalpy_mole()", "unimplemented");
return (0.0);
}
/**
* Calculate the entropy in mks units of
* J kmol-1 K-1
*/
doublereal HKFT_PDSS::
entropy_mole() const {
throw CanteraError("HKFT_PDSS::entropy_mole()", "unimplemented");
return (0.0);
}
/**
* Calculate the Gibbs free energy in mks units of
* J kmol-1 K-1.
*/
doublereal HKFT_PDSS::gibbs_mole() const {
double val = deltaG();
return (m_Mu0_tr_pr + val);
}
/**
* Calculate the constant pressure heat capacity
* in mks units of J kmol-1 K-1
*/
doublereal HKFT_PDSS::
cp_mole() const {
throw CanteraError("HKFT_PDSS::cp_mole()", "unimplemented");
return (0.0);
}
/**
* Calculate the constant volume heat capacity
* in mks units of J kmol-1 K-1
*/
doublereal HKFT_PDSS::
cv_mole() const {
throw CanteraError("HKFT_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 HKFT_PDSS::
enthalpyDelp_mole() const {
throw CanteraError("HKFT_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 HKFT_PDSS::
intEnergyDelp_mole() const {
throw CanteraError("HKFT_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 HKFT_PDSS::
entropyDelp_mole() const {
throw CanteraError("HKFT_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 HKFT_PDSS::
gibbsDelp_mole() const {
throw CanteraError("HKFT_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 HKFT_PDSS::
cpDelp_mole() const {
throw CanteraError("HKFT_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 HKFT_PDSS::
cvDelp_mole() const {
throw CanteraError("HKFT_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 HKFT_PDSS::
pressure() const {
throw CanteraError("HKFT_PDSS::pressure()", "unimplemented");
return (0.0);
}
void HKFT_PDSS::
setPressure(doublereal p) {
throw CanteraError("HKFT_PDSS::pressure()", "unimplemented");
}
/// critical temperature
doublereal HKFT_PDSS::critTemperature() const {
throw CanteraError("HKFT_PDSS::critTemperature()", "unimplemented");
return (0.0);
}
/// critical pressure
doublereal HKFT_PDSS::critPressure() const {
throw CanteraError("HKFT_PDSS::critPressure()", "unimplemented");
return (0.0);
}
/// critical density
doublereal HKFT_PDSS::critDensity() const {
throw CanteraError("HKFT_PDSS::critDensity()", "unimplemented");
return (0.0);
}
void HKFT_PDSS::setDensity(double dens) {
m_dens = dens;
}
double HKFT_PDSS::density() const {
return m_dens;
}
double HKFT_PDSS::temperature() const {
return m_temp;
}
void HKFT_PDSS::setTemperature(double temp) {
m_temp = temp;
}
doublereal HKFT_PDSS::molecularWeight() const {
return m_mw;
}
void HKFT_PDSS::setMolecularWeight(double mw) {
m_mw = mw;
}
void HKFT_PDSS::setState_TP(double temp, double pres) {
throw CanteraError("HKFT_PDSS::setState_TP()", "unimplemented");
}
/// saturation pressure
doublereal HKFT_PDSS::satPressure(doublereal t){
throw CanteraError("HKFT_PDSS::satPressure()", "unimplemented");
return (0.0);
}
double HKFT_PDSS::deltaG() const {
double pbar = m_pres * 1.0E-5;
double m_presR_bar = OneAtm * 1.0E-5;
double sterm = - m_Entrop_tr_pr * (m_temp - 298.15);
double c1term = -m_c1 * (m_temp * log(m_temp/298.15) - (m_temp - 298.15));
double a1term = m_a1 * (pbar - m_presR_bar);
double a2term = m_a2 * log((2600. + pbar)/(2600. + m_presR_bar));
double c2term = -m_c2 * (( 1.0/(m_temp - 228.) - 1.0/(298.15 - 228.) ) * (228. - m_temp)/228.
- m_temp / (228.*228.) * log( (298.15*(m_temp-228.)) / (m_temp*(298.15-228.)) ));
double a3term = m_a3 / (m_temp - 228.) * (pbar - m_presR_bar);
double a4term = m_a4 / (m_temp - 228.) * log((2600. + pbar)/(2600. + m_presR_bar));
double nu = 166027;
double r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
double gval = gstar(m_temp, m_pres, 0);
double r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
double omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval) );
double relepsilon = m_wprops->relEpsilon(m_temp, m_pres, 0);
double Z = -1.0 / relepsilon;
double wterm = - omega_j * (Z + 1.0);
double wrterm = m_omega_pr_tr * (m_Z_pr_tr + 1.0);
double yterm = m_omega_pr_tr * m_Y_pr_tr * (m_temp - 298.15);
double deltaG_calgmol = sterm + c1term + a1term + a2term + c2term + a3term + a4term + wterm + wrterm + yterm;
// Convert to Joules / kmol
double deltaG = deltaG_calgmol * 1.0E3 * 4.184;
return deltaG;
}
double HKFT_PDSS::electrostatic_radii_calc() {
return 0.0;
}
//! Internal formula for the calculation of a_g()
/*
* The output of this is in units of Angstroms
*/
double HKFT_PDSS::ag(const double temp, const int ifunc) const {
static double ag_coeff[3] = { -2.037662, 5.747000E-3, -6.557892E-6};
if (ifunc == 0) {
double t2 = temp * temp;
double val = ag_coeff[0] + ag_coeff[1] * temp + ag_coeff[2] * t2;
return val;
} else if (ifunc == 1) {
return ag_coeff[1] + ag_coeff[2] * 2.0 * temp;
}
if (ifunc != 2) {
return 0.0;
}
return ag_coeff[2] * 2.0;;
}
//! Internal formula for the calculation of b_g()
/*
* the output of this is unitless
*/
double HKFT_PDSS::bg(const double temp, const int ifunc) const {
static double bg_coeff[3] = { 6.107361, -1.074377E-2, 1.268348E-5};
if (ifunc == 0) {
double t2 = temp * temp;
double val = bg_coeff[0] + bg_coeff[1] * temp + bg_coeff[2] * t2;
return val;
} else if (ifunc == 1) {
return bg_coeff[1] + bg_coeff[2] * 2.0 * temp;
}
if (ifunc != 2) {
return 0.0;
}
return bg_coeff[2] * 2.0;
}
double HKFT_PDSS::f(const double temp, const double pres, const int ifunc) const {
static double af_coeff[3] = { 3.666666E1, -0.1504956E-9, 0.5107997E-13};
double TC = temp - 273.15;
double presBar = pres / 1.0E5;
if (TC < 155.0) return 0.0;
if (TC > 355.0) TC = 355.0;
if (presBar > 1000.) return 0.0;
double T1 = (TC-155.0)/300.;
double fac1;
double p2 = presBar * presBar;
double p3 = presBar * p2;
double p4 = p2 * p2;
double fac2 = af_coeff[1] * p3 + af_coeff[2] * p4;
if (ifunc == 0) {
fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
return fac1 * fac2;
} else if (ifunc == 1) {
fac1 = (4.8 * pow(T1,3.8) + 16.0 * af_coeff[0] * pow(T1, 15.0)) / 300.;
return fac1 * fac2;
} else if (ifunc == 2) {
fac1 = (4.8 * 3.8 * pow(T1,2.8) + 16.0 * 15.0 * af_coeff[0] * pow(T1, 14.0)) / (300. * 300.);
return fac1 * fac2;
} else if (ifunc == 3) {
fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
fac2 = (3.0 * af_coeff[1] * p2 + 4.0 * af_coeff[2] * p3 )/ 1.0E5;
return fac1 * fac2;
} else {
throw CanteraError("HKFT_PDSS::gg", "unimplemented");
}
return 0.0;
}
double HKFT_PDSS::g(const double temp, const double pres, const int ifunc) const {
double afunc = ag(temp, 0);
double bfunc = bg(temp, 0);
m_waterSS->setState_TP(temp, pres);
m_densWaterSS = m_waterSS->density();
// density in gm cm-3
double dens = m_densWaterSS * 1.0E-3;
double gval = afunc * pow((1.0-dens), bfunc);
if (dens >= 1.0) {
return 0.0;
}
if (ifunc == 0) {
return gval;
} else if (ifunc == 1 || ifunc == 2) {
double afuncdT = ag(temp, 1);
double bfuncdT = bg(temp, 1);
double alpha = m_waterSS->thermalExpansionCoeff();
double fac1 = afuncdT * gval / afunc;
double fac2 = bfuncdT * gval * log(1.0 - dens);
double fac3 = gval * alpha * bfunc * dens / (1.0 - dens);
double dgdt = fac1 + fac2 + fac3;
if (ifunc == 1) {
return dgdt;
}
double afuncdT2 = ag(temp, 2);
double bfuncdT2 = bg(temp, 2);
double dfac1dT = dgdt * afuncdT / afunc + afuncdT2 * gval / afunc
- afuncdT * afuncdT * gval / (afunc * afunc);
double ddensdT = - alpha * dens;
double dfac2dT = bfuncdT2 * gval * log(1.0 - dens)
+ bfuncdT * dgdt * log(1.0 - dens)
- bfuncdT * gval /(1.0 - dens) * ddensdT;
double dalphadT = m_waterSS->dthermalExpansionCoeffdT();
double dfac3dT = dgdt * alpha * bfunc * dens / (1.0 - dens)
+ gval * dalphadT * bfunc * dens / (1.0 - dens)
+ gval * alpha * bfuncdT * dens / (1.0 - dens)
+ gval * alpha * bfunc * ddensdT / (1.0 - dens)
- gval * alpha * bfunc * dens / ((1.0 - dens) * (1.0 - dens)) * ddensdT;
return dfac1dT + dfac2dT + dfac3dT;
} else {
throw CanteraError("HKFT_PDSS::gg", "unimplemented");
}
return 0.0;
}
double HKFT_PDSS::gstar(const double temp, const double pres, const int ifunc) const {
double gval = g(temp, pres, ifunc);
double fval = f(temp, pres, ifunc);
return gval - fval;
}
/* awData structure */
/**
* Database for atomic molecular weights
*
* Values are taken from the 1989 Standard Atomic Weights, CRC
*
* awTable[] is a static function with scope limited to this file.
* It can only be referenced via the static Elements class function,
* LookupWtElements().
*
* units = kg / kg-mol (or equivalently gm / gm-mol)
*
* (note: this structure was picked because it's simple, compact,
* and extensible).
*
*/
struct GeData {
char name[4]; ///< Null Terminated name, First letter capitalized
double GeValue; /// < Gibbs free energies of elements J kmol-1
};
//! Values of G_elements(T=298.15,1atm)
/*!
* all units are Joules kmol-1
*/
static struct GeData geDataTable[] = {
{"H", -19.48112E6}, // NIST Webbook - Cox, Wagman 1984
{"Na", -15.29509E6}, // NIST Webbook - Cox, Wagman 1984
{"O", -30.58303E6}, // NIST Webbook - Cox, Wagman 1984
{"Cl", -33.25580E6}, // NIST Webbook - Cox, Wagman 1984
{"Si", -5.61118E6}, // Janaf
{"C", -1.71138E6}, // barin, Knack, NBS Bulletin 1971
{"S", -9.55690E6}, // Yellow - webbook
{"Al", -8.42870E6}, // Webbook polynomial
{"K", -19.26943E6} // Webbook
};
//! Static function to look up Element Free Energies
/*!
*
* This static function looks up the argument string in the
* database above and returns the associated Gibbs Free energies.
*
* @param ElemName String. Only the first 3 characters are significant
*
* @return
* Return value contains the Gibbs free energy for that element
*
* @exception CanteraError
* If a match is not found, a CanteraError is thrown as well
*/
double HKFT_PDSS::LookupGe(const std::string& s) {
int num = sizeof(geDataTable) / sizeof(struct GeData);
string s3 = s.substr(0,3);
for (int i = 0; i < num; i++) {
//if (!std::strncmp(s.c_str(), aWTable[i].name, 3)) {
if (s3 == geDataTable[i].name) {
return (geDataTable[i].GeValue);
}
}
throw CanteraError("LookupGe", "element not found");
return -1.0;
}
void HKFT_PDSS::convertDGFormation() {
/*
* Ok let's get the element compositions and conversion factors.
*/
int ne = m_tp->nElements();
double na;
double ge;
string ename;
double totalSum = 0.0;
for (int m = 0; m < ne; m++) {
na = m_tp->nAtoms(m_spindex, m);
if (na > 0.0) {
ename = m_tp->elementName(m);
ge = LookupGe(ename);
totalSum += na * ge;
}
}
// Add in the charge
if (m_charge_j != 0.0) {
ename = "H";
ge = LookupGe(ename);
totalSum -= m_charge_j * ge;
}
// Ok, now do the calculation. Convert to joules kmol-1
double dg = m_deltaG_formation_tr_pr * 4.184 * 1.0E3;
//! Store the result into an internal variable.
m_Mu0_tr_pr = dg + totalSum;
}
}

260
Cantera/src/thermo/HKFT_PDSS.h
View file

@ -1,260 +0,0 @@
/**
* @file HKFT_PDSS.h
*
* Declares class PDSS pressure dependent standard state
* for a single species
*/
/* $Author$
* $Date$
* $Revision$
*
*
*/
#ifndef CT_HKFT_PDSS_H
#define CT_HKFT_PDSS_H
#include "ct_defs.h"
class XML_Node;
class ThermoPhase;
class WaterPropsIAPWS;
#include "PDSS.h"
namespace Cantera {
class WaterProps;
class WaterPDSS;
/**
* Class for pressure dependent standard states corresponding to
* ionic solutes in electrolyte water.
*
* NOTE: This is largely not done or not complete.
*
*/
class HKFT_PDSS : public PDSS {
public:
/**
* Basic list of constructors and duplicators
*/
HKFT_PDSS(ThermoPhase *tp, int spindex);
HKFT_PDSS(const HKFT_PDSS &b);
HKFT_PDSS& operator=(const HKFT_PDSS&b);
HKFT_PDSS(ThermoPhase *tp, int spindex,
std::string inputFile, std::string id = "");
HKFT_PDSS(ThermoPhase *tp, int spindex,
XML_Node& phaseRef, std::string id = "");
virtual ~HKFT_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);
void constructHKFT_PDSS(ThermoPhase *tp, int spindex);
void constructHKFT_PDSSFile(ThermoPhase *tp, int spindex,
std::string inputFile, std::string id);
void constructHKFT_PDSSXML(ThermoPhase *tp, int spindex,
XML_Node& phaseNode, std::string id);
virtual void initThermoXML(XML_Node& eosdata, std::string id);
virtual void initThermo();
virtual void setParametersFromXML(const XML_Node& eosdata);
private:
//! Main routine that actually calculates the gibbs free energy difference
//! between the reference state at Tr, Pr and T,P
/*!
* This is eEqn. 59 in Johnson et al. (1992).
*
*/
double deltaG() const;
double electrostatic_radii_calc();
double ag(const double temp, const int ifunc = 0) const;
double bg(const double temp, const int ifunc = 0) const;
double g(const double temp, const double pres, const int ifunc = 0) const;
double f(const double temp, const double pres, const int ifunc = 0) const;
double gstar(const double temp, const double pres, const int ifunc = 0) const;
double LookupGe(const std::string& s);
void convertDGFormation();
protected:
private:
//! Water standard state calculator
/*!
* derived from the equation of state for water.
*/
WaterPDSS *m_waterSS;
//! Current value of the pressure for this object
doublereal m_pres;
//! density of standard-state water
/*!
* internal temporary variable
*/
mutable double m_densWaterSS;
/**
* Pointer to the water property calculator
*/
WaterProps *m_waterProps;
//! Born coefficient for the current ion or species
doublereal m_born_coeff_j;
//! Electrostatic radii
doublereal r_e_j;
//! Value of deltaG of Formation at Tr and Pr (cal gmol-1)
/*!
* Tr = 298.15 Pr = 1 atm
*
* This is the delta G for the formation reaction of the
* ion from elements in their stable state at Tr, Pr.
*/
doublereal m_deltaG_formation_tr_pr;
//! Value of the Absolute Gibbs Free Energy NIST scale
/*!
* J kmol-1
*/
doublereal m_Mu0_tr_pr;
//! Value of S_j at Tr and Pr (cal gmol-1 K-1)
/*!
* Tr = 298.15 Pr = 1 atm
*/
doublereal m_Entrop_tr_pr;
//! a1 coefficient (cal gmol-1 bar-1)
doublereal m_a1;
//! a2 coefficient (cal gmol-1)
doublereal m_a2;
//! c1 coefficient (cal gmol-1 K-1)
doublereal m_c1;
//! c2 coefficient (cal K gmol-1)
doublereal m_c2;
//! a3 coefficient (cal K gmol-1 bar-1)
doublereal m_a3;
//! a4 coefficient (cal K gmol-1)
doublereal m_a4;
//! omega_pr_tr coefficient(cal gmol-1)
doublereal m_omega_pr_tr;
//! y = dZdT = 1/(esp*esp) desp/dT
double m_Y_pr_tr;
//double m_Y_pr_tr = -5.799E-5;
double m_Z_pr_tr;
//double m_Z_pr_tr = -0.0127803;
//! Reference pressure is 1 atm in units of bar= 1.0132
doublereal m_presR_bar;
//! Charge of the ion
doublereal m_charge_j;
WaterProps *m_wprops;
};
}
#endif

11
Cantera/src/thermo/HMWSoln_input.cpp
View file

@ -974,16 +974,19 @@ namespace Cantera {
"\" is not allowed");
}
} else {
if (modelString != "constant_incompressible") {
if (modelString != "constant_incompressible" && modelString != "hkft") {
throw CanteraError("HMWSoln::initThermoXML",
"Solute SS Model \"" + modelStringa +
"\" is not known");
}
m_speciesSize[k] = getFloat(*ss, "molarVolume", "-");
if (modelString == "constant_incompressible" ) {
m_speciesSize[k] = getFloat(*ss, "molarVolume", "-");
#ifdef DEBUG_HKM_NOT
cout << "species " << sss[k] << " has volume " <<
m_speciesSize[k] << endl;
cout << "species " << sss[k] << " has volume " <<
m_speciesSize[k] << endl;
#endif
}
// HKM Note, have to fill up m_speciesSize[] for HKFT species
}
}

514
Cantera/src/thermo/PDSS_HKFT.cpp
View file

@ -5,6 +5,8 @@
#include "xml.h"
#include "ctml.h"
#include "PDSS_HKFT.h"
#include "WaterProps.h"
#include "PDSS_Water.h"
#include "VPStandardStateTP.h"
@ -18,20 +20,89 @@ namespace Cantera {
PDSS_HKFT::PDSS_HKFT(VPStandardStateTP *tp, int spindex) :
PDSS(tp, spindex)
PDSS(tp, spindex),
m_waterSS(0),
m_pres(OneAtm),
m_densWaterSS(-1.0),
m_waterProps(0),
m_born_coeff_j(-1.0),
m_r_e_j(-1.0),
m_deltaG_formation_tr_pr(0.0),
m_deltaH_formation_tr_pr(0.0),
m_Mu0_tr_pr(0.0),
m_Entrop_tr_pr(0.0),
m_a1(0.0),
m_a2(0.0),
m_a3(0.0),
m_a4(0.0),
m_c1(0.0),
m_c2(0.0),
m_omega_pr_tr(0.0),
m_Y_pr_tr(0.0),
m_Z_pr_tr(0.0),
m_presR_bar(0.0),
m_charge_j(0.0)
{
m_pdssType = cPDSS_MOLAL_HKFT;
}
PDSS_HKFT::PDSS_HKFT(VPStandardStateTP *tp, int spindex, std::string inputFile, std::string id) :
PDSS(tp, spindex)
PDSS(tp, spindex),
m_waterSS(0),
m_pres(OneAtm),
m_densWaterSS(-1.0),
m_waterProps(0),
m_born_coeff_j(-1.0),
m_r_e_j(-1.0),
m_deltaG_formation_tr_pr(0.0),
m_deltaH_formation_tr_pr(0.0),
m_Mu0_tr_pr(0.0),
m_Entrop_tr_pr(0.0),
m_a1(0.0),
m_a2(0.0),
m_a3(0.0),
m_a4(0.0),
m_c1(0.0),
m_c2(0.0),
m_omega_pr_tr(0.0),
m_Y_pr_tr(0.0),
m_Z_pr_tr(0.0),
m_presR_bar(0.0),
m_charge_j(0.0)
{
m_pdssType = cPDSS_MOLAL_HKFT;
constructPDSSFile(tp, spindex, inputFile, id);
}
PDSS_HKFT::PDSS_HKFT(VPStandardStateTP *tp, int spindex, const XML_Node& speciesNode,
const XML_Node& phaseRoot, bool spInstalled) :
PDSS(tp, spindex)
PDSS(tp, spindex),
m_waterSS(0),
m_pres(OneAtm),
m_densWaterSS(-1.0),
m_waterProps(0),
m_born_coeff_j(-1.0),
m_r_e_j(-1.0),
m_deltaG_formation_tr_pr(0.0),
m_deltaH_formation_tr_pr(0.0),
m_Mu0_tr_pr(0.0),
m_Entrop_tr_pr(0.0),
m_a1(0.0),
m_a2(0.0),
m_a3(0.0),
m_a4(0.0),
m_c1(0.0),
m_c2(0.0),
m_omega_pr_tr(0.0),
m_Y_pr_tr(0.0),
m_Z_pr_tr(0.0),
m_presR_bar(0.0),
m_charge_j(0.0)
{
m_pdssType = cPDSS_MOLAL_HKFT;
// We have to read the info from here
constructPDSSXML(tp, spindex, speciesNode, phaseRoot, spInstalled);
}
PDSS_HKFT::PDSS_HKFT(const PDSS_HKFT &b) :
@ -61,11 +132,14 @@ namespace Cantera {
* Destructor for the PDSS_HKFT class
*/
PDSS_HKFT::~PDSS_HKFT() {
delete m_waterProps;
}
//! Duplicator
PDSS* PDSS_HKFT::duplMyselfAsPDSS() const {
PDSS_HKFT * idg = new PDSS_HKFT(*this);
return (PDSS *) idg;
}
/**
* Return the molar enthalpy in units of J kmol-1
@ -105,12 +179,11 @@ namespace Cantera {
/**
* Calculate the Gibbs free energy in mks units of
* J kmol-1 K-1.
* J kmol-1
*/
doublereal
PDSS_HKFT::gibbs_mole() const {
throw CanteraError("PDSS_HKFT::gibbs_mole()", "unimplemented");
return (0.0);
doublereal PDSS_HKFT::gibbs_mole() const {
double delG = deltaG();
return (m_Mu0_tr_pr + delG);
}
/**
@ -176,13 +249,12 @@ namespace Cantera {
*/
doublereal
PDSS_HKFT::pressure() const {
throw CanteraError("PDSS_HKFT::pressure()", "unimplemented");
return (0.0);
return m_pres;
}
void
PDSS_HKFT::setPressure(doublereal p) {
throw CanteraError("PDSS_HKFT::pressure()", "unimplemented");
m_pres = p;
}
void PDSS_HKFT::setTemperature(doublereal temp) {
@ -225,7 +297,30 @@ namespace Cantera {
PDSS::initThermo();
SpeciesThermo &sp = m_tp->speciesThermo();
m_p0 = sp.refPressure(m_spindex);
m_waterSS = (PDSS_Water *) m_tp->providePDSS(0);
/*
* Section to initialize m_Z_pr_tr and m_Y_pr_tr
*/
double temp = 273.15 + 25.;
double pres = OneAtm;
double relepsilon = m_waterProps->relEpsilon(temp, pres, 0);
m_waterSS->setState_TP(temp, pres);
m_densWaterSS = m_waterSS->density();
m_Z_pr_tr = -1.0 / relepsilon;
//double m_Z_pr_tr = -0.0127803;
//printf("m_Z_pr_tr = %20.10g\n", m_Z_pr_tr );
double drelepsilondT = m_waterProps->relEpsilon(temp, pres, 1);
//double m_Y_pr_tr = -5.799E-5;
m_Y_pr_tr = drelepsilondT / (relepsilon * relepsilon);
//printf("m_Y_pr_tr = %20.10g\n", m_Y_pr_tr );
m_presR_bar = OneAtm / 1.0E5;
m_charge_j = m_tp->charge(m_spindex);
convertDGFormation();
m_waterProps = new WaterProps(m_waterSS);
}
@ -244,19 +339,100 @@ namespace Cantera {
throw CanteraError("PDSS_HKFT::constructPDSSXML", "spInstalled false not handled");
}
const XML_Node *tn = speciesNode.findByName("thermo");
if (!tn) {
throw CanteraError("PDSS_HKFT::constructPDSSXML",
"no thermo Node for species " + speciesNode.name());
}
std::string model = lowercase((*tn)["model"]);
if (model != "hkft") {
throw CanteraError("PDSS_HKFT::initThermoXML",
"thermo model for species isn't hkft: "
+ speciesNode.name());
}
const XML_Node *hh = tn->findByName("HKFT");
if (!hh) {
throw CanteraError("PDSS_HKFT::constructPDSSXML",
"no Thermo::HKFT Node for species " + speciesNode.name());
}
if (hh->hasChild("DG0_f_Pr_Tr")) {
double val = getFloat(*hh, "DG0_f_Pr_Tr");
m_deltaG_formation_tr_pr = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing DG0_f_Pr_Tr field");
}
if (hh->hasChild("DH0_f_Pr_Tr")) {
double val = getFloat(*hh, "DH0_f_Pr_Tr");
m_deltaH_formation_tr_pr = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing DH0_f_Pr_Tr field");
}
if (hh->hasChild("S0_Pr_Tr")) {
double val = getFloat(*hh, "S0_Pr_Tr");
m_Entrop_tr_pr= val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing S0_Pr_Tr field");
}
const XML_Node *ss = speciesNode.findByName("standardState");
if (!ss) {
throw CanteraError("PDSS_HKFT::constructPDSSXML",
"no standardState Node for species " + speciesNode.name());
}
std::string model = (*ss)["model"];
if (model != "constant_incompressible") {
model = lowercase((*ss)["model"]);
if (model != "hkft") {
throw CanteraError("PDSS_HKFT::initThermoXML",
"standardState model for species isn't constant_incompressible: "
"standardState model for species isn't hkft: "
+ speciesNode.name());
}
if (ss->hasChild("a1")) {
double val = getFloat(*ss, "a1");
m_a1 = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a1 field");
}
if (ss->hasChild("a2")) {
double val = getFloat(*ss, "a2");
m_a2 = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a2 field");
}
if (ss->hasChild("a3")) {
double val = getFloat(*ss, "a3");
m_a3 = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a3 field");
}
if (ss->hasChild("a4")) {
double val = getFloat(*ss, "a4");
m_a4 = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a4 field");
}
if (ss->hasChild("c1")) {
double val = getFloat(*ss, "c1");
m_c1 = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing c1 field");
}
if (ss->hasChild("c2")) {
double val = getFloat(*ss, "c2");
m_c2 = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing c2 field");
}
if (ss->hasChild("omega_Pr_Tr")) {
double val = getFloat(*ss, "omega_Pr_Tr");
m_omega_pr_tr = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing omega_Pr_Tr field");
}
std::string id = "";
initThermoXML(phaseNode, id);
}
@ -298,4 +474,302 @@ namespace Cantera {
delete fxml;
}
double PDSS_HKFT::deltaG() const {
double pbar = m_pres * 1.0E-5;
double m_presR_bar = OneAtm * 1.0E-5;
double sterm = - m_Entrop_tr_pr * (m_temp - 298.15);
double c1term = -m_c1 * (m_temp * log(m_temp/298.15) - (m_temp - 298.15));
double a1term = m_a1 * (pbar - m_presR_bar);
double a2term = m_a2 * log((2600. + pbar)/(2600. + m_presR_bar));
double c2term = -m_c2 * (( 1.0/(m_temp - 228.) - 1.0/(298.15 - 228.) ) * (228. - m_temp)/228.
- m_temp / (228.*228.) * log( (298.15*(m_temp-228.)) / (m_temp*(298.15-228.)) ));
double a3term = m_a3 / (m_temp - 228.) * (pbar - m_presR_bar);
double a4term = m_a4 / (m_temp - 228.) * log((2600. + pbar)/(2600. + m_presR_bar));
double nu = 166027;
double r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
double gval = gstar(m_temp, m_pres, 0);
double r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
double omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval) );
double relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
double Z = -1.0 / relepsilon;
double wterm = - omega_j * (Z + 1.0);
double wrterm = m_omega_pr_tr * (m_Z_pr_tr + 1.0);
double yterm = m_omega_pr_tr * m_Y_pr_tr * (m_temp - 298.15);
double deltaG_calgmol = sterm + c1term + a1term + a2term + c2term + a3term + a4term + wterm + wrterm + yterm;
// Convert to Joules / kmol
double deltaG = deltaG_calgmol * 1.0E3 * 4.184;
return deltaG;
}
double PDSS_HKFT::electrostatic_radii_calc() {
return 0.0;
}
//! Internal formula for the calculation of a_g()
/*
* The output of this is in units of Angstroms
*/
double PDSS_HKFT::ag(const double temp, const int ifunc) const {
static double ag_coeff[3] = { -2.037662, 5.747000E-3, -6.557892E-6};
if (ifunc == 0) {
double t2 = temp * temp;
double val = ag_coeff[0] + ag_coeff[1] * temp + ag_coeff[2] * t2;
return val;
} else if (ifunc == 1) {
return ag_coeff[1] + ag_coeff[2] * 2.0 * temp;
}
if (ifunc != 2) {
return 0.0;
}
return ag_coeff[2] * 2.0;;
}
//! Internal formula for the calculation of b_g()
/*
* the output of this is unitless
*/
double PDSS_HKFT::bg(const double temp, const int ifunc) const {
static double bg_coeff[3] = { 6.107361, -1.074377E-2, 1.268348E-5};
if (ifunc == 0) {
double t2 = temp * temp;
double val = bg_coeff[0] + bg_coeff[1] * temp + bg_coeff[2] * t2;
return val;
} else if (ifunc == 1) {
return bg_coeff[1] + bg_coeff[2] * 2.0 * temp;
}
if (ifunc != 2) {
return 0.0;
}
return bg_coeff[2] * 2.0;
}
double PDSS_HKFT::f(const double temp, const double pres, const int ifunc) const {
static double af_coeff[3] = { 3.666666E1, -0.1504956E-9, 0.5107997E-13};
double TC = temp - 273.15;
double presBar = pres / 1.0E5;
if (TC < 155.0) return 0.0;
if (TC > 355.0) TC = 355.0;
if (presBar > 1000.) return 0.0;
double T1 = (TC-155.0)/300.;
double fac1;
double p2 = presBar * presBar;
double p3 = presBar * p2;
double p4 = p2 * p2;
double fac2 = af_coeff[1] * p3 + af_coeff[2] * p4;
if (ifunc == 0) {
fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
return fac1 * fac2;
} else if (ifunc == 1) {
fac1 = (4.8 * pow(T1,3.8) + 16.0 * af_coeff[0] * pow(T1, 15.0)) / 300.;
return fac1 * fac2;
} else if (ifunc == 2) {
fac1 = (4.8 * 3.8 * pow(T1,2.8) + 16.0 * 15.0 * af_coeff[0] * pow(T1, 14.0)) / (300. * 300.);
return fac1 * fac2;
} else if (ifunc == 3) {
fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
fac2 = (3.0 * af_coeff[1] * p2 + 4.0 * af_coeff[2] * p3 )/ 1.0E5;
return fac1 * fac2;
} else {
throw CanteraError("HKFT_PDSS::gg", "unimplemented");
}
return 0.0;
}
double PDSS_HKFT::g(const double temp, const double pres, const int ifunc) const {
double afunc = ag(temp, 0);
double bfunc = bg(temp, 0);
m_waterSS->setState_TP(temp, pres);
m_densWaterSS = m_waterSS->density();
// density in gm cm-3
double dens = m_densWaterSS * 1.0E-3;
double gval = afunc * pow((1.0-dens), bfunc);
if (dens >= 1.0) {
return 0.0;
}
if (ifunc == 0) {
return gval;
} else if (ifunc == 1 || ifunc == 2) {
double afuncdT = ag(temp, 1);
double bfuncdT = bg(temp, 1);
double alpha = m_waterSS->thermalExpansionCoeff();
double fac1 = afuncdT * gval / afunc;
double fac2 = bfuncdT * gval * log(1.0 - dens);
double fac3 = gval * alpha * bfunc * dens / (1.0 - dens);
double dgdt = fac1 + fac2 + fac3;
if (ifunc == 1) {
return dgdt;
}
double afuncdT2 = ag(temp, 2);
double bfuncdT2 = bg(temp, 2);
double dfac1dT = dgdt * afuncdT / afunc + afuncdT2 * gval / afunc
- afuncdT * afuncdT * gval / (afunc * afunc);
double ddensdT = - alpha * dens;
double dfac2dT = bfuncdT2 * gval * log(1.0 - dens)
+ bfuncdT * dgdt * log(1.0 - dens)
- bfuncdT * gval /(1.0 - dens) * ddensdT;
double dalphadT = m_waterSS->dthermalExpansionCoeffdT();
double dfac3dT = dgdt * alpha * bfunc * dens / (1.0 - dens)
+ gval * dalphadT * bfunc * dens / (1.0 - dens)
+ gval * alpha * bfuncdT * dens / (1.0 - dens)
+ gval * alpha * bfunc * ddensdT / (1.0 - dens)
- gval * alpha * bfunc * dens / ((1.0 - dens) * (1.0 - dens)) * ddensdT;
return dfac1dT + dfac2dT + dfac3dT;
} else {
throw CanteraError("HKFT_PDSS::gg", "unimplemented");
}
return 0.0;
}
double PDSS_HKFT::gstar(const double temp, const double pres, const int ifunc) const {
double gval = g(temp, pres, ifunc);
double fval = f(temp, pres, ifunc);
return gval - fval;
}
/* awData structure */
/**
* Database for atomic molecular weights
*
* Values are taken from the 1989 Standard Atomic Weights, CRC
*
* awTable[] is a static function with scope limited to this file.
* It can only be referenced via the static Elements class function,
* LookupWtElements().
*
* units = kg / kg-mol (or equivalently gm / gm-mol)
*
* (note: this structure was picked because it's simple, compact,
* and extensible).
*
*/
struct GeData {
char name[4]; ///< Null Terminated name, First letter capitalized
double GeValue; /// < Gibbs free energies of elements J kmol-1
};
//! Values of G_elements(T=298.15,1atm)
/*!
* all units are Joules kmol-1
*/
static struct GeData geDataTable[] = {
{"H", -19.48112E6}, // NIST Webbook - Cox, Wagman 1984
{"Na", -15.29509E6}, // NIST Webbook - Cox, Wagman 1984
{"O", -30.58303E6}, // NIST Webbook - Cox, Wagman 1984
{"Cl", -33.25580E6}, // NIST Webbook - Cox, Wagman 1984
{"Si", -5.61118E6}, // Janaf
{"C", -1.71138E6}, // barin, Knack, NBS Bulletin 1971
{"S", -9.55690E6}, // Yellow - webbook
{"Al", -8.42870E6}, // Webbook polynomial
{"K", -19.26943E6}, // Webbook
{"E", 0.0} // Don't overcount
};
//! Static function to look up Element Free Energies
/*!
*
* This static function looks up the argument string in the
* database above and returns the associated Gibbs Free energies.
*
* @param ElemName String. Only the first 3 characters are significant
*
* @return
* Return value contains the Gibbs free energy for that element
*
* @exception CanteraError
* If a match is not found, a CanteraError is thrown as well
*/
double PDSS_HKFT::LookupGe(const std::string& s) {
int num = sizeof(geDataTable) / sizeof(struct GeData);
string s3 = s.substr(0,3);
for (int i = 0; i < num; i++) {
//if (!std::strncmp(s.c_str(), aWTable[i].name, 3)) {
if (s3 == geDataTable[i].name) {
return (geDataTable[i].GeValue);
}
}
throw CanteraError("LookupGe", "element " + s + " not found");
return -1.0;
}
void PDSS_HKFT::convertDGFormation() {
/*
* Ok let's get the element compositions and conversion factors.
*/
int ne = m_tp->nElements();
double na;
double ge;
string ename;
double totalSum = 0.0;
for (int m = 0; m < ne; m++) {
na = m_tp->nAtoms(m_spindex, m);
if (na > 0.0) {
ename = m_tp->elementName(m);
ge = LookupGe(ename);
totalSum += na * ge;
}
}
// Add in the charge
if (m_charge_j != 0.0) {
ename = "H";
ge = LookupGe(ename);
totalSum -= m_charge_j * ge;
}
// Ok, now do the calculation. Convert to joules kmol-1
double dg = m_deltaG_formation_tr_pr * 4.184 * 1.0E3;
//! Store the result into an internal variable.
m_Mu0_tr_pr = dg + totalSum;
}
}

122
Cantera/src/thermo/PDSS_HKFT.h
View file

@ -24,6 +24,8 @@ class WaterPropsIAPWS;
namespace Cantera {
class XML_Node;
class VPStandardState;
class PDSS_Water;
class WaterProps;
//! Class for pressure dependent standard states corresponding to
@ -403,8 +405,126 @@ namespace Cantera {
//@}
protected:
private:
//! Main routine that actually calculates the gibbs free energy difference
//! between the reference state at Tr, Pr and T,P
/*!
* This is eEqn. 59 in Johnson et al. (1992).
*
*/
double deltaG() const;
double electrostatic_radii_calc();
double ag(const double temp, const int ifunc = 0) const;
double bg(const double temp, const int ifunc = 0) const;
double g(const double temp, const double pres, const int ifunc = 0) const;
double f(const double temp, const double pres, const int ifunc = 0) const;
double gstar(const double temp, const double pres, const int ifunc = 0) const;
double LookupGe(const std::string& s);
void convertDGFormation();
private:
//! Water standard state calculator
/*!
* derived from the equation of state for water.
*/
PDSS_Water *m_waterSS;
//! Current value of the pressure for this object
doublereal m_pres;
//! density of standard-state water
/*!
* internal temporary variable
*/
mutable double m_densWaterSS;
/**
* Pointer to the water property calculator
*/
WaterProps *m_waterProps;
//! Born coefficient for the current ion or species
doublereal m_born_coeff_j;
//! Electrostatic radii
doublereal m_r_e_j;
//! Value of deltaG of Formation at Tr and Pr (cal gmol-1)
/*!
* Tr = 298.15 Pr = 1 atm
*
* This is the delta G for the formation reaction of the
* ion from elements in their stable state at Tr, Pr.
*/
doublereal m_deltaG_formation_tr_pr;
//! Value of deltaH of Formation at Tr and Pr (cal gmol-1)
/*!
* Tr = 298.15 Pr = 1 atm
*
* This is the delta H for the formation reaction of the
* ion from elements in their stable state at Tr, Pr.
*/
doublereal m_deltaH_formation_tr_pr;
//! Value of the Absolute Gibbs Free Energy NIST scale at tr and pr
/*!
* this is the NIST scale value of Gibbs free energy at T_r = 298.15
* and P_r = 1 atm.
*
* J kmol-1
*/
doublereal m_Mu0_tr_pr;
//! Value of S_j at Tr and Pr (cal gmol-1 K-1)
/*!
* Tr = 298.15 Pr = 1 atm
*/
doublereal m_Entrop_tr_pr;
//! a1 coefficient (cal gmol-1 bar-1)
doublereal m_a1;
//! a2 coefficient (cal gmol-1)
doublereal m_a2;
//! a3 coefficient (cal K gmol-1 bar-1)
doublereal m_a3;
//! a4 coefficient (cal K gmol-1)
doublereal m_a4;
//! c1 coefficient (cal gmol-1 K-1)
doublereal m_c1;
//! c2 coefficient (cal K gmol-1)
doublereal m_c2;
//! omega_pr_tr coefficient(cal gmol-1)
doublereal m_omega_pr_tr;
//! y = dZdT = 1/(esp*esp) desp/dT
double m_Y_pr_tr;
//double m_Y_pr_tr = -5.799E-5;
double m_Z_pr_tr;
//double m_Z_pr_tr = -0.0127803;
//! Reference pressure is 1 atm in units of bar= 1.0132
doublereal m_presR_bar;
//! Charge of the ion
doublereal m_charge_j;
};

22
Cantera/src/thermo/PDSS_Water.cpp
View file

@ -428,6 +428,28 @@ namespace Cantera {
m_pres = p;
}
// Return the volumetric thermal expansion coefficient. Units: 1/K.
/*
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
doublereal PDSS_Water::thermalExpansionCoeff() const {
doublereal pres = pressure();
doublereal val = m_sub->coeffThermExp(m_temp, pres);
return val;
}
doublereal PDSS_Water::dthermalExpansionCoeffdT() const {
doublereal pres = pressure();
double tt = m_temp - 0.04;
doublereal vald = m_sub->coeffThermExp(tt, pres);
doublereal val2 = m_sub->coeffThermExp(m_temp, pres);
doublereal val = (val2 - vald) / 0.04;
return val;
}
/// critical temperature
doublereal PDSS_Water::critTemperature() const { return m_sub->Tcrit(); }

18
Cantera/src/thermo/PDSS_Water.h
View file

@ -317,6 +317,24 @@ namespace Cantera {
*/
doublereal density() const;
//! Return the volumetric thermal expansion coefficient. Units: 1/K.
/*!
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
virtual doublereal thermalExpansionCoeff() const;
//! Return the derivative of the volumetric thermal expansion coefficient. Units: 1/K2.
/*!
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
virtual doublereal dthermalExpansionCoeffdT() const;
/**
* @}
* @name Miscellaneous properties of the standard state

673
Cantera/src/thermo/SpeciesThermoFactory.cpp
View file

@ -61,8 +61,8 @@ namespace Cantera {
* @todo Make sure that spDadta_node is species Data XML node by checking its name is speciesData
*/
static void getSpeciesThermoTypes(XML_Node* spData_node,
int& has_nasa, int& has_shomate, int& has_simple,
int &has_other) {
int& has_nasa, int& has_shomate, int& has_simple,
int &has_other) {
const XML_Node& sparray = *spData_node;
std::vector<XML_Node*> sp;
@ -110,7 +110,7 @@ namespace Cantera {
return new GeneralSpeciesThermo();
}
return newSpeciesThermo(NASA*inasa
+ SHOMATE*ishomate + SIMPLE*isimple);
+ SHOMATE*ishomate + SIMPLE*isimple);
}
SpeciesThermo* SpeciesThermoFactory::
@ -142,17 +142,17 @@ namespace Cantera {
int inasa = 0, ishomate = 0, isimple = 0, iother = 0;
for (int j = 0; j < n; j++) {
try {
getSpeciesThermoTypes(nodes[j], inasa, ishomate, isimple, iother);
getSpeciesThermoTypes(nodes[j], inasa, ishomate, isimple, iother);
} catch (UnknownSpeciesThermoModel) {
iother = 1;
popError();
iother = 1;
popError();
}
}
if (iother) {
return new GeneralSpeciesThermo();
}
return newSpeciesThermo(NASA*inasa
+ SHOMATE*ishomate + SIMPLE*isimple);
+ SHOMATE*ishomate + SIMPLE*isimple);
}
@ -174,7 +174,7 @@ namespace Cantera {
return new SpeciesThermoDuo<ShomateThermo, SimpleThermo>;
default:
throw UnknownSpeciesThermo(
"SpeciesThermoFactory::newSpeciesThermo",type);
"SpeciesThermoFactory::newSpeciesThermo",type);
return 0;
}
}
@ -185,359 +185,359 @@ namespace Cantera {
* temperature.
*/
void NasaThermo::checkContinuity(std::string name, double tmid, const doublereal* clow,
doublereal* chigh) {
doublereal* chigh) {
// heat capacity
doublereal cplow = poly4(tmid, clow);
doublereal cphigh = poly4(tmid, chigh);
doublereal delta = cplow - cphigh;
if (fabs(delta/(fabs(cplow)+1.0E-4)) > 0.001) {
writelog("\n\n**** WARNING ****\nFor species "+name+
", discontinuity in cp/R detected at Tmid = "
+fp2str(tmid)+"\n");
writelog("\tValue computed using low-temperature polynomial: "
+fp2str(cplow)+".\n");
writelog("\tValue computed using high-temperature polynomial: "
+fp2str(cphigh)+".\n");
}
// enthalpy
doublereal hrtlow = enthalpy_RT(tmid, clow);
doublereal hrthigh = enthalpy_RT(tmid, chigh);
delta = hrtlow - hrthigh;
if (fabs(delta/(fabs(hrtlow)+cplow*tmid)) > 0.001) {
writelog("\n\n**** WARNING ****\nFor species "+name+
", discontinuity in h/RT detected at Tmid = "
+fp2str(tmid)+"\n");
writelog("\tValue computed using low-temperature polynomial: "
+fp2str(hrtlow)+".\n");
writelog("\tValue computed using high-temperature polynomial: "
+fp2str(hrthigh)+".\n");
}
// entropy
doublereal srlow = entropy_R(tmid, clow);
doublereal srhigh = entropy_R(tmid, chigh);
delta = srlow - srhigh;
if (fabs(delta/(fabs(srlow)+cplow)) > 0.001) {
writelog("\n\n**** WARNING ****\nFor species "+name+
", discontinuity in s/R detected at Tmid = "
+fp2str(tmid)+"\n");
writelog("\tValue computed using low-temperature polynomial: "
+fp2str(srlow)+".\n");
writelog("\tValue computed using high-temperature polynomial: "
+fp2str(srhigh)+".\n");
}
// heat capacity
doublereal cplow = poly4(tmid, clow);
doublereal cphigh = poly4(tmid, chigh);
doublereal delta = cplow - cphigh;
if (fabs(delta/(fabs(cplow)+1.0E-4)) > 0.001) {
writelog("\n\n**** WARNING ****\nFor species "+name+
", discontinuity in cp/R detected at Tmid = "
+fp2str(tmid)+"\n");
writelog("\tValue computed using low-temperature polynomial: "
+fp2str(cplow)+".\n");
writelog("\tValue computed using high-temperature polynomial: "
+fp2str(cphigh)+".\n");
}
/**
* Install a NASA polynomial thermodynamic property
* parameterization for species k into a SpeciesThermo instance.
* This is called by method installThermoForSpecies if a NASA
* block is found in the XML input.
*/
static void installNasaThermoFromXML(std::string speciesName,
SpeciesThermo& sp, int k,
const XML_Node* f0ptr, const XML_Node* f1ptr) {
doublereal tmin0, tmax0, tmin1, tmax1, tmin, tmid, tmax;
const XML_Node& f0 = *f0ptr;
// default to a single temperature range
bool dualRange = false;
// but if f1ptr is suppled, then it is a two-range
// parameterization
if (f1ptr) {dualRange = true;}
tmin0 = fpValue(f0["Tmin"]);
tmax0 = fpValue(f0["Tmax"]);
tmin1 = tmax0;
tmax1 = tmin1 + 0.0001;
if (dualRange) {
tmin1 = fpValue((*f1ptr)["Tmin"]);
tmax1 = fpValue((*f1ptr)["Tmax"]);
}
vector_fp c0, c1;
if (fabs(tmax0 - tmin1) < 0.01) {
// f0 has the lower T data, and f1 the higher T data
tmin = tmin0;
tmid = tmax0;
tmax = tmax1;
getFloatArray(f0.child("floatArray"), c0, false);
if (dualRange)
getFloatArray(f1ptr->child("floatArray"), c1, false);
else {
// if there is no higher range data, then copy c0 to c1.
c1.resize(7,0.0);
copy(c0.begin(), c0.end(), c1.begin());
}
}
else if (fabs(tmax1 - tmin0) < 0.01) {
// f1 has the lower T data, and f0 the higher T data
tmin = tmin1;
tmid = tmax1;
tmax = tmax0;
getFloatArray(f1ptr->child("floatArray"), c0, false);
getFloatArray(f0.child("floatArray"), c1, false);
}
else {
throw CanteraError("installNasaThermo",
"non-continuous temperature ranges.");
}
// The NasaThermo species property manager expects the
// coefficients in a different order, so rearrange them.
array_fp c(15);
c[0] = tmid;
doublereal p0 = OneAtm;
c[1] = c0[5];
c[2] = c0[6];
copy(c0.begin(), c0.begin()+5, c.begin() + 3);
c[8] = c1[5];
c[9] = c1[6];
copy(c1.begin(), c1.begin()+5, c.begin() + 10);
sp.install(speciesName, k, NASA, &c[0], tmin, tmax, p0);
// enthalpy
doublereal hrtlow = enthalpy_RT(tmid, clow);
doublereal hrthigh = enthalpy_RT(tmid, chigh);
delta = hrtlow - hrthigh;
if (fabs(delta/(fabs(hrtlow)+cplow*tmid)) > 0.001) {
writelog("\n\n**** WARNING ****\nFor species "+name+
", discontinuity in h/RT detected at Tmid = "
+fp2str(tmid)+"\n");
writelog("\tValue computed using low-temperature polynomial: "
+fp2str(hrtlow)+".\n");
writelog("\tValue computed using high-temperature polynomial: "
+fp2str(hrthigh)+".\n");
}
// entropy
doublereal srlow = entropy_R(tmid, clow);
doublereal srhigh = entropy_R(tmid, chigh);
delta = srlow - srhigh;
if (fabs(delta/(fabs(srlow)+cplow)) > 0.001) {
writelog("\n\n**** WARNING ****\nFor species "+name+
", discontinuity in s/R detected at Tmid = "
+fp2str(tmid)+"\n");
writelog("\tValue computed using low-temperature polynomial: "
+fp2str(srlow)+".\n");
writelog("\tValue computed using high-temperature polynomial: "
+fp2str(srhigh)+".\n");
}
}
/**
* Install a NASA polynomial thermodynamic property
* parameterization for species k into a SpeciesThermo instance.
* This is called by method installThermoForSpecies if a NASA
* block is found in the XML input.
*/
static void installNasaThermoFromXML(std::string speciesName,
SpeciesThermo& sp, int k,
const XML_Node* f0ptr, const XML_Node* f1ptr) {
doublereal tmin0, tmax0, tmin1, tmax1, tmin, tmid, tmax;
const XML_Node& f0 = *f0ptr;
// default to a single temperature range
bool dualRange = false;
// but if f1ptr is suppled, then it is a two-range
// parameterization
if (f1ptr) {dualRange = true;}
tmin0 = fpValue(f0["Tmin"]);
tmax0 = fpValue(f0["Tmax"]);
tmin1 = tmax0;
tmax1 = tmin1 + 0.0001;
if (dualRange) {
tmin1 = fpValue((*f1ptr)["Tmin"]);
tmax1 = fpValue((*f1ptr)["Tmax"]);
}
vector_fp c0, c1;
if (fabs(tmax0 - tmin1) < 0.01) {
// f0 has the lower T data, and f1 the higher T data
tmin = tmin0;
tmid = tmax0;
tmax = tmax1;
getFloatArray(f0.child("floatArray"), c0, false);
if (dualRange)
getFloatArray(f1ptr->child("floatArray"), c1, false);
else {
// if there is no higher range data, then copy c0 to c1.
c1.resize(7,0.0);
copy(c0.begin(), c0.end(), c1.begin());
}
}
else if (fabs(tmax1 - tmin0) < 0.01) {
// f1 has the lower T data, and f0 the higher T data
tmin = tmin1;
tmid = tmax1;
tmax = tmax0;
getFloatArray(f1ptr->child("floatArray"), c0, false);
getFloatArray(f0.child("floatArray"), c1, false);
}
else {
throw CanteraError("installNasaThermo",
"non-continuous temperature ranges.");
}
// The NasaThermo species property manager expects the
// coefficients in a different order, so rearrange them.
array_fp c(15);
c[0] = tmid;
doublereal p0 = OneAtm;
c[1] = c0[5];
c[2] = c0[6];
copy(c0.begin(), c0.begin()+5, c.begin() + 3);
c[8] = c1[5];
c[9] = c1[6];
copy(c1.begin(), c1.begin()+5, c.begin() + 10);
sp.install(speciesName, k, NASA, &c[0], tmin, tmax, p0);
}
#ifdef INCL_NASA96
/**
* Install a NASA96 polynomial thermodynamic property
* parameterization for species k into a SpeciesThermo instance.
*/
static void installNasa96ThermoFromXML(std::string speciesName,
SpeciesThermo& sp, int k,
const XML_Node* f0ptr, const XML_Node* f1ptr) {
doublereal tmin0, tmax0, tmin1, tmax1, tmin, tmid, tmax;
/**
* Install a NASA96 polynomial thermodynamic property
* parameterization for species k into a SpeciesThermo instance.
*/
static void installNasa96ThermoFromXML(std::string speciesName,
SpeciesThermo& sp, int k,
const XML_Node* f0ptr, const XML_Node* f1ptr) {
doublereal tmin0, tmax0, tmin1, tmax1, tmin, tmid, tmax;
const XML_Node& f0 = *f0ptr;
bool dualRange = false;
if (f1ptr) {dualRange = true;}
tmin0 = fpValue(f0["Tmin"]);
tmax0 = fpValue(f0["Tmax"]);
tmin1 = tmax0;
tmax1 = tmin1 + 0.0001;
if (dualRange) {
tmin1 = fpValue((*f1ptr)["Tmin"]);
tmax1 = fpValue((*f1ptr)["Tmax"]);
}
vector_fp c0, c1;
if (fabs(tmax0 - tmin1) < 0.01) {
tmin = tmin0;
tmid = tmax0;
tmax = tmax1;
getFloatArray(f0.child("floatArray"), c0, false);
if (dualRange)
getFloatArray(f1ptr->child("floatArray"), c1, false);
else {
c1.resize(7,0.0);
copy(c0.begin(), c0.end(), c1.begin());
}
}
else if (fabs(tmax1 - tmin0) < 0.01) {
tmin = tmin1;
tmid = tmax1;
tmax = tmax0;
getFloatArray(f1ptr->child("floatArray"), c0, false);
getFloatArray(f0.child("floatArray"), c1, false);
}
else {
throw CanteraError("installNasaThermo",
"non-continuous temperature ranges.");
}
array_fp c(15);
c[0] = tmid;
doublereal p0 = OneAtm;
c[1] = c0[5];
c[2] = c0[6];
copy(c0.begin(), c0.begin()+5, c.begin() + 3);
c[8] = c1[5];
c[9] = c1[6];
copy(c1.begin(), c1.begin()+5, c.begin() + 10);
sp.install(speciesName, k, NASA, &c[0], tmin, tmax, p0);
const XML_Node& f0 = *f0ptr;
bool dualRange = false;
if (f1ptr) {dualRange = true;}
tmin0 = fpValue(f0["Tmin"]);
tmax0 = fpValue(f0["Tmax"]);
tmin1 = tmax0;
tmax1 = tmin1 + 0.0001;
if (dualRange) {
tmin1 = fpValue((*f1ptr)["Tmin"]);
tmax1 = fpValue((*f1ptr)["Tmax"]);
}
vector_fp c0, c1;
if (fabs(tmax0 - tmin1) < 0.01) {
tmin = tmin0;
tmid = tmax0;
tmax = tmax1;
getFloatArray(f0.child("floatArray"), c0, false);
if (dualRange)
getFloatArray(f1ptr->child("floatArray"), c1, false);
else {
c1.resize(7,0.0);
copy(c0.begin(), c0.end(), c1.begin());
}
}
else if (fabs(tmax1 - tmin0) < 0.01) {
tmin = tmin1;
tmid = tmax1;
tmax = tmax0;
getFloatArray(f1ptr->child("floatArray"), c0, false);
getFloatArray(f0.child("floatArray"), c1, false);
}
else {
throw CanteraError("installNasaThermo",
"non-continuous temperature ranges.");
}
array_fp c(15);
c[0] = tmid;
doublereal p0 = OneAtm;
c[1] = c0[5];
c[2] = c0[6];
copy(c0.begin(), c0.begin()+5, c.begin() + 3);
c[8] = c1[5];
c[9] = c1[6];
copy(c1.begin(), c1.begin()+5, c.begin() + 10);
sp.install(speciesName, k, NASA, &c[0], tmin, tmax, p0);
}
#endif
/**
* Install a Shomate polynomial thermodynamic property
* parameterization for species k.
*/
static void installShomateThermoFromXML(std::string speciesName,
SpeciesThermo& sp, int k,
const XML_Node* f0ptr, const XML_Node* f1ptr) {
doublereal tmin0, tmax0, tmin1, tmax1, tmin, tmid, tmax;
const XML_Node& f0 = *f0ptr;
bool dualRange = false;
if (f1ptr) {dualRange = true;}
tmin0 = fpValue(f0["Tmin"]);
tmax0 = fpValue(f0["Tmax"]);
tmin1 = tmax0;
tmax1 = tmin1 + 0.0001;
if (dualRange) {
tmin1 = fpValue((*f1ptr)["Tmin"]);
tmax1 = fpValue((*f1ptr)["Tmax"]);
}
vector_fp c0, c1;
if (fabs(tmax0 - tmin1) < 0.01) {
tmin = tmin0;
tmid = tmax0;
tmax = tmax1;
getFloatArray(f0.child("floatArray"), c0, false);
if (dualRange)
getFloatArray(f1ptr->child("floatArray"), c1, false);
else {
c1.resize(7,0.0);
copy(c0.begin(), c0.begin()+7, c1.begin());
}
}
else if (fabs(tmax1 - tmin0) < 0.01) {
tmin = tmin1;
tmid = tmax1;
tmax = tmax0;
getFloatArray(f1ptr->child("floatArray"), c0, false);
getFloatArray(f0.child("floatArray"), c1, false);
}
else {
throw CanteraError("installShomateThermoFromXML",
"non-continuous temperature ranges.");
}
array_fp c(15);
c[0] = tmid;
doublereal p0 = OneAtm;
copy(c0.begin(), c0.begin()+7, c.begin() + 1);
copy(c1.begin(), c1.begin()+7, c.begin() + 8);
sp.install(speciesName, k, SHOMATE, &c[0], tmin, tmax, p0);
}
/**
* Install a constant-cp thermodynamic property
* parameterization for species k.
*/
static void installSimpleThermoFromXML(std::string speciesName,
SpeciesThermo& sp, int k,
const XML_Node& f) {
doublereal tmin, tmax;
tmin = fpValue(f["Tmin"]);
tmax = fpValue(f["Tmax"]);
if (tmax == 0.0) tmax = 1.0e30;
vector_fp c(4);
c[0] = getFloat(f, "t0", "-");
c[1] = getFloat(f, "h0", "-");
c[2] = getFloat(f, "s0", "-");
c[3] = getFloat(f, "cp0", "-");
doublereal p0 = OneAtm;
sp.install(speciesName, k, SIMPLE, &c[0], tmin, tmax, p0);
}
/**
* Install a NASA9 polynomial thermodynamic property
* parameterization for species k into a SpeciesThermo instance.
* This is called by method installThermoForSpecies if a NASA9
* block is found in the XML input.
*/
static void installNasa9ThermoFromXML(std::string speciesName,
/**
* Install a Shomate polynomial thermodynamic property
* parameterization for species k.
*/
static void installShomateThermoFromXML(std::string speciesName,
SpeciesThermo& sp, int k,
const std::vector<XML_Node*>& tp)
{
const XML_Node * fptr = tp[0];
int nRegTmp = tp.size();
int nRegions = 0;
vector_fp cPoly;
Nasa9Poly1 *np_ptr = 0;
std::vector<Nasa9Poly1 *> regionPtrs;
doublereal tmin, tmax, pref;
// Loop over all of the possible temperature regions
for (int i = 0; i < nRegTmp; i++) {
fptr = tp[i];
if (fptr) {
if (fptr->name() == "NASA9") {
if (fptr->hasChild("floatArray")) {
const XML_Node* f0ptr, const XML_Node* f1ptr) {
doublereal tmin0, tmax0, tmin1, tmax1, tmin, tmid, tmax;
tmin = fpValue((*fptr)["Tmin"]);
tmax = fpValue((*fptr)["Tmax"]);
pref = fpValue((*fptr)["P0"]);
const XML_Node& f0 = *f0ptr;
bool dualRange = false;
if (f1ptr) {dualRange = true;}
tmin0 = fpValue(f0["Tmin"]);
tmax0 = fpValue(f0["Tmax"]);
tmin1 = tmax0;
tmax1 = tmin1 + 0.0001;
if (dualRange) {
tmin1 = fpValue((*f1ptr)["Tmin"]);
tmax1 = fpValue((*f1ptr)["Tmax"]);
}
getFloatArray(fptr->child("floatArray"), cPoly, false);
if (cPoly.size() != 9) {
throw CanteraError("installNasa9ThermoFromXML",
"Expected 9 coeff polynomial");
}
np_ptr = new Nasa9Poly1(k, tmin, tmax, pref,
DATA_PTR(cPoly));
regionPtrs.push_back(np_ptr);
nRegions++;
}
}
vector_fp c0, c1;
if (fabs(tmax0 - tmin1) < 0.01) {
tmin = tmin0;
tmid = tmax0;
tmax = tmax1;
getFloatArray(f0.child("floatArray"), c0, false);
if (dualRange)
getFloatArray(f1ptr->child("floatArray"), c1, false);
else {
c1.resize(7,0.0);
copy(c0.begin(), c0.begin()+7, c1.begin());
}
}
else if (fabs(tmax1 - tmin0) < 0.01) {
tmin = tmin1;
tmid = tmax1;
tmax = tmax0;
getFloatArray(f1ptr->child("floatArray"), c0, false);
getFloatArray(f0.child("floatArray"), c1, false);
}
else {
throw CanteraError("installShomateThermoFromXML",
"non-continuous temperature ranges.");
}
array_fp c(15);
c[0] = tmid;
doublereal p0 = OneAtm;
copy(c0.begin(), c0.begin()+7, c.begin() + 1);
copy(c1.begin(), c1.begin()+7, c.begin() + 8);
sp.install(speciesName, k, SHOMATE, &c[0], tmin, tmax, p0);
}
/**
* Install a constant-cp thermodynamic property
* parameterization for species k.
*/
static void installSimpleThermoFromXML(std::string speciesName,
SpeciesThermo& sp, int k,
const XML_Node& f) {
doublereal tmin, tmax;
tmin = fpValue(f["Tmin"]);
tmax = fpValue(f["Tmax"]);
if (tmax == 0.0) tmax = 1.0e30;
vector_fp c(4);
c[0] = getFloat(f, "t0", "-");
c[1] = getFloat(f, "h0", "-");
c[2] = getFloat(f, "s0", "-");
c[3] = getFloat(f, "cp0", "-");
doublereal p0 = OneAtm;
sp.install(speciesName, k, SIMPLE, &c[0], tmin, tmax, p0);
}
/**
* Install a NASA9 polynomial thermodynamic property
* parameterization for species k into a SpeciesThermo instance.
* This is called by method installThermoForSpecies if a NASA9
* block is found in the XML input.
*/
static void installNasa9ThermoFromXML(std::string speciesName,
SpeciesThermo& sp, int k,
const std::vector<XML_Node*>& tp)
{
const XML_Node * fptr = tp[0];
int nRegTmp = tp.size();
int nRegions = 0;
vector_fp cPoly;
Nasa9Poly1 *np_ptr = 0;
std::vector<Nasa9Poly1 *> regionPtrs;
doublereal tmin, tmax, pref;
// Loop over all of the possible temperature regions
for (int i = 0; i < nRegTmp; i++) {
fptr = tp[i];
if (fptr) {
if (fptr->name() == "NASA9") {
if (fptr->hasChild("floatArray")) {
tmin = fpValue((*fptr)["Tmin"]);
tmax = fpValue((*fptr)["Tmax"]);
pref = fpValue((*fptr)["P0"]);
getFloatArray(fptr->child("floatArray"), cPoly, false);
if (cPoly.size() != 9) {
throw CanteraError("installNasa9ThermoFromXML",
"Expected 9 coeff polynomial");
}
np_ptr = new Nasa9Poly1(k, tmin, tmax, pref,
DATA_PTR(cPoly));
regionPtrs.push_back(np_ptr);
nRegions++;
}
}
}
if (nRegions == 0) {
throw UnknownSpeciesThermoModel("installThermoForSpecies",
speciesName, " " );
} else if (nRegions == 1) {
sp.install_STIT(np_ptr);
} else {
Nasa9PolyMultiTempRegion* npMulti_ptr = new Nasa9PolyMultiTempRegion(regionPtrs);
sp.install_STIT(npMulti_ptr);
}
}
if (nRegions == 0) {
throw UnknownSpeciesThermoModel("installThermoForSpecies",
speciesName, " " );
} else if (nRegions == 1) {
sp.install_STIT(np_ptr);
} else {
Nasa9PolyMultiTempRegion* npMulti_ptr = new Nasa9PolyMultiTempRegion(regionPtrs);
sp.install_STIT(npMulti_ptr);
}
}
/**
* Install an Adsorbate thermodynamic property
* parameterization for species k into a SpeciesThermo instance.
* This is called by method installThermoForSpecies if a NASA9
* block is found in the XML input.
*/
/**
* Install an Adsorbate thermodynamic property
* parameterization for species k into a SpeciesThermo instance.
* This is called by method installThermoForSpecies if a NASA9
* block is found in the XML input.
*/
#ifdef WITH_ADSORBATE
static void installAdsorbateThermoFromXML(std::string speciesName,
SpeciesThermo& sp, int k,
const XML_Node& f) {
vector_fp freqs;
doublereal tmin, tmax, pref;
int nfreq = 0;
tmin = fpValue(f["Tmin"]);
tmax = fpValue(f["Tmax"]);
pref = fpValue(f["P0"]);
if (tmax == 0.0) tmax = 1.0e30;
static void installAdsorbateThermoFromXML(std::string speciesName,
SpeciesThermo& sp, int k,
const XML_Node& f) {
vector_fp freqs;
doublereal tmin, tmax, pref;
int nfreq = 0;
tmin = fpValue(f["Tmin"]);
tmax = fpValue(f["Tmax"]);
pref = fpValue(f["P0"]);
if (tmax == 0.0) tmax = 1.0e30;
if (f.hasChild("floatArray")) {
getFloatArray(f.child("floatArray"), freqs, false);
nfreq = freqs.size();
}
for (int n = 0; n < nfreq; n++) {
freqs[n] *= 3.0e10;
}
vector_fp coeffs(nfreq + 2);
coeffs[0] = nfreq;
coeffs[1] = getFloat(f, "binding_energy", "-");
copy(freqs.begin(), freqs.end(), coeffs.begin() + 2);
//posc = new Adsorbate(k, tmin, tmax, pref,
// DATA_PTR(coeffs));
(&sp)->install(speciesName, k, ADSORBATE, &coeffs[0], tmin, tmax, pref);
if (f.hasChild("floatArray")) {
getFloatArray(f.child("floatArray"), freqs, false);
nfreq = freqs.size();
}
for (int n = 0; n < nfreq; n++) {
freqs[n] *= 3.0e10;
}
vector_fp coeffs(nfreq + 2);
coeffs[0] = nfreq;
coeffs[1] = getFloat(f, "binding_energy", "-");
copy(freqs.begin(), freqs.end(), coeffs.begin() + 2);
//posc = new Adsorbate(k, tmin, tmax, pref,
// DATA_PTR(coeffs));
(&sp)->install(speciesName, k, ADSORBATE, &coeffs[0], tmin, tmax, pref);
}
#endif
/**
* Install a species thermodynamic property parameterization
* for one species into a species thermo manager.
* @param k species number
* @param s XML node specifying species
* @param spthermo species thermo manager
* @param phaseNode_ptr Optional Pointer to the XML phase
* information for the phase in which the species
* resides
*/
/**
* Install a species thermodynamic property parameterization
* for one species into a species thermo manager.
* @param k species number
* @param s XML node specifying species
* @param spthermo species thermo manager
* @param phaseNode_ptr Optional Pointer to the XML phase
* information for the phase in which the species
* resides
*/
void SpeciesThermoFactory::
installThermoForSpecies(int k, const XML_Node& s,
SpeciesThermo& spthermo,
@ -572,6 +572,9 @@ namespace Cantera {
else if (f->name() == "NASA9") {
installNasa9ThermoFromXML(s["name"], spthermo, k, tp);
}
// else if (f->name() == "HKFT") {
// installHKFTThermoFromXML(s["name"], spthermo, k, tp);
//}
#ifdef WITH_ADSORBATE
else if (f->name() == "adsorbate") {
installAdsorbateThermoFromXML(s["name"], spthermo, k, *f);

59
Cantera/src/thermo/VPSSMgrFactory.cpp
View file

@ -84,29 +84,62 @@ namespace Cantera {
sparray.getChildren("species",sp);
size_t n, ns = sp.size();
for (n = 0; n < ns; n++) {
bool ifound = false;
XML_Node* spNode = sp[n];
if (spNode->hasChild("standardState")) {
const XML_Node& ssN = sp[n]->child("standardState");
string mm = ssN["model"];
if (mm == "waterIAPWS" || mm == "waterPDSS") {
has_water++;
ifound = true;
}
} else {
if (mm == "HKFT") {
has_hptx++;
ifound = true;
}
}
if (!ifound) {
if (spNode->hasChild("thermo")) {
const XML_Node& th = sp[n]->child("thermo");
if (th.hasChild("NASA")) has_nasa = 1;
if (th.hasChild("Shomate")) has_shomate = 1;
if (th.hasChild("const_cp")) has_simple = 1;
if (th.hasChild("poly")) {
if (th.child("poly")["order"] == "1") has_simple = 1;
else throw CanteraError("newSpeciesThermo",
"poly with order > 1 not yet supported");
if (th.hasChild("NASA")) {
has_nasa++;
ifound = true;
}
if (th.hasChild("Shomate")) {
has_shomate++;
ifound = true;
}
if (th.hasChild("const_cp")){
has_simple = 1;
ifound = true;
}
if (th.hasChild("poly")) {
if (th.child("poly")["order"] == "1") {
has_simple = 1;
ifound = true;
} else throw CanteraError("newSpeciesThermo",
"poly with order > 1 not yet supported");
}
if (th.hasChild("Mu0")) {
has_other++;
ifound = true;
}
if (th.hasChild("NASA9")) {
has_other++;
ifound = true;
}
if (th.hasChild("NASA9MULTITEMP")) {
has_other++;
ifound = true;
}
if (th.hasChild("adsorbate")) {
has_other++;
ifound = true;
}
if (th.hasChild("HKFT")) {
has_hptx++;
ifound = true;
}
if (th.hasChild("Mu0")) has_other = 1;
if (th.hasChild("NASA9")) has_other = 1;
if (th.hasChild("NASA9MULTITEMP")) has_other = 1;
if (th.hasChild("adsorbate")) has_other = 1;
} else {
throw UnknownVPSSMgrModel("getVPSSMgrTypes:",
spNode->attrib("name"));

4
Cantera/src/thermo/VPSSMgr_General.cpp
View file

@ -30,6 +30,7 @@
#include "PDSS_IdealGas.h"
#include "PDSS_Water.h"
#include "PDSS_ConstVol.h"
#include "PDSS_HKFT.h"
using namespace std;
@ -133,6 +134,9 @@ namespace Cantera {
} else if (model == "waterIAPWS" || model == "waterPDSS") {
doST = false;
kPDSS = new PDSS_Water();
} else if (model == "HKFT") {
doST = false;
kPDSS = new PDSS_HKFT(m_vptp_ptr, k, speciesNode, *phaseNode_ptr, true);
} else {
throw CanteraError("VPSSMgr_General::returnPDSS_ptr",
"unknown");

63
Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp
View file

@ -28,6 +28,7 @@
#include "xml.h"
#include "VPStandardStateTP.h"
#include "PDSS_Water.h"
#include "PDSS_HKFT.h"
using namespace std;
@ -46,7 +47,7 @@ namespace Cantera {
VPSSMgr_Water_HKFT::~VPSSMgr_Water_HKFT()
{
delete m_waterSS;
// m_waterSS is owned by VPStandardState
}
VPSSMgr_Water_HKFT::VPSSMgr_Water_HKFT(const VPSSMgr_Water_HKFT &right) :
@ -167,14 +168,6 @@ namespace Cantera {
void VPSSMgr_Water_HKFT::_updateStandardStateThermo() {
doublereal RT = GasConstant * m_tlast;
doublereal del_pRT = (m_plast - m_p0) / (RT);
for (int k = 1; k < m_kk; k++) {
m_hss_RT[k] = m_h0_RT[k] + del_pRT * m_Vss[k];
m_cpss_R[k] = m_cp0_R[k];
m_sss_R[k] = m_s0_R[k];
m_gss_RT[k] = m_hss_RT[k] - m_sss_R[k];
// m_Vss[k] constant
}
// Do the water
m_waterSS->setState_TP(m_tlast, m_plast);
m_hss_RT[0] = (m_waterSS->enthalpy_mole())/ RT;
@ -182,6 +175,18 @@ namespace Cantera {
m_cpss_R[0] = (m_waterSS->cp_mole()) / GasConstant;
m_gss_RT[0] = (m_hss_RT[0] - m_sss_R[0]);
m_Vss[0] = (m_waterSS->density()) / m_vptp_ptr->molecularWeight(0);
for (int k = 1; k < m_kk; k++) {
PDSS_HKFT *ps = (PDSS_HKFT *) m_vptp_ptr->providePDSS(k);
ps->setState_TP(m_tlast, m_plast);
m_hss_RT[k] = m_h0_RT[k] + del_pRT * m_Vss[k];
m_cpss_R[k] = m_cp0_R[k];
m_sss_R[k] = m_s0_R[k];
m_gss_RT[k] = ps->gibbs_mole() / RT;
}
}
void VPSSMgr_Water_HKFT::initThermo() {
@ -214,20 +219,25 @@ namespace Cantera {
throw CanteraError("VPSSMgr_Water_HKFT::initThermoXML",
"no standardState Node for species " + s->name());
}
std::string model = (*ss)["model"];
if (model != "constant_incompressible") {
std::string model = lowercase((*ss)["model"]);
if (model != "hkft") {
throw CanteraError("VPSSMgr_Water_HKFT::initThermoXML",
"standardState model for species isn't constant_incompressible: " + s->name());
"standardState model for species isn't hkft: " + s->name());
}
m_Vss[k] = getFloat(*ss, "molarVolume", "-");
// m_Vss[k] = getFloat(*ss, "molarVolume", "-");
}
}
PDSS *
VPSSMgr_Water_HKFT::createInstallPDSS(int k, const XML_Node& speciesNode,
const XML_Node *phaseNode_ptr) {
VPSSMgr::installSTSpecies(k, speciesNode, phaseNode_ptr);
PDSS *kPDSS = 0;
const XML_Node *ss = speciesNode.findByName("standardState");
if (!ss) {
throw CanteraError("VPSSMgr_Water_HKFT::installSpecies",
"no standardState Node for species " + speciesNode.name());
}
// Will have to do something for water
// -> make sure it's species 0
// -> make sure it's designated as a real water EOS
@ -237,33 +247,28 @@ namespace Cantera {
throw CanteraError("VPSSMgr_Water_HKFT::installSpecies",
"h2o wrong name: " + xn);
}
const XML_Node *ss = speciesNode.findByName("standardState");
std::string model = (*ss)["model"];
if (model != "waterIAPSS" && model != "waterPDSS") {
if (model != "waterIAPWS" && model != "waterPDSS") {
throw CanteraError("VPSSMgr_Water_HKFT::installSpecies",
"wrong SS mode: " + model);
}
VPSSMgr::installSTSpecies(k, speciesNode, phaseNode_ptr);
if (m_waterSS) delete m_waterSS;
m_waterSS = new PDSS_Water(m_vptp_ptr, 0);
kPDSS = m_waterSS;
} else {
const XML_Node *ss = speciesNode.findByName("standardState");
if (!ss) {
throw CanteraError("VPSSMgr_Water_HKFT::installSpecies",
"no standardState Node for species " + speciesNode.name());
}
std::string model = (*ss)["model"];
if (model != "constant_incompressible") {
if (model != "HKFT") {
throw CanteraError("VPSSMgr_Water_HKFT::initThermoXML",
"standardState model for species isn't "
"constant_incompressible: " + speciesNode.name());
"HKFT: " + speciesNode.name());
}
if ((int) m_Vss.size() < k+1) {
m_Vss.resize(k+1, 0.0);
}
m_Vss[k] = getFloat(*ss, "molarVolume", "-");
kPDSS = new PDSS_HKFT(m_vptp_ptr, k, speciesNode, *phaseNode_ptr, true);
}
return 0;
return kPDSS;
}
PDSS_enumType VPSSMgr_Water_HKFT::reportPDSSType(int k) const {

465
Cantera/src/thermo/WaterPDSS.cpp
View file

@ -1,465 +0,0 @@
/**
* @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 "ThermoFactory.h"
#include <math.h>
#include "ThermoPhase.h"
namespace Cantera {
/**
* Basic list of constructors and duplicators
*/
WaterPDSS::WaterPDSS() :
PDSS(),
m_sub(0),
m_temp(0.0),
m_dens(0.0),
m_iState(-3000),
EW_Offset(0.0),
SW_Offset(0.0),
m_verbose(0),
m_allowGasPhase(false)
{
m_sub = new WaterPropsIAPWS();
m_spthermo = 0;
}
WaterPDSS::WaterPDSS(ThermoPhase *tp, int spindex) :
PDSS(tp, spindex),
m_sub(0),
m_temp(0.0),
m_dens(0.0),
m_iState(-3000),
EW_Offset(0.0),
SW_Offset(0.0),
m_verbose(0),
m_allowGasPhase(false)
{
m_sub = new WaterPropsIAPWS();
constructPDSS(tp, spindex);
m_spthermo = 0;
}
WaterPDSS::WaterPDSS(ThermoPhase *tp, int spindex,
std::string inputFile, std::string id) :
PDSS(tp, spindex),
m_sub(0),
m_temp(0.0),
m_dens(0.0),
m_iState(-3000),
EW_Offset(0.0),
SW_Offset(0.0),
m_verbose(0),
m_allowGasPhase(false)
{
m_sub = new WaterPropsIAPWS();
constructPDSSFile(tp, spindex, inputFile, id);
m_spthermo = 0;
}
WaterPDSS::WaterPDSS(ThermoPhase *tp, int spindex,
XML_Node& phaseRoot, std::string id) :
PDSS(tp, spindex),
m_sub(0),
m_temp(0.0),
m_dens(0.0),
m_iState(-3000),
EW_Offset(0.0),
SW_Offset(0.0),
m_verbose(0),
m_allowGasPhase(false)
{
m_sub = new WaterPropsIAPWS();
constructPDSSXML(tp, spindex, phaseRoot, id) ;
m_spthermo = 0;
}
WaterPDSS::WaterPDSS(const WaterPDSS &b) :
PDSS(),
m_sub(0),
m_temp(0.0),
m_dens(0.0),
m_iState(-3000),
EW_Offset(b.EW_Offset),
SW_Offset(b.SW_Offset),
m_verbose(b.m_verbose),
m_allowGasPhase(b.m_allowGasPhase)
{
m_sub = new WaterPropsIAPWS();
/*
* 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;
/*
* Call the base class operator
*/
PDSS::operator=(b);
m_sub->operator=(*(b.m_sub));
m_temp = b.m_temp;
m_dens = b.m_dens;
m_iState = b.m_iState;
EW_Offset = b.EW_Offset;
SW_Offset = b.SW_Offset;
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, std::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,
std::string inputFile, std::string id) {
if (inputFile.size() == 0) {
throw CanteraError("WaterTp::initThermo",
"input file is null");
}
std::string path = findInputFile(inputFile);
std::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, std::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) {
std::string stateString = "T = " +
fp2str(T) + " K and p = " + fp2str(p) + " Pa";
throw CanteraError("WaterPDSS:setPressure()",
"Failed to set water SS state: " + stateString);
}
m_dens = dd;
}
// Return the volumetric thermal expansion coefficient. Units: 1/K.
/*
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
doublereal WaterPDSS::thermalExpansionCoeff() const {
doublereal pres = pressure();
doublereal val = m_sub->coeffThermExp(m_temp, pres);
return val;
}
doublereal WaterPDSS::dthermalExpansionCoeffdT() const {
doublereal pres = pressure();
double tt = m_temp - 0.04;
doublereal vald = m_sub->coeffThermExp(tt, pres);
doublereal val2 = m_sub->coeffThermExp(m_temp, pres);
doublereal val = (val2 - vald) / 0.04;
return val;
}
/// critical temperature
doublereal WaterPDSS::critTemperature() const { return m_sub->Tcrit(); }
/// 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;
}
}

224
Cantera/src/thermo/WaterPDSS.h
View file

@ -1,224 +0,0 @@
/**
* @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();
WaterPDSS(ThermoPhase *tp, int spindex);
WaterPDSS(const WaterPDSS &b);
WaterPDSS& operator=(const WaterPDSS&b);
WaterPDSS(ThermoPhase *tp, int spindex,
std::string inputFile, std::string id = "");
WaterPDSS(ThermoPhase *tp, int spindex,
XML_Node& phaseRef, std::string id = "");
virtual ~WaterPDSS();
/**
*
* @name Utilities
* @{
*/
virtual int pdssType() const { return -1; }
/**
* @}
* @name 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);
//! Return the volumetric thermal expansion coefficient. Units: 1/K.
/*!
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
virtual doublereal thermalExpansionCoeff() const;
//! Return the derivative of the volumetric thermal expansion coefficient. Units: 1/K2.
/*!
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
virtual doublereal dthermalExpansionCoeffdT() const;
//@}
/// @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,
std::string inputFile, std::string id);
virtual void constructPDSSXML(ThermoPhase *tp, int spindex,
XML_Node& phaseNode, std::string id);
virtual void initThermoXML(XML_Node& eosdata, std::string id);
virtual void initThermo();
virtual void setParametersFromXML(const XML_Node& eosdata);
WaterPropsIAPWS *getWater() const {
return m_sub;
}
protected:
private:
mutable WaterPropsIAPWS *m_sub;
/**
* 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;
/**
* 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