From 34bae00f15944ceca3aca86591c3c9d5d6972ea3 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Sun, 24 Aug 2008 01:56:24 +0000 Subject: [PATCH] Interrum update for standing up HKFT standard state Gibbs free energies are working and checked against excel spread sheet! --- Cantera/src/thermo/HKFT_PDSS.cpp | 646 ------------------- Cantera/src/thermo/HKFT_PDSS.h | 260 -------- Cantera/src/thermo/HMWSoln_input.cpp | 11 +- Cantera/src/thermo/PDSS_HKFT.cpp | 514 ++++++++++++++- Cantera/src/thermo/PDSS_HKFT.h | 122 +++- Cantera/src/thermo/PDSS_Water.cpp | 22 + Cantera/src/thermo/PDSS_Water.h | 18 + Cantera/src/thermo/SpeciesThermoFactory.cpp | 673 ++++++++++---------- Cantera/src/thermo/VPSSMgrFactory.cpp | 59 +- Cantera/src/thermo/VPSSMgr_General.cpp | 4 + Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp | 63 +- Cantera/src/thermo/WaterPDSS.cpp | 465 -------------- Cantera/src/thermo/WaterPDSS.h | 224 ------- 13 files changed, 1084 insertions(+), 1997 deletions(-) delete mode 100644 Cantera/src/thermo/HKFT_PDSS.cpp delete mode 100644 Cantera/src/thermo/HKFT_PDSS.h delete mode 100644 Cantera/src/thermo/WaterPDSS.cpp delete mode 100644 Cantera/src/thermo/WaterPDSS.h diff --git a/Cantera/src/thermo/HKFT_PDSS.cpp b/Cantera/src/thermo/HKFT_PDSS.cpp deleted file mode 100644 index a038543af..000000000 --- a/Cantera/src/thermo/HKFT_PDSS.cpp +++ /dev/null @@ -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; - } - -} diff --git a/Cantera/src/thermo/HKFT_PDSS.h b/Cantera/src/thermo/HKFT_PDSS.h deleted file mode 100644 index b9c87c5d0..000000000 --- a/Cantera/src/thermo/HKFT_PDSS.h +++ /dev/null @@ -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 - - - diff --git a/Cantera/src/thermo/HMWSoln_input.cpp b/Cantera/src/thermo/HMWSoln_input.cpp index d198a3b65..c4cb26275 100644 --- a/Cantera/src/thermo/HMWSoln_input.cpp +++ b/Cantera/src/thermo/HMWSoln_input.cpp @@ -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 } } diff --git a/Cantera/src/thermo/PDSS_HKFT.cpp b/Cantera/src/thermo/PDSS_HKFT.cpp index 2cecefc21..5d6e3a44f 100644 --- a/Cantera/src/thermo/PDSS_HKFT.cpp +++ b/Cantera/src/thermo/PDSS_HKFT.cpp @@ -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; + } + + } diff --git a/Cantera/src/thermo/PDSS_HKFT.h b/Cantera/src/thermo/PDSS_HKFT.h index 136dba4c5..498f5a983 100644 --- a/Cantera/src/thermo/PDSS_HKFT.h +++ b/Cantera/src/thermo/PDSS_HKFT.h @@ -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; }; diff --git a/Cantera/src/thermo/PDSS_Water.cpp b/Cantera/src/thermo/PDSS_Water.cpp index b336f42ca..0bc11c0ed 100644 --- a/Cantera/src/thermo/PDSS_Water.cpp +++ b/Cantera/src/thermo/PDSS_Water.cpp @@ -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(); } diff --git a/Cantera/src/thermo/PDSS_Water.h b/Cantera/src/thermo/PDSS_Water.h index e229a0432..90734d398 100644 --- a/Cantera/src/thermo/PDSS_Water.h +++ b/Cantera/src/thermo/PDSS_Water.h @@ -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 diff --git a/Cantera/src/thermo/SpeciesThermoFactory.cpp b/Cantera/src/thermo/SpeciesThermoFactory.cpp index a2100229c..b490adf46 100755 --- a/Cantera/src/thermo/SpeciesThermoFactory.cpp +++ b/Cantera/src/thermo/SpeciesThermoFactory.cpp @@ -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 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; 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& tp) - { - const XML_Node * fptr = tp[0]; - int nRegTmp = tp.size(); - int nRegions = 0; - vector_fp cPoly; - Nasa9Poly1 *np_ptr = 0; - std::vector 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& tp) + { + const XML_Node * fptr = tp[0]; + int nRegTmp = tp.size(); + int nRegions = 0; + vector_fp cPoly; + Nasa9Poly1 *np_ptr = 0; + std::vector 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); diff --git a/Cantera/src/thermo/VPSSMgrFactory.cpp b/Cantera/src/thermo/VPSSMgrFactory.cpp index 9ea97c19b..d08c15970 100644 --- a/Cantera/src/thermo/VPSSMgrFactory.cpp +++ b/Cantera/src/thermo/VPSSMgrFactory.cpp @@ -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")); diff --git a/Cantera/src/thermo/VPSSMgr_General.cpp b/Cantera/src/thermo/VPSSMgr_General.cpp index f0f00ed52..9ea1ca28f 100644 --- a/Cantera/src/thermo/VPSSMgr_General.cpp +++ b/Cantera/src/thermo/VPSSMgr_General.cpp @@ -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"); diff --git a/Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp b/Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp index aa5a50b05..1bf185370 100644 --- a/Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp +++ b/Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp @@ -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 { diff --git a/Cantera/src/thermo/WaterPDSS.cpp b/Cantera/src/thermo/WaterPDSS.cpp deleted file mode 100644 index 81479a8f7..000000000 --- a/Cantera/src/thermo/WaterPDSS.cpp +++ /dev/null @@ -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 - - -#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; - } - - - -} diff --git a/Cantera/src/thermo/WaterPDSS.h b/Cantera/src/thermo/WaterPDSS.h deleted file mode 100644 index 40673013a..000000000 --- a/Cantera/src/thermo/WaterPDSS.h +++ /dev/null @@ -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 - - -