From 4c630fc59211b5ebe878afef99baca7637c81a3f Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 9 Aug 2017 16:37:11 -0400 Subject: [PATCH] [Thermo] Make PDSS_HKFT constructible without XML Also fixes PDSS_HKFT to actually make use of the units specified for each input constant. --- include/cantera/thermo/PDSS_HKFT.h | 20 ++- src/thermo/PDSS_HKFT.cpp | 79 +++++----- test/thermo/phaseConstructors.cpp | 139 ++++++++++++++---- test/thermo/thermoParameterizations.cpp | 2 + .../simpleTransport/HMW_NaCl_pdss.xml | 24 +-- 5 files changed, 181 insertions(+), 83 deletions(-) diff --git a/include/cantera/thermo/PDSS_HKFT.h b/include/cantera/thermo/PDSS_HKFT.h index 7d03f3585..91decd014 100644 --- a/include/cantera/thermo/PDSS_HKFT.h +++ b/include/cantera/thermo/PDSS_HKFT.h @@ -13,11 +13,11 @@ #define CT_PDSS_HKFT_H #include "PDSS.h" +#include "WaterProps.h" namespace Cantera { class PDSS_Water; -class WaterProps; //! Class for pressure dependent standard states corresponding to //! ionic solutes in electrolyte water. @@ -87,6 +87,24 @@ public: virtual bool useSTITbyPDSS() const { return true; } virtual void initThermo(); + //! Set enthalpy of formation at Pr, Tr [J/kmol] + void setDeltaH0(double dh0); + + //! Set Gibbs free energy of formation at Pr, Tr [J/kmol] + void setDeltaG0(double dg0); + + //! Set entropy of formation at Pr, Tr [J/kmol/K] + void setS0(double s0); + + //! Set "a" coefficients (array of 4 elements). Units of each coefficient + //! are [J/kmol/Pa, J/kmol, J*K/kmol/Pa, J*K/kmol] + void set_a(double* a); + + //! Set "c" coefficients (array of 2 elements). Units of each coefficient + //! are [J/kmol/K, J*K/kmol] + void set_c(double* c); + void setOmega(double omega); //!< Set omega [J/kmol] + void setParametersFromXML(const XML_Node& speciesNode); //! This utility function reports back the type of parameterization and diff --git a/src/thermo/PDSS_HKFT.cpp b/src/thermo/PDSS_HKFT.cpp index 2b48b7337..c38a19d0e 100644 --- a/src/thermo/PDSS_HKFT.cpp +++ b/src/thermo/PDSS_HKFT.cpp @@ -316,6 +316,34 @@ void PDSS_HKFT::initThermo() } } +void PDSS_HKFT::setDeltaH0(double dh0) { + m_deltaH_formation_tr_pr = dh0 / toSI("cal/gmol"); +} + +void PDSS_HKFT::setDeltaG0(double dg0) { + m_deltaG_formation_tr_pr = dg0 / toSI("cal/gmol"); +} + +void PDSS_HKFT::setS0(double s0) { + m_Entrop_tr_pr = s0 / toSI("cal/gmol/K"); +} + +void PDSS_HKFT::set_a(double* a) { + m_a1 = a[0] / toSI("cal/gmol/bar"); + m_a2 = a[1] / toSI("cal/gmol"); + m_a3 = a[2] / toSI("cal-K/gmol/bar"); + m_a4 = a[3] / toSI("cal-K/gmol"); +} + +void PDSS_HKFT::set_c(double* c) { + m_c1 = c[0] / toSI("cal/gmol/K"); + m_c2 = c[1] / toSI("cal-K/gmol"); +} + +void PDSS_HKFT::setOmega(double omega) { + m_omega_pr_tr = omega / toSI("cal/gmol"); +} + void PDSS_HKFT::setParametersFromXML(const XML_Node& speciesNode) { PDSS::setParametersFromXML(speciesNode); @@ -357,20 +385,17 @@ void PDSS_HKFT::setParametersFromXML(const XML_Node& speciesNode) } if (hh->hasChild("DG0_f_Pr_Tr")) { - doublereal val = getFloat(*hh, "DG0_f_Pr_Tr"); - m_deltaG_formation_tr_pr = val; + setDeltaG0(getFloat(*hh, "DG0_f_Pr_Tr", "toSI")); hasDGO = 1; } if (hh->hasChild("DH0_f_Pr_Tr")) { - doublereal val = getFloat(*hh, "DH0_f_Pr_Tr"); - m_deltaH_formation_tr_pr = val; + setDeltaH0(getFloat(*hh, "DH0_f_Pr_Tr", "toSI")); hasDHO = 1; } if (hh->hasChild("S0_Pr_Tr")) { - doublereal val = getFloat(*hh, "S0_Pr_Tr"); - m_Entrop_tr_pr= val; + setS0(getFloat(*hh, "S0_Pr_Tr", "toSI")); hasSO = 1; } @@ -384,42 +409,14 @@ void PDSS_HKFT::setParametersFromXML(const XML_Node& speciesNode) "standardState model for species isn't hkft: " + speciesNode.name()); } - if (ss->hasChild("a1")) { - m_a1 = getFloat(*ss, "a1"); - } else { - throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a1 field"); - } - if (ss->hasChild("a2")) { - m_a2 = getFloat(*ss, "a2"); - } else { - throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a2 field"); - } - if (ss->hasChild("a3")) { - m_a3 = getFloat(*ss, "a3"); - } else { - throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a3 field"); - } - if (ss->hasChild("a4")) { - m_a4 = getFloat(*ss, "a4"); - } else { - throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a4 field"); - } + double a[4] = {getFloat(*ss, "a1", "toSI"), getFloat(*ss, "a2", "toSI"), + getFloat(*ss, "a3", "toSI"), getFloat(*ss, "a4", "toSI")}; + set_a(a); - if (ss->hasChild("c1")) { - m_c1 = getFloat(*ss, "c1"); - } else { - throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing c1 field"); - } - if (ss->hasChild("c2")) { - m_c2 = getFloat(*ss, "c2"); - } else { - throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing c2 field"); - } - if (ss->hasChild("omega_Pr_Tr")) { - m_omega_pr_tr = getFloat(*ss, "omega_Pr_Tr"); - } else { - throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing omega_Pr_Tr field"); - } + double c[2] = {getFloat(*ss, "c1", "toSI"), getFloat(*ss, "c2", "toSI")}; + set_c(c); + + setOmega(getFloat(*ss, "omega_Pr_Tr", "toSI")); int isum = hasDGO + hasDHO + hasSO; if (isum < 2) { diff --git a/test/thermo/phaseConstructors.cpp b/test/thermo/phaseConstructors.cpp index c094b75c3..1849ba73f 100644 --- a/test/thermo/phaseConstructors.cpp +++ b/test/thermo/phaseConstructors.cpp @@ -18,6 +18,7 @@ #include "cantera/thermo/LatticeSolidPhase.h" #include "cantera/thermo/IdealSolidSolnPhase.h" #include "cantera/thermo/HMWSoln.h" +#include "cantera/thermo/PDSS_HKFT.h" #include "cantera/thermo/NasaPoly2.h" #include "cantera/thermo/ConstCpPoly.h" @@ -551,6 +552,38 @@ TEST(IdealSolidSolnPhase, fromScratch) EXPECT_NEAR(p.gibbs_mole(), -313642293.1654253, 1e-4); } +static void set_hmw_interactions(HMWSoln& p) { + double beta0_nacl[] = {0.0765, 0.008946, -3.3158E-6, -777.03, -4.4706}; + double beta1_nacl[] = {0.2664, 6.1608E-5, 1.0715E-6, 0.0, 0.0}; + double beta2_nacl[] = {0.0, 0.0, 0.0, 0.0, 0.0}; + double cphi_nacl[] = {0.00127, -4.655E-5, 0.0, 33.317, 0.09421}; + p.setBinarySalt("Na+", "Cl-", 5, beta0_nacl, beta1_nacl, beta2_nacl, + cphi_nacl, 2.0, 0.0); + + double beta0_hcl[] = {0.1775, 0.0, 0.0, 0.0, 0.0}; + double beta1_hcl[] = {0.2945, 0.0, 0.0, 0.0, 0.0}; + double beta2_hcl[] = {0.0, 0.0, 0.0, 0.0, 0.0}; + double cphi_hcl[] = {0.0008, 0.0, 0.0, 0.0, 0.0}; + p.setBinarySalt("H+", "Cl-", 5, beta0_hcl, beta1_hcl, beta2_hcl, + cphi_hcl, 2.0, 0.0); + + double beta0_naoh[] = {0.0864, 0.0, 0.0, 0.0, 0.0}; + double beta1_naoh[] = {0.253, 0.0, 0.0, 0.0, 0.0}; + double beta2_naoh[] = {0.0, 0.0, 0.0, 0.0, 0.0}; + double cphi_naoh[] = {0.0044, 0.0, 0.0, 0.0, 0.0}; + p.setBinarySalt("Na+", "OH-", 5, beta0_naoh, beta1_naoh, beta2_naoh, + cphi_naoh, 2.0, 0.0); + + double theta_cloh[] = {-0.05, 0.0, 0.0, 0.0, 0.0}; + double psi_nacloh[] = {-0.006, 0.0, 0.0, 0.0, 0.0}; + double theta_nah[] = {0.036, 0.0, 0.0, 0.0, 0.0}; + double psi_clnah[] = {-0.004, 0.0, 0.0, 0.0, 0.0}; + p.setTheta("Cl-", "OH-", 5, theta_cloh); + p.setPsi("Na+", "Cl-", "OH-", 5, psi_nacloh); + p.setTheta("Na+", "H+", 5, theta_nah); + p.setPsi("Cl-", "Na+", "H+", 5, psi_clnah); +} + TEST(HMWSoln, fromScratch) { // Regression test based on HMW_test_3 @@ -583,35 +616,7 @@ TEST(HMWSoln, fromScratch) p.setA_Debye(1.175930); p.initThermo(); - double beta0_nacl[] = {0.0765, 0.008946, -3.3158E-6, -777.03, -4.4706}; - double beta1_nacl[] = {0.2664, 6.1608E-5, 1.0715E-6, 0.0, 0.0}; - double beta2_nacl[] = {0.0, 0.0, 0.0, 0.0, 0.0}; - double cphi_nacl[] = {0.00127, -4.655E-5, 0.0, 33.317, 0.09421}; - p.setBinarySalt("Na+", "Cl-", 5, beta0_nacl, beta1_nacl, beta2_nacl, - cphi_nacl, 2.0, 0.0); - - double beta0_hcl[] = {0.1775, 0.0, 0.0, 0.0, 0.0}; - double beta1_hcl[] = {0.2945, 0.0, 0.0, 0.0, 0.0}; - double beta2_hcl[] = {0.0, 0.0, 0.0, 0.0, 0.0}; - double cphi_hcl[] = {0.0008, 0.0, 0.0, 0.0, 0.0}; - p.setBinarySalt("H+", "Cl-", 5, beta0_hcl, beta1_hcl, beta2_hcl, - cphi_hcl, 2.0, 0.0); - - double beta0_naoh[] = {0.0864, 0.0, 0.0, 0.0, 0.0}; - double beta1_naoh[] = {0.253, 0.0, 0.0, 0.0, 0.0}; - double beta2_naoh[] = {0.0, 0.0, 0.0, 0.0, 0.0}; - double cphi_naoh[] = {0.0044, 0.0, 0.0, 0.0, 0.0}; - p.setBinarySalt("Na+", "OH-", 5, beta0_naoh, beta1_naoh, beta2_naoh, - cphi_naoh, 2.0, 0.0); - - double theta_cloh[] = {-0.05, 0.0, 0.0, 0.0, 0.0}; - double psi_nacloh[] = {-0.006, 0.0, 0.0, 0.0, 0.0}; - double theta_nah[] = {0.036, 0.0, 0.0, 0.0, 0.0}; - double psi_clnah[] = {-0.004, 0.0, 0.0, 0.0, 0.0}; - p.setTheta("Cl-", "OH-", 5, theta_cloh); - p.setPsi("Na+", "Cl-", "OH-", 5, psi_nacloh); - p.setTheta("Na+", "H+", 5, theta_nah); - p.setPsi("Cl-", "Na+", "H+", 5, psi_clnah); + set_hmw_interactions(p); p.setMolalitiesByName("Na+:6.0997 Cl-:6.0996986044628 H+:2.1628E-9 OH-:1.3977E-6"); p.setState_TP(150 + 273.15, 101325); @@ -638,4 +643,80 @@ TEST(HMWSoln, fromScratch) } } +TEST(HMWSoln, fromScratch_HKFT) +{ + HMWSoln p; + p.addUndefinedElements(); + auto sH2O = make_species("H2O(l)", "H:2, O:1", h2oliq_nasa_coeffs); + auto sNa = make_species("Na+", "Na:1, E:-1", 0.0, + 298.15, -125.5213, 333.15, -125.5213, 1e5); + sNa->charge = 1; + auto sCl = make_species("Cl-", "Cl:1, E:1", 0.0, + 298.15, -52.8716, 333.15, -52.8716, 1e5); + sCl->charge = -1; + auto sH = make_species("H+", "H:1, E:-1", 0.0, 298.15, 0.0, 333.15, 0.0, 1e5); + sH->charge = 1; + auto sOH = make_species("OH-", "O:1, H:1, E:1", 0.0, + 298.15, -91.523, 333.15, -91.523, 1e5); + sOH->charge = -1; + for (auto& s : {sH2O, sNa, sCl, sH, sOH}) { + p.addSpecies(s); + } + double h0[] = {-57433, Undef, 0.0, -54977}; + double g0[] = {Undef, -31379, 0.0, -37595}; + double s0[] = {13.96, 13.56, Undef, -2.56}; + double a[][4] = {{0.1839, -228.5, 3.256, -27260}, + {0.4032, 480.1, 5.563, -28470}, + {0.0, 0.0, 0.0, 0.0}, + {0.12527, 7.38, 1.8423, -27821}}; + double c[][2] = {{18.18, -29810}, {-4.4, -57140}, {0.0, 0.0}, {4.15, -103460}}; + double omega[] = {33060, 145600, 0.0, 172460}; + + std::unique_ptr ss(new PDSS_Water()); + p.installPDSS(0, std::move(ss)); + for (size_t k = 0; k < 4; k++) { + std::unique_ptr ss(new PDSS_HKFT()); + if (h0[k] != Undef) { + ss->setDeltaH0(h0[k] * toSI("cal/gmol")); + } + if (g0[k] != Undef) { + ss->setDeltaG0(g0[k] * toSI("cal/gmol")); + } + if (s0[k] != Undef) { + ss->setS0(s0[k] * toSI("cal/gmol/K")); + } + a[k][0] *= toSI("cal/gmol/bar"); + a[k][1] *= toSI("cal/gmol"); + a[k][2] *= toSI("cal-K/gmol/bar"); + a[k][3] *= toSI("cal-K/gmol"); + c[k][0] *= toSI("cal/gmol/K"); + c[k][1] *= toSI("cal-K/gmol"); + ss->set_a(a[k]); + ss->set_c(c[k]); + ss->setOmega(omega[k] * toSI("cal/gmol")); + p.installPDSS(k+1, std::move(ss)); + } + p.setPitzerTempModel("complex"); + p.setA_Debye(-1); + p.initThermo(); + + set_hmw_interactions(p); + p.setMolalitiesByName("Na+:6.0954 Cl-:6.0954 H+:2.1628E-9 OH-:1.3977E-6"); + p.setState_TP(50 + 273.15, 101325); + + size_t N = p.nSpecies(); + vector_fp mv(N), h(N), mu(N), ac(N), acoeff(N); + p.getPartialMolarVolumes(mv.data()); + p.getPartialMolarEnthalpies(h.data()); + p.getChemPotentials(mu.data()); + p.getActivities(ac.data()); + p.getActivityCoefficients(acoeff.data()); + + double mvRef[] = {0.01815224, 0.00157182, 0.01954605, 0.00173137, -0.0020266}; + + for (size_t k = 0; k < N; k++) { + EXPECT_NEAR(mv[k], mvRef[k], 2e-8); + } +} + } // namespace Cantera diff --git a/test/thermo/thermoParameterizations.cpp b/test/thermo/thermoParameterizations.cpp index 48649a7c7..93a2254db 100644 --- a/test/thermo/thermoParameterizations.cpp +++ b/test/thermo/thermoParameterizations.cpp @@ -4,8 +4,10 @@ #include "cantera/thermo/ConstCpPoly.h" #include "cantera/thermo/NasaPoly2.h" #include "cantera/thermo/ShomatePoly.h" +#include "cantera/thermo/PDSS_HKFT.h" #include "cantera/base/stringUtils.h" #include "thermo_data.h" +#include using namespace Cantera; diff --git a/test_problems/simpleTransport/HMW_NaCl_pdss.xml b/test_problems/simpleTransport/HMW_NaCl_pdss.xml index 3f08a1ff3..02b4de49e 100644 --- a/test_problems/simpleTransport/HMW_NaCl_pdss.xml +++ b/test_problems/simpleTransport/HMW_NaCl_pdss.xml @@ -141,10 +141,10 @@ 0.1839 -228.5 - 3.256 - -27260. + 3.256 + -27260. 18.18 - -29810. + -29810. 33060. @@ -168,10 +168,10 @@ 0.4032 480.1 - 5.563 - -28470. + 5.563 + -28470. -4.4 - -57140. + -57140. 145600. @@ -195,10 +195,10 @@ 0.0 0.0 - 0.0 - 0.0 + 0.0 + 0.0 0.0 - 0.0 + 0.0 0.0 @@ -223,10 +223,10 @@ 0.12527 7.38 - 1.8423 - -27821 + 1.8423 + -27821 4.15 - -103460. + -103460. 172460.