From 03cef87ff848d72a68e0300d0e910cb1e64b7234 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Fri, 22 Aug 2008 02:54:28 +0000 Subject: [PATCH] Incremental update --- Cantera/src/thermo/HKFT_PDSS.cpp | 122 +++++++++++++++++++++++++++---- Cantera/src/thermo/HKFT_PDSS.h | 38 ++++++++-- Cantera/src/thermo/PDSS.cpp | 3 +- 3 files changed, 138 insertions(+), 25 deletions(-) diff --git a/Cantera/src/thermo/HKFT_PDSS.cpp b/Cantera/src/thermo/HKFT_PDSS.cpp index e1c3f4a52..a038543af 100644 --- a/Cantera/src/thermo/HKFT_PDSS.cpp +++ b/Cantera/src/thermo/HKFT_PDSS.cpp @@ -201,10 +201,9 @@ namespace Cantera { * Calculate the Gibbs free energy in mks units of * J kmol-1 K-1. */ - doublereal HKFT_PDSS:: - gibbs_mole() const { - throw CanteraError("HKFT_PDSS::gibbs_mole()", "unimplemented"); - return (0.0); + doublereal HKFT_PDSS::gibbs_mole() const { + double val = deltaG(); + return (m_Mu0_tr_pr + val); } /** @@ -363,19 +362,19 @@ namespace Cantera { } - double HKFT_PDSS::deltaG() { + 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 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. + 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); @@ -393,9 +392,6 @@ namespace Cantera { double relepsilon = m_wprops->relEpsilon(m_temp, m_pres, 0); - //double Y_pr_tr = -5.799E-5; - //double Z_pr_tr = -0.0127803; - double Z = -1.0 / relepsilon; double wterm = - omega_j * (Z + 1.0); @@ -404,9 +400,9 @@ namespace Cantera { double yterm = m_omega_pr_tr * m_Y_pr_tr * (m_temp - 298.15); - double deltaG_calgmol = m_deltaG_tr_pr + sterm + c1term + a1term + a2term + c2term + a3term + a4term + wterm + wrterm + yterm; + double deltaG_calgmol = sterm + c1term + a1term + a2term + c2term + a3term + a4term + wterm + wrterm + yterm; - // Convert to Joules / kg + // Convert to Joules / kmol double deltaG = deltaG_calgmol * 1.0E3 * 4.184; return deltaG; } @@ -453,7 +449,7 @@ namespace Cantera { return bg_coeff[2] * 2.0; } - double HKFT_PDSS::f(const double temp, const double pres, const int ifunc) { + 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; @@ -490,7 +486,7 @@ namespace Cantera { return 0.0; } - double HKFT_PDSS::g(const double temp, const double pres, const int ifunc) { + 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); @@ -545,10 +541,106 @@ namespace Cantera { return 0.0; } - double HKFT_PDSS::gstar(const double temp, const double pres, const int ifunc) { + 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 index 0b304aeda..b9c87c5d0 100644 --- a/Cantera/src/thermo/HKFT_PDSS.h +++ b/Cantera/src/thermo/HKFT_PDSS.h @@ -138,16 +138,28 @@ namespace Cantera { 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 deltaG(); double electrostatic_radii_calc(); - private: + 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); - double f(const double temp, const double pres, const int ifunc = 0); - double gstar(const double temp, const double pres, const int ifunc = 0); + 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: @@ -165,7 +177,7 @@ namespace Cantera { /*! * internal temporary variable */ - double m_densWaterSS; + mutable double m_densWaterSS; /** * Pointer to the water property calculator @@ -181,11 +193,21 @@ namespace Cantera { doublereal r_e_j; - //! Value of deltaG at Tr and Pr (cal gmol-1) + //! 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_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) /*! diff --git a/Cantera/src/thermo/PDSS.cpp b/Cantera/src/thermo/PDSS.cpp index e16b55e39..257f3dcff 100644 --- a/Cantera/src/thermo/PDSS.cpp +++ b/Cantera/src/thermo/PDSS.cpp @@ -234,8 +234,7 @@ namespace Cantera { void PDSS::initThermo() { } - void PDSS:: - setParametersFromXML(const XML_Node& eosdata) { + void PDSS::setParametersFromXML(const XML_Node& eosdata) { } // Return the molar enthalpy in units of J kmol-1