diff --git a/include/cantera/thermo/VPSSMgr.h b/include/cantera/thermo/VPSSMgr.h index da785f04b..7155a4d2b 100644 --- a/include/cantera/thermo/VPSSMgr.h +++ b/include/cantera/thermo/VPSSMgr.h @@ -677,8 +677,8 @@ protected: //! properties were calculated at. mutable doublereal m_plast; - //! Reference pressure (Pa) must be the same for all species - defaults to 1 atm. - mutable doublereal m_p0; + //! Reference pressure (Pa) for each species + mutable vector_fp m_p0; //! minimum temperature for the standard state calculations doublereal m_minTemp; diff --git a/include/cantera/thermo/VPStandardStateTP.h b/include/cantera/thermo/VPStandardStateTP.h index a6d0a4133..c70c1b4bc 100644 --- a/include/cantera/thermo/VPStandardStateTP.h +++ b/include/cantera/thermo/VPStandardStateTP.h @@ -287,12 +287,6 @@ protected: //! were calculated at. mutable doublereal m_Plast_ss; - /*! - * Reference pressure (Pa) must be the same for all species - * - defaults to OneAtm - */ - doublereal m_P0; - // -> suggest making this private! //! Pointer to the VPSS manager that calculates all of the standard state //! info efficiently. diff --git a/src/thermo/VPSSMgr.cpp b/src/thermo/VPSSMgr.cpp index 6a6f732aa..98e871f19 100644 --- a/src/thermo/VPSSMgr.cpp +++ b/src/thermo/VPSSMgr.cpp @@ -28,7 +28,6 @@ VPSSMgr::VPSSMgr(VPStandardStateTP* vptp_ptr, MultiSpeciesThermo* spthermo) : m_spthermo(spthermo), m_tlast(-1.0), m_plast(-1.0), - m_p0(-1.0), m_minTemp(-1.0), m_maxTemp(1.0E8), m_useTmpRefStateStorage(false), @@ -266,10 +265,10 @@ void VPSSMgr::initLengths() void VPSSMgr::initThermoXML(XML_Node& phaseNode, const std::string& id) { - const PDSS* kPDSS = m_vptp_ptr->providePDSS(0); - m_p0 = kPDSS->refPressure(); for (size_t i = 0; i < m_kk; i++) { const PDSS* kPDSS = m_vptp_ptr->providePDSS(i); + m_p0.resize(std::max(m_p0.size(), i+1)); + m_p0[i] = kPDSS->refPressure(); m_minTemp = std::max(m_minTemp, kPDSS->minTemp()); m_maxTemp = std::min(m_maxTemp, kPDSS->maxTemp()); } @@ -281,9 +280,8 @@ void VPSSMgr::installSTSpecies(size_t k, const XML_Node& s, shared_ptr stit(newSpeciesThermoInterpType(s.child("thermo"))); stit->validate(s["name"]); m_spthermo->install_STIT(k, stit); - if (m_p0 < 0.0) { - m_p0 = m_spthermo->refPressure(k); - } + m_p0.resize(std::max(m_p0.size(), k+1)); + m_p0[k] = m_spthermo->refPressure(k); } PDSS* VPSSMgr::createInstallPDSS(size_t k, const XML_Node& s, @@ -315,7 +313,7 @@ doublereal VPSSMgr::refPressure(size_t k) const if (k != npos) { return m_vptp_ptr->providePDSS(k)->refPressure(); } - return m_p0; + return m_p0[0]; } } diff --git a/src/thermo/VPSSMgr_ConstVol.cpp b/src/thermo/VPSSMgr_ConstVol.cpp index 69c9b136d..c6edd0895 100644 --- a/src/thermo/VPSSMgr_ConstVol.cpp +++ b/src/thermo/VPSSMgr_ConstVol.cpp @@ -34,10 +34,9 @@ VPSSMgr_ConstVol::VPSSMgr_ConstVol(VPStandardStateTP* vp_ptr, MultiSpeciesThermo */ void VPSSMgr_ConstVol::_updateStandardStateThermo() { - doublereal del_pRT = (m_plast - m_p0) / (GasConstant * m_tlast); - for (size_t k = 0; k < m_kk; k++) { - m_hss_RT[k] = m_h0_RT[k] + del_pRT * m_Vss[k]; + m_hss_RT[k] = m_h0_RT[k] + + (m_plast - m_p0[k]) / (GasConstant * m_tlast) * 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]; diff --git a/src/thermo/VPSSMgr_General.cpp b/src/thermo/VPSSMgr_General.cpp index 6edf0e4bc..3605d5e52 100644 --- a/src/thermo/VPSSMgr_General.cpp +++ b/src/thermo/VPSSMgr_General.cpp @@ -146,10 +146,8 @@ PDSS* VPSSMgr_General::createInstallPDSS(size_t k, const XML_Node& speciesNode, m_kk = std::max(m_kk, k+1); m_minTemp = std::max(m_minTemp, kPDSS->minTemp()); m_maxTemp = std::min(m_maxTemp, kPDSS->maxTemp()); - doublereal p0 = kPDSS->refPressure(); - if (k == 0) { - m_p0 = p0; - } + m_p0.resize(std::max(m_p0.size(), k+1)); + m_p0[k] = kPDSS->refPressure(); return kPDSS; } diff --git a/src/thermo/VPSSMgr_IdealGas.cpp b/src/thermo/VPSSMgr_IdealGas.cpp index 40d184256..7423f2d3c 100644 --- a/src/thermo/VPSSMgr_IdealGas.cpp +++ b/src/thermo/VPSSMgr_IdealGas.cpp @@ -43,13 +43,12 @@ void VPSSMgr_IdealGas::getStandardVolumes(doublereal* vol) const void VPSSMgr_IdealGas::_updateStandardStateThermo() { - doublereal pp = log(m_plast / m_p0); doublereal v = temperature() *GasConstant /m_plast; for (size_t k = 0; k < m_kk; k++) { m_hss_RT[k] = m_h0_RT[k]; m_cpss_R[k] = m_cp0_R[k]; - m_sss_R[k] = m_s0_R[k] - pp; + m_sss_R[k] = m_s0_R[k] - log(m_plast / m_p0[k]); m_gss_RT[k] = m_hss_RT[k] - m_sss_R[k]; m_Vss[k] = v; } @@ -76,7 +75,7 @@ PDSS* VPSSMgr_IdealGas::createInstallPDSS(size_t k, const XML_Node& speciesNode, PDSS* kPDSS = new PDSS_IdealGas(m_vptp_ptr, k, speciesNode, *phaseNode_ptr, true); - m_p0 = m_spthermo->refPressure(k); + m_p0[k] = m_spthermo->refPressure(k); return kPDSS; } diff --git a/src/thermo/VPSSMgr_Water_ConstVol.cpp b/src/thermo/VPSSMgr_Water_ConstVol.cpp index a68dcd725..ed934a33f 100644 --- a/src/thermo/VPSSMgr_Water_ConstVol.cpp +++ b/src/thermo/VPSSMgr_Water_ConstVol.cpp @@ -33,9 +33,9 @@ VPSSMgr_Water_ConstVol::VPSSMgr_Water_ConstVol(VPStandardStateTP* vp_ptr, void VPSSMgr_Water_ConstVol::getEnthalpy_RT_ref(doublereal* hrt) const { // Everything should be OK except for the water SS - m_p0 = m_waterSS->pref_safe(m_tlast); - if (m_p0 != m_plast) { - m_waterSS->setState_TP(m_tlast, m_p0); + m_p0[0] = m_waterSS->pref_safe(m_tlast); + if (m_p0[0] != m_plast) { + m_waterSS->setState_TP(m_tlast, m_p0[0]); m_h0_RT[0] = (m_waterSS->enthalpy_mole()) / (GasConstant * m_tlast); m_waterSS->setState_TP(m_tlast, m_plast); } else { @@ -47,9 +47,9 @@ void VPSSMgr_Water_ConstVol::getEnthalpy_RT_ref(doublereal* hrt) const void VPSSMgr_Water_ConstVol::getGibbs_RT_ref(doublereal* grt) const { // Everything should be OK except for the water SS - m_p0 = m_waterSS->pref_safe(m_tlast); - if (m_p0 != m_plast) { - m_waterSS->setState_TP(m_tlast, m_p0); + m_p0[0] = m_waterSS->pref_safe(m_tlast); + if (m_p0[0] != m_plast) { + m_waterSS->setState_TP(m_tlast, m_p0[0]); m_g0_RT[0] = (m_waterSS->gibbs_mole()) / (GasConstant * m_tlast); m_waterSS->setState_TP(m_tlast, m_plast); } else { @@ -69,9 +69,9 @@ void VPSSMgr_Water_ConstVol::getGibbs_ref(doublereal* g) const void VPSSMgr_Water_ConstVol::getEntropy_R_ref(doublereal* sr) const { // Everything should be OK except for the water SS - m_p0 = m_waterSS->pref_safe(m_tlast); - if (m_p0 != m_plast) { - m_waterSS->setState_TP(m_tlast, m_p0); + m_p0[0] = m_waterSS->pref_safe(m_tlast); + if (m_p0[0] != m_plast) { + m_waterSS->setState_TP(m_tlast, m_p0[0]); m_s0_R[0] = (m_waterSS->entropy_mole()) / GasConstant; m_waterSS->setState_TP(m_tlast, m_plast); } else { @@ -83,9 +83,9 @@ void VPSSMgr_Water_ConstVol::getEntropy_R_ref(doublereal* sr) const void VPSSMgr_Water_ConstVol::getCp_R_ref(doublereal* cpr) const { // Everything should be OK except for the water SS - m_p0 = m_waterSS->pref_safe(m_tlast); - if (m_p0 != m_plast) { - m_waterSS->setState_TP(m_tlast, m_p0); + m_p0[0] = m_waterSS->pref_safe(m_tlast); + if (m_p0[0] != m_plast) { + m_waterSS->setState_TP(m_tlast, m_p0[0]); m_cp0_R[0] = (m_waterSS->cp_mole()) / GasConstant; m_waterSS->setState_TP(m_tlast, m_plast); } else { @@ -97,9 +97,9 @@ void VPSSMgr_Water_ConstVol::getCp_R_ref(doublereal* cpr) const void VPSSMgr_Water_ConstVol::getStandardVolumes_ref(doublereal* vol) const { // Everything should be OK except for the water SS - m_p0 = m_waterSS->pref_safe(m_tlast); - if (m_p0 != m_plast) { - m_waterSS->setState_TP(m_tlast, m_p0); + m_p0[0] = m_waterSS->pref_safe(m_tlast); + if (m_p0[0] != m_plast) { + m_waterSS->setState_TP(m_tlast, m_p0[0]); m_V0[0] = m_vptp_ptr->molecularWeight(0) / m_waterSS->density(); m_waterSS->setState_TP(m_tlast, m_plast); } else { @@ -110,27 +110,25 @@ void VPSSMgr_Water_ConstVol::getStandardVolumes_ref(doublereal* vol) const void VPSSMgr_Water_ConstVol::_updateRefStateThermo() const { - m_p0 = m_waterSS->pref_safe(m_tlast); + m_p0[0] = m_waterSS->pref_safe(m_tlast); m_spthermo->update(m_tlast, &m_cp0_R[0], &m_h0_RT[0], &m_s0_R[0]); for (size_t k = 0; k < m_kk; k++) { m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; m_vptp_ptr->providePDSS(k)->setTemperature(m_tlast); } - m_waterSS->setState_TP(m_tlast, m_p0); + m_waterSS->setState_TP(m_tlast, m_p0[0]); m_h0_RT[0] = (m_waterSS->enthalpy_mole()) / (GasConstant * m_tlast); m_s0_R[0] = (m_waterSS->entropy_mole()) / GasConstant; m_cp0_R[0] = (m_waterSS->cp_mole()) / GasConstant; - m_g0_RT[0] = (m_hss_RT[0] - m_sss_R[0]); - m_V0[0] = m_vptp_ptr->molecularWeight(0) / (m_waterSS->density()); + m_g0_RT[0] = (m_h0_RT[0] - m_s0_R[0]); + m_V0[0] = m_waterSS->molarVolume(); m_waterSS->setState_TP(m_tlast, m_plast); } void VPSSMgr_Water_ConstVol::_updateStandardStateThermo() { - doublereal del_pRT = (m_plast - OneAtm) / (GasConstant * m_tlast); - for (size_t k = 1; k < m_kk; k++) { - m_hss_RT[k] = m_h0_RT[k] + del_pRT * m_Vss[k]; + m_hss_RT[k] = m_h0_RT[k] + (m_plast - m_p0[k]) / (GasConstant * m_tlast) * 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]; @@ -144,7 +142,7 @@ void VPSSMgr_Water_ConstVol::_updateStandardStateThermo() m_sss_R[0] = (m_waterSS->entropy_mole()) / GasConstant; 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_vptp_ptr->molecularWeight(0) / m_waterSS->density()); + m_Vss[0] = m_waterSS->molarVolume(); } void VPSSMgr_Water_ConstVol::initThermoXML(XML_Node& phaseNode, @@ -225,6 +223,7 @@ PDSS* VPSSMgr_Water_ConstVol::createInstallPDSS(size_t k, // instantiate a new kPDSS object kPDSS = new PDSS_ConstVol(m_vptp_ptr, k, speciesNode, *phaseNode_ptr, true); + m_p0[k] = kPDSS->refPressure(); } return kPDSS; } diff --git a/src/thermo/VPSSMgr_Water_HKFT.cpp b/src/thermo/VPSSMgr_Water_HKFT.cpp index c3a4fe066..a26dab75c 100644 --- a/src/thermo/VPSSMgr_Water_HKFT.cpp +++ b/src/thermo/VPSSMgr_Water_HKFT.cpp @@ -110,8 +110,8 @@ void VPSSMgr_Water_HKFT::updateRefStateThermo() const void VPSSMgr_Water_HKFT::_updateRefStateThermo() const { - m_p0 = m_waterSS->pref_safe(m_tlast); - m_waterSS->setState_TP(m_tlast, m_p0); + m_p0[0] = m_waterSS->pref_safe(m_tlast); + m_waterSS->setState_TP(m_tlast, m_p0[0]); m_h0_RT[0] = (m_waterSS->enthalpy_mole()) / (GasConstant * m_tlast); m_s0_R[0] = (m_waterSS->entropy_mole()) / GasConstant; m_cp0_R[0] = (m_waterSS->cp_mole()) / GasConstant; @@ -120,7 +120,7 @@ void VPSSMgr_Water_HKFT::_updateRefStateThermo() const PDSS* ps; for (size_t k = 1; k < m_kk; k++) { ps = m_vptp_ptr->providePDSS(k); - ps->setState_TP(m_tlast, m_p0); + ps->setState_TP(m_tlast, m_p0[0]); m_cp0_R[k] = ps->cp_R(); m_s0_R[k] = ps->entropy_mole() / GasConstant; m_g0_RT[k] = ps->gibbs_RT(); diff --git a/src/thermo/VPStandardStateTP.cpp b/src/thermo/VPStandardStateTP.cpp index a25d58fca..6373e097a 100644 --- a/src/thermo/VPStandardStateTP.cpp +++ b/src/thermo/VPStandardStateTP.cpp @@ -20,8 +20,7 @@ namespace Cantera VPStandardStateTP::VPStandardStateTP() : m_Pcurrent(OneAtm), m_Tlast_ss(-1.0), - m_Plast_ss(-1.0), - m_P0(OneAtm) + m_Plast_ss(-1.0) { } diff --git a/test/data/HMW_NaCl.xml b/test/data/HMW_NaCl.xml new file mode 100644 index 000000000..3d9e100f6 --- /dev/null +++ b/test/data/HMW_NaCl.xml @@ -0,0 +1,252 @@ + + + + + H2O(L) Cl- H+ Na+ OH- + + + 298.15 + 101325.0 + + Na+:6.0954 + Cl-:6.0954 + H+:2.1628E-9 + OH-:1.3977E-6 + + + + + + + 1.175930 + 3.28640E9 + + + + 0.0765 + 0.2664 + 0.0 + 0.00127 + 2.0 + + + + 0.1775 + 0.2945 + 0.0 + 0.0008 + 2.0 + + + + 0.0864 + 0.253 + 0.0 + 0.0044 + 2.0 + + + + -0.05 + + + + -0.05 + -0.006 + + + + 0.036 + + + + 0.036 + -0.004 + + + + H2O(L) + + O H C E Fe Si N Na Cl + + + + + H2O(L) Cl- H+ Na+ OH- + + + 298.15 + 101325.0 + + Na+:6.0954 + Cl-:6.0954 + H+:2.1628E-9 + OH-:1.3977E-6 + + + + + + + + + 1.175930 + + 3.28640E9 + + + + 0.0765 + 0.2664 + 0.0 + 0.00127 + 2.0 + + + + 0.1775 + 0.2945 + 0.0 + 0.0008 + 2.0 + + + + 0.0864 + 0.253 + 0.0 + 0.0044 + 2.0 + + + + -0.05 + + + + -0.05 + -0.006 + + + + 0.036 + + + + 0.036 + -0.004 + + + + H2O(L) + + O H C E Fe Si N Na Cl + + + + + + + H:2 O:1 + + + + 7.255750050E+01, -6.624454020E-01, 2.561987460E-03, -4.365919230E-06, + 2.781789810E-09, -4.188654990E+04, -2.882801370E+02 + + + + + + + + + Na:1 E:-1 + +1 + + + 0.0 + 2 + + -125.5213, -125.5213 + + + 298.15, 333.15 + + + + + 1.3 + + + + + Cl:1 E:1 + -1 + + 1.3 + + + + 0.0 + 2 + + -52.8716 , -52.8716 + + + 298.15, 333.15 + + + + + + + H:1 E:-1 + +1 + + 1.3 + + + + 0.0 + 2 + + 0.0 , 0.0 + + + 298.15, 333.15 + + + + + + + O:1 H:1 E:1 + -1 + + 1.3 + + + + 0.0 + 2 + + -91.523 , -91.523 + + + 298.15, 333.15 + + + + + + + + diff --git a/test/thermo/standardStateManagers.cpp b/test/thermo/standardStateManagers.cpp new file mode 100644 index 000000000..754539f8c --- /dev/null +++ b/test/thermo/standardStateManagers.cpp @@ -0,0 +1,59 @@ +#include "gtest/gtest.h" +#include "cantera/thermo/HMWSoln.h" + +using namespace Cantera; + +TEST(HMW, VPSSMgrGeneral_vs_VPSSMgrWater_ConstVol) +{ + // Calculations should give the same result using either the generic + // VPSSMgr_General class or one of the more specialized classes such as + // VPSSMgr_Water_ConstVol. + HMWSoln p1("../data/HMW_NaCl.xml", "water_constvol"); + HMWSoln p2("../data/HMW_NaCl.xml", "general"); + size_t n = p1.nSpecies(); + vector_fp molalities(n); + p1.getMolalities(molalities.data()); + molalities[2] = 2.1628E-9; + molalities[3] = 6.0997; + molalities[4] = 1.3977E-6; + molalities[1] = molalities[2] + molalities[3] - molalities[4]; + p1.setMolalities(molalities.data()); + p2.setMolalities(molalities.data()); + p1.setState_TP(310.15, 201325); + p2.setState_TP(310.15, 201325); + + vector_fp v1(n); + vector_fp v2(n); + p1.getStandardVolumes(v1.data()); + p2.getStandardVolumes(v2.data()); + for (size_t i = 0; i < n; i++) { + EXPECT_NEAR(v1[i], v2[i], 1e-9) << p1.speciesName(i); + } + + p1.getCp_R(v1.data()); + p2.getCp_R(v2.data()); + for (size_t i = 0; i < n; i++) { + EXPECT_NEAR(v1[i], v2[i], 1e-10) << p1.speciesName(i); + } + + p1.getEntropy_R(v1.data()); + p2.getEntropy_R(v2.data()); + for (size_t i = 0; i < n; i++) { + EXPECT_NEAR(v1[i], v2[i], 1e-10) << p1.speciesName(i); + } + + p1.getEnthalpy_RT(v1.data()); + p2.getEnthalpy_RT(v2.data()); + for (size_t i = 0; i < n; i++) { + EXPECT_NEAR(v1[i], v2[i], 1e-10) << p1.speciesName(i); + } + + p1.getChemPotentials_RT(v1.data()); + p2.getChemPotentials_RT(v2.data()); + for (size_t i = 0; i < n; i++) { + EXPECT_NEAR(v1[i], v2[i], 1e-10) << p1.speciesName(i); + } + + EXPECT_NEAR(p1.entropy_mole(), p2.entropy_mole(), 1e-7); + EXPECT_NEAR(p1.enthalpy_mole(), p2.enthalpy_mole(), 1e-4); +} diff --git a/test_problems/cathermo/HMW_test_1/output_noD_blessed.txt b/test_problems/cathermo/HMW_test_1/output_noD_blessed.txt index 0d2dff721..879760ab0 100644 --- a/test_problems/cathermo/HMW_test_1/output_noD_blessed.txt +++ b/test_problems/cathermo/HMW_test_1/output_noD_blessed.txt @@ -37,22 +37,22 @@ a2 = 3.04284e-10 Species Standard chemical potentials (kJ/gmol) ------------------------------------------------------------ H2O(L) -306.685728 - Cl- -131.066416 - H+ 0 - Na+ -311.16189 - OH- -226.88157 + Cl- -131.064694 + H+ 0.0017225 + Na+ -311.160167 + OH- -226.879848 ------------------------------------------------------------ Some DeltaSS values: Delta(mu_0) - NaCl(S): Na+ + Cl- -> NaCl(S): 9.597906 kJ/gmol - : 3.871747 (dimensionless) - : 1.681478 (dimensionless/ln10) + NaCl(S): Na+ + Cl- -> NaCl(S): 9.594461 kJ/gmol + : 3.870358 (dimensionless) + : 1.680875 (dimensionless/ln10) G0(NaCl(S)) = -432.6304 (fixed) - G0(Na+) = -311.1619 - G0(Cl-) = -131.0664 - OH-: H2O(L) - H+ -> OH-: 79.80416 kJ/gmol - : 32.1926 (dimensionless) - : 13.98107 (dimensionless/ln10) - G0(OH-) = -226.8816 - G0(H+) = 0 + G0(Na+) = -311.1602 + G0(Cl-) = -131.0647 + OH-: H2O(L) - H+ -> OH-: 79.8076 kJ/gmol + : 32.19399 (dimensionless) + : 13.98167 (dimensionless/ln10) + G0(OH-) = -226.8798 + G0(H+) = 0.0017225 G0(H2O(L)) = -306.6857 ------------------------------------------------------------ diff --git a/test_problems/cathermo/HMW_test_3/output_noD_blessed.txt b/test_problems/cathermo/HMW_test_3/output_noD_blessed.txt index 4e19d7b6a..1c8217aa3 100644 --- a/test_problems/cathermo/HMW_test_3/output_noD_blessed.txt +++ b/test_problems/cathermo/HMW_test_3/output_noD_blessed.txt @@ -38,16 +38,16 @@ Index Name MoleF MolalityCropped Charge Species Standard chemical potentials (kJ/gmol) ------------------------------------------------------------ H2O(L) -317.175788 - Cl- -186.016281 - H+ 0 - Na+ -441.617151 - OH- -322.002134 + Cl- -186.014558 + H+ 0.0017225 + Na+ -441.615429 + OH- -322.000412 ------------------------------------------------------------ Some DeltaSS values: Delta(mu_0) NaCl(S): Na+ + Cl- -> NaCl(S): 195 kJ/gmol - : 78.663 (dimensionless) - : 34.163 (dimensionless/ln10) - OH-: H2O(L) - H+ -> OH-: -4.8263 kJ/gmol - : -1.9469 (dimensionless) - : -0.84554 (dimensionless/ln10) + : 78.662 (dimensionless) + : 34.162 (dimensionless/ln10) + OH-: H2O(L) - H+ -> OH-: -4.8229 kJ/gmol + : -1.9455 (dimensionless) + : -0.84493 (dimensionless/ln10) ------------------------------------------------------------