diff --git a/Cantera/src/thermo/HKFT_PDSS.cpp b/Cantera/src/thermo/HKFT_PDSS.cpp index 6e4f79482..e1c3f4a52 100644 --- a/Cantera/src/thermo/HKFT_PDSS.cpp +++ b/Cantera/src/thermo/HKFT_PDSS.cpp @@ -6,6 +6,7 @@ #include "ctml.h" #include "HKFT_PDSS.h" #include "WaterPDSS.h" +#include "WaterProps.h" #include "ThermoPhase.h" @@ -145,6 +146,21 @@ namespace Cantera { } 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:: @@ -345,7 +361,60 @@ namespace Cantera { throw CanteraError("HKFT_PDSS::satPressure()", "unimplemented"); return (0.0); } + + + double HKFT_PDSS::deltaG() { + + 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 Y_pr_tr = -5.799E-5; + //double Z_pr_tr = -0.0127803; + + 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 = m_deltaG_tr_pr + sterm + c1term + a1term + a2term + c2term + a3term + a4term + wterm + wrterm + yterm; + + // Convert to Joules / kg + 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 diff --git a/Cantera/src/thermo/HKFT_PDSS.h b/Cantera/src/thermo/HKFT_PDSS.h index 4a30d73db..0b304aeda 100644 --- a/Cantera/src/thermo/HKFT_PDSS.h +++ b/Cantera/src/thermo/HKFT_PDSS.h @@ -139,6 +139,9 @@ namespace Cantera { virtual void initThermo(); virtual void setParametersFromXML(const XML_Node& eosdata); + 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; @@ -155,6 +158,9 @@ namespace Cantera { */ WaterPDSS *m_waterSS; + //! Current value of the pressure for this object + doublereal m_pres; + //! density of standard-state water /*! * internal temporary variable @@ -167,6 +173,60 @@ namespace Cantera { 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 at Tr and Pr (cal gmol-1) + /*! + * Tr = 298.15 Pr = 1 atm + */ + doublereal m_deltaG_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; }; diff --git a/Cantera/src/thermo/WaterProps.h b/Cantera/src/thermo/WaterProps.h index 7bf9601a1..f486ba8d6 100644 --- a/Cantera/src/thermo/WaterProps.h +++ b/Cantera/src/thermo/WaterProps.h @@ -135,7 +135,7 @@ namespace Cantera { * value at 25C and 1 atm, relEps = 78.38 * */ - static double relEpsilon(double T, double P_pascal, int ifunc = 0); + double relEpsilon(double T, double P_pascal, int ifunc = 0); //! ADebye calculates the value of A_Debye as a function