diff --git a/Cantera/src/thermo/PDSS.cpp b/Cantera/src/thermo/PDSS.cpp index 1d464d586..8409882fe 100644 --- a/Cantera/src/thermo/PDSS.cpp +++ b/Cantera/src/thermo/PDSS.cpp @@ -195,13 +195,16 @@ namespace Cantera { void PDSS::initThermoXML(const XML_Node& phaseNode, std::string& id) { AssertThrow(m_tp != 0, "PDSS::initThermoXML()"); - m_vpssmgr_ptr = m_tp->provideVPSSMgr(); + m_p0 = m_vpssmgr_ptr->refPressure(m_spindex); + m_minTemp = m_vpssmgr_ptr->minTemp(m_spindex); + m_maxTemp = m_vpssmgr_ptr->maxTemp(m_spindex); } void PDSS::initThermo() { AssertThrow(m_tp != 0, "PDSS::initThermo()"); m_vpssmgr_ptr = m_tp->provideVPSSMgr(); initPtrs(); + m_mw = m_tp->molecularWeight(m_spindex); } void PDSS::initAllPtrs(VPStandardStateTP *tp, VPSSMgr *vpssmgr_ptr, diff --git a/Cantera/src/thermo/PDSS_ConstVol.cpp b/Cantera/src/thermo/PDSS_ConstVol.cpp index d568da921..5eff50510 100644 --- a/Cantera/src/thermo/PDSS_ConstVol.cpp +++ b/Cantera/src/thermo/PDSS_ConstVol.cpp @@ -122,7 +122,7 @@ namespace Cantera { m_constMolarVolume = getFloat(*ss, "molarVolume", "toSI"); std::string id = ""; - initThermoXML(phaseNode, id); + // initThermoXML(phaseNode, id); } diff --git a/Cantera/src/thermo/PDSS_HKFT.cpp b/Cantera/src/thermo/PDSS_HKFT.cpp index df6c9165f..cf704673b 100644 --- a/Cantera/src/thermo/PDSS_HKFT.cpp +++ b/Cantera/src/thermo/PDSS_HKFT.cpp @@ -21,7 +21,6 @@ namespace Cantera { PDSS_HKFT::PDSS_HKFT(VPStandardStateTP *tp, int spindex) : PDSS(tp, spindex), m_waterSS(0), - m_pres(OneAtm), m_densWaterSS(-1.0), m_waterProps(0), m_born_coeff_j(-1.0), @@ -43,6 +42,7 @@ namespace Cantera { m_domega_jdT_prtr(0.0), m_charge_j(0.0) { + m_pres = OneAtm; m_pdssType = cPDSS_MOLAL_HKFT; m_presR_bar = OneAtm * 1.0E-5; } @@ -51,7 +51,6 @@ namespace Cantera { PDSS_HKFT::PDSS_HKFT(VPStandardStateTP *tp, int spindex, std::string inputFile, std::string id) : PDSS(tp, spindex), m_waterSS(0), - m_pres(OneAtm), m_densWaterSS(-1.0), m_waterProps(0), m_born_coeff_j(-1.0), @@ -73,6 +72,7 @@ namespace Cantera { m_domega_jdT_prtr(0.0), m_charge_j(0.0) { + m_pres = OneAtm; m_pdssType = cPDSS_MOLAL_HKFT; m_presR_bar = OneAtm * 1.0E-5; constructPDSSFile(tp, spindex, inputFile, id); @@ -82,7 +82,6 @@ namespace Cantera { const XML_Node& phaseRoot, bool spInstalled) : PDSS(tp, spindex), m_waterSS(0), - m_pres(OneAtm), m_densWaterSS(-1.0), m_waterProps(0), m_born_coeff_j(-1.0), @@ -104,6 +103,7 @@ namespace Cantera { m_domega_jdT_prtr(0.0), m_charge_j(0.0) { + m_pres = OneAtm; m_pdssType = cPDSS_MOLAL_HKFT; m_presR_bar = OneAtm * 1.0E-5; // We have to read the info from here @@ -113,7 +113,6 @@ namespace Cantera { PDSS_HKFT::PDSS_HKFT(const PDSS_HKFT &b) : PDSS(b), m_waterSS(0), - m_pres(OneAtm), m_densWaterSS(-1.0), m_waterProps(0), m_born_coeff_j(-1.0), @@ -156,7 +155,6 @@ namespace Cantera { //! Need to call initAllPtrs AFTER, to get the correct m_waterSS m_waterSS = 0; - m_pres = b.m_pres; m_densWaterSS = b.m_densWaterSS; //! Need to call initAllPtrs AFTER, to get the correct m_waterProps if (m_waterProps) { @@ -211,6 +209,14 @@ namespace Cantera { doublereal GG = gibbs_mole(); doublereal SS = entropy_mole(); doublereal h = GG + m_temp * SS; + +#ifdef DEBUG_MODE_NOT + doublereal h2 = enthalpy_mole2(); + if (fabs(h - h2) > 1.0E-1) { + printf("we are here, h = %g, h2 = %g, k = %d, T = %g, P = %g p0 = %g\n", h, h2, m_spindex, m_temp, m_pres, + m_p0); + } +#endif return h; } @@ -221,6 +227,15 @@ namespace Cantera { return hh / RT; } + doublereal + PDSS_HKFT::enthalpy_mole2() const { + doublereal delH = deltaH(); + double enthTRPR = m_Mu0_tr_pr + 298.15*m_Entrop_tr_pr * 1.0E3 * 4.184; + double res = delH + enthTRPR; + return res; + } + + /** * Calculate the internal energy in mks units of * J kmol-1 @@ -262,9 +277,9 @@ namespace Cantera { doublereal c2term = m_c2 / (m_temp - 228.) / (m_temp - 228.); - doublereal a3term = m_a3 / (m_temp - 228.) / (m_temp - 228.) / (m_temp - 228.) * 2.0 * m_temp * (m_pres - OneAtm); + doublereal a3term = -m_a3 / (m_temp - 228.) / (m_temp - 228.) / (m_temp - 228.) * 2.0 * m_temp * (pbar - m_presR_bar); - doublereal a4term = m_a4 / (m_temp - 228.) / (m_temp - 228.) / (m_temp - 228.) * 2.0 * m_temp + doublereal a4term = -m_a4 / (m_temp - 228.) / (m_temp - 228.) / (m_temp - 228.) * 2.0 * m_temp * log((2600. + pbar)/(2600. + m_presR_bar)); doublereal omega_j; @@ -299,8 +314,8 @@ namespace Cantera { domega_jdT = nu * (-(charge2 / r_e_j2 * dr_e_jdT) +(m_charge_j / r_e_H2 * dgvaldT )); - d2omega_jdT2 = nu * ( 2.0*charge2*dr_e_jdT/(r_e_j2*r_e_j) - charge2*d2r_e_jdT2/r_e_j2 - -2.0 *dgvaldT /(r_e_H2*r_e_H) + charge2*d2gvaldT2 /r_e_H2); + d2omega_jdT2 = nu * ( 2.0*charge2*dr_e_jdT*dr_e_jdT/(r_e_j2*r_e_j) - charge2*d2r_e_jdT2/r_e_j2 + -2.0*m_charge_j*dgvaldT*dgvaldT/(r_e_H2*r_e_H) + m_charge_j*d2gvaldT2 /r_e_H2); } doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0); @@ -310,6 +325,15 @@ namespace Cantera { doublereal d2relepsilondT2 = m_waterProps->relEpsilon(m_temp, m_pres, 2); +#ifdef DEBUG_MODE_NOT + doublereal d1 = m_waterProps->relEpsilon(m_temp, m_pres, 1); + doublereal d2 = m_waterProps->relEpsilon(m_temp + 0.0001, m_pres, 1); + doublereal d3 = (d2 - d1) / 0.0001; + if (fabs ( d2relepsilondT2 - d3) > 1.0E-6) { + printf("we are here\n"); + } +#endif + doublereal X = d2relepsilondT2 / (relepsilon* relepsilon) - 2.0 * relepsilon * Y * Y; doublereal Z = -1.0 / relepsilon; @@ -319,11 +343,26 @@ namespace Cantera { doublereal xterm = omega_j * m_temp * X; doublereal otterm = m_temp * d2omega_jdT2 * (Z + 1.0); + + doublereal rterm = - m_domega_jdT_prtr * (m_Z_pr_tr + 1.0); - doublereal Cp_calgmol = c1term + c2term + a3term + a4term + yterm + xterm + otterm; + doublereal Cp_calgmol = c1term + c2term + a3term + a4term + yterm + xterm + otterm + rterm; // Convert to Joules / kmol doublereal Cp = Cp_calgmol * 1.0E3 * 4.184; + + +#ifdef DEBUG_MODE_NOT + double e1 = enthalpy_mole(); + m_temp = m_temp - 0.001; + double e2 = enthalpy_mole(); + m_temp = m_temp + 0.001; + double cpd = (e1 - e2) / 0.001; + if (fabs(Cp - cpd) > 10.0) { + printf("Cp difference : raw: %g, delta: %g, k = %d, T = %g, m_pres = %g\n", + Cp, cpd, m_spindex, m_temp, m_pres); + } +#endif return Cp; } @@ -490,9 +529,7 @@ namespace Cantera { void PDSS_HKFT::initThermo() { 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 @@ -572,7 +609,7 @@ namespace Cantera { void PDSS_HKFT::constructPDSSXML(VPStandardStateTP *tp, int spindex, const XML_Node& speciesNode, const XML_Node& phaseNode, bool spInstalled) { - PDSS::initThermo(); + //PDSS::initThermo(); // m_p0 = OneAtm; @@ -729,7 +766,69 @@ namespace Cantera { delete fxml; } +#ifdef DEBUG_MODE + doublereal PDSS_HKFT::deltaH() const { + + doublereal pbar = m_pres * 1.0E-5; + doublereal c1term = m_c1 * (m_temp - 298.15); + + doublereal a1term = m_a1 * (pbar - m_presR_bar); + + doublereal a2term = m_a2 * log((2600. + pbar)/(2600. + m_presR_bar)); + + doublereal c2term = -m_c2 * ( 1.0/(m_temp - 228.) - 1.0/(298.15 - 228.) ); + + double a3tmp = (2.0 * m_temp - 228.)/ (m_temp - 228.) /(m_temp - 228.); + + doublereal a3term = m_a3 * a3tmp * (pbar - m_presR_bar); + + doublereal a4term = m_a4 * a3tmp * log((2600. + pbar)/(2600. + m_presR_bar)); + + doublereal omega_j; + doublereal domega_jdT; + if (m_charge_j == 0.0) { + omega_j = m_omega_pr_tr; + domega_jdT = 0.0; + } else { + doublereal nu = 166027; + doublereal r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082); + doublereal gval = gstar(m_temp, m_pres, 0); + doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval; + + doublereal dgvaldT = gstar(m_temp, m_pres, 1); + doublereal dr_e_jdT = fabs(m_charge_j) * dgvaldT; + + omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval) ); + + domega_jdT = - nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdT) + + nu * m_charge_j / (3.082 + gval) / (3.082 + gval) * dgvaldT; + } + + doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0); + doublereal drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1); + + doublereal Y = drelepsilondT / (relepsilon * relepsilon); + + doublereal Z = -1.0 / relepsilon; + + doublereal yterm = m_temp * omega_j * Y; + doublereal yrterm = - 298.15 * m_omega_pr_tr * m_Y_pr_tr; + + doublereal wterm = - omega_j * (Z + 1.0); + doublereal wrterm = + m_omega_pr_tr * (m_Z_pr_tr + 1.0); + + doublereal otterm = m_temp * domega_jdT * (Z + 1.0); + doublereal otrterm = - m_temp * m_domega_jdT_prtr * (m_Z_pr_tr + 1.0); + + doublereal deltaH_calgmol = c1term + a1term + a2term + c2term + a3term + a4term + + yterm + yrterm + wterm + wrterm + otterm + otrterm; + + // Convert to Joules / kmol + doublereal deltaH = deltaH_calgmol * 1.0E3 * 4.184; + return deltaH; + } +#endif doublereal PDSS_HKFT::deltaG() const { @@ -854,7 +953,7 @@ namespace Cantera { if (ifunc != 2) { return 0.0; } - return ag_coeff[2] * 2.0;; + return ag_coeff[2] * 2.0; } @@ -961,7 +1060,7 @@ namespace Cantera { + 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; + + gval * alpha * bfunc * dens / ((1.0 - dens) * (1.0 - dens)) * ddensdT; return dfac1dT + dfac2dT + dfac3dT; @@ -981,7 +1080,25 @@ namespace Cantera { doublereal PDSS_HKFT::gstar(const doublereal temp, const doublereal pres, const int ifunc) const { doublereal gval = g(temp, pres, ifunc); doublereal fval = f(temp, pres, ifunc); - return gval - fval; + double res = gval - fval; +#ifdef DEBUG_MODE_NOT + if (ifunc == 2) { + double gval1 = g(temp, pres, 1); + double fval1 = f(temp, pres, 1); + double gval2 = g(temp + 0.001, pres, 1); + double fval2 = f(temp + 0.001, pres, 1); + double gvalT = (gval2 - gval1) / 0.001; + double fvalT = (fval2 - fval1) / 0.001; + if (fabs(gvalT - gval) > 1.0E-9) { + printf("we are here\n"); + } + if (fabs(fvalT - fval) > 1.0E-9) { + printf("we are here\n"); + } + // return gvalT - fvalT; + } +#endif + return res; } diff --git a/Cantera/src/thermo/PDSS_HKFT.h b/Cantera/src/thermo/PDSS_HKFT.h index b54fb244d..feec62513 100644 --- a/Cantera/src/thermo/PDSS_HKFT.h +++ b/Cantera/src/thermo/PDSS_HKFT.h @@ -148,7 +148,9 @@ namespace Cantera { * @return returns the species standard state enthalpy in J kmol-1 */ virtual doublereal enthalpy_mole() const; - +#ifdef DEBUG_MODE + virtual doublereal enthalpy_mole2() const; +#endif //! Return the standard state molar enthalpy divided by RT /*! * Returns the species standard state enthalpy divided by RT at the @@ -486,6 +488,10 @@ namespace Cantera { */ doublereal deltaS() const; +#ifdef DEBUG_MODE + doublereal deltaH() const; +#endif + //! Internal formula for the calculation of a_g() /*! * The output of this is in units of Angstroms @@ -589,9 +595,6 @@ namespace Cantera { */ PDSS_Water *m_waterSS; - //! Current value of the pressure for this object - mutable doublereal m_pres; - //! density of standard-state water /*! * internal temporary variable diff --git a/Cantera/src/thermo/PDSS_IdealGas.cpp b/Cantera/src/thermo/PDSS_IdealGas.cpp index 558305b12..48be63302 100644 --- a/Cantera/src/thermo/PDSS_IdealGas.cpp +++ b/Cantera/src/thermo/PDSS_IdealGas.cpp @@ -110,8 +110,8 @@ namespace Cantera { */ void PDSS_IdealGas::constructPDSSXML(VPStandardStateTP *tp, int spindex, const XML_Node& phaseNode, std::string id) { - initThermo(); - initThermoXML(phaseNode, id); + //initThermo(); + //initThermoXML(phaseNode, id); } @@ -119,13 +119,13 @@ namespace Cantera { std::string inputFile, std::string id) { if (inputFile.size() == 0) { - throw CanteraError("PDSS_IdealGas::initThermo", + throw CanteraError("PDSS_IdealGas::constructPDSSFile", "input file is null"); } std::string path = findInputFile(inputFile); ifstream fin(path.c_str()); if (!fin) { - throw CanteraError("PDSS_IdealGas::initThermo","could not open " + throw CanteraError("PDSS_IdealGas::constructPDSSFile","could not open " +path+" for reading."); } /* @@ -137,7 +137,7 @@ namespace Cantera { fxml->build(fin); XML_Node *fxml_phase = findXMLPhase(fxml, id); if (!fxml_phase) { - throw CanteraError("PDSS_IdealGas::initThermo", + throw CanteraError("PDSS_IdealGas::constructPDSSFile", "ERROR: Can not find phase named " + id + " in file named " + inputFile); } @@ -155,7 +155,6 @@ namespace Cantera { m_p0 = sp.refPressure(m_spindex); m_minTemp = m_spthermo->minTemp(m_spindex); m_maxTemp = m_spthermo->maxTemp(m_spindex); - m_mw = m_tp->molecularWeight(m_spindex); } /** diff --git a/Cantera/src/thermo/SpeciesThermoInterpType.cpp b/Cantera/src/thermo/SpeciesThermoInterpType.cpp index da2ecf503..7b7fa5365 100644 --- a/Cantera/src/thermo/SpeciesThermoInterpType.cpp +++ b/Cantera/src/thermo/SpeciesThermoInterpType.cpp @@ -156,7 +156,8 @@ namespace Cantera { doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const { - m_vpssmgr_ptr->setState_T(temp); + //m_vpssmgr_ptr->setState_T(temp); + m_PDSS_ptr->setTemperature(temp); AssertThrowMsg(m_speciesIndex >= 0, "STITbyPDSS::updatePropertiesTemp", "object was probably not installed correctly"); h_RT[m_speciesIndex] = m_PDSS_ptr->enthalpy_RT_ref(); diff --git a/Cantera/src/thermo/VPSSMgr.cpp b/Cantera/src/thermo/VPSSMgr.cpp index ef8650cba..a53982384 100644 --- a/Cantera/src/thermo/VPSSMgr.cpp +++ b/Cantera/src/thermo/VPSSMgr.cpp @@ -299,12 +299,18 @@ namespace Cantera { } void VPSSMgr::updateStandardStateThermo() { + _updateStandardStateThermo(); } void VPSSMgr::updateRefStateThermo() const { + _updateRefStateThermo(); } void VPSSMgr::_updateStandardStateThermo() { + for (int k = 0; k < m_kk; k++) { + PDSS *kPDSS = m_vptp_ptr->providePDSS(k); + kPDSS->setState_TP(m_tlast, m_plast); + } err("_updateStandardStateThermo()"); } @@ -354,8 +360,21 @@ namespace Cantera { } void VPSSMgr::initThermoXML(XML_Node& phaseNode, std::string id) { - // Add a check to see that all references pressures are the same + const PDSS *kPDSS = m_vptp_ptr->providePDSS(0); + m_p0 = kPDSS->refPressure(); + for (int i = 0; i < m_kk; i++) { + const PDSS *kPDSS = m_vptp_ptr->providePDSS(i); + doublereal mint = kPDSS->minTemp(); + if (mint > m_minTemp) { + m_minTemp = mint; + } + mint = kPDSS->maxTemp(); + if (mint < m_maxTemp) { + m_maxTemp = mint; + } + } #ifdef DEBUG_MODE + // Add a check to see that all references pressures are the same double m_p0_k; if (m_spthermo) { for (int k = 0; k < m_kk; k++) { @@ -385,6 +404,8 @@ namespace Cantera { } } #endif + + } void VPSSMgr::installSTSpecies(int k, const XML_Node& s, @@ -406,14 +427,26 @@ namespace Cantera { /*****************************************************************/ doublereal VPSSMgr::minTemp(int k) const { + if (k >= 0) { + const PDSS *kPDSS = m_vptp_ptr->providePDSS(k); + return kPDSS->minTemp(); + } return m_minTemp; } doublereal VPSSMgr::maxTemp(int k) const { + if (k >= 0) { + const PDSS *kPDSS = m_vptp_ptr->providePDSS(k); + return kPDSS->maxTemp(); + } return m_maxTemp; } - doublereal VPSSMgr::refPressure() const { + doublereal VPSSMgr::refPressure(int k) const { + if (k >= 0) { + const PDSS *kPDSS = m_vptp_ptr->providePDSS(k); + return kPDSS->refPressure(); + } return m_p0; } diff --git a/Cantera/src/thermo/VPSSMgr.h b/Cantera/src/thermo/VPSSMgr.h index 171018bc1..e3b555c80 100644 --- a/Cantera/src/thermo/VPSSMgr.h +++ b/Cantera/src/thermo/VPSSMgr.h @@ -689,7 +689,7 @@ namespace Cantera { * in the same phase have the same reference state pressures. * */ - virtual doublereal refPressure() const ; + virtual doublereal refPressure(int k = -1) const ; //@} @@ -805,7 +805,7 @@ namespace Cantera { * Reference pressure (Pa) must be the same for all species * - defaults to 1 atm. */ - doublereal m_p0; + mutable doublereal m_p0; //! minimum temperature for the standard state calculations doublereal m_minTemp; diff --git a/Cantera/src/thermo/VPSSMgr_Water_ConstVol.cpp b/Cantera/src/thermo/VPSSMgr_Water_ConstVol.cpp index cb42dca4b..fd9793589 100644 --- a/Cantera/src/thermo/VPSSMgr_Water_ConstVol.cpp +++ b/Cantera/src/thermo/VPSSMgr_Water_ConstVol.cpp @@ -87,6 +87,7 @@ namespace Cantera { 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) { doublereal RT = GasConstant * m_tlast; m_waterSS->setState_TP(m_tlast, m_p0); @@ -101,6 +102,7 @@ namespace Cantera { 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) { doublereal RT = GasConstant * m_tlast; m_waterSS->setState_TP(m_tlast, m_p0); @@ -124,6 +126,7 @@ namespace Cantera { 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_s0_R[0] = (m_waterSS->entropy_mole()) / GasConstant; @@ -137,6 +140,7 @@ namespace Cantera { 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_cp0_R[0] = (m_waterSS->cp_mole()) / GasConstant; @@ -150,6 +154,7 @@ namespace Cantera { 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_V0[0] = m_vptp_ptr->molecularWeight(0) / m_waterSS->density(); @@ -160,9 +165,14 @@ namespace Cantera { copy(m_V0.begin(), m_V0.end(), vol); } - - void VPSSMgr_Water_ConstVol::updateRefStateThermo() const { - // Fix up the water + void VPSSMgr_Water_ConstVol::_updateRefStateThermo() const { + m_p0 = m_waterSS->pref_safe(m_tlast); + m_spthermo->update(m_tlast, &m_cp0_R[0], &m_h0_RT[0], &m_s0_R[0]); + for (int k = 0; k < m_kk; k++) { + m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; + PDSS *kPDSS = m_vptp_ptr->providePDSS(k); + kPDSS->setTemperature(m_tlast); + } doublereal RT = GasConstant * m_tlast; m_waterSS->setState_TP(m_tlast, m_p0); m_h0_RT[0] = (m_waterSS->enthalpy_mole())/ RT; @@ -173,9 +183,12 @@ namespace Cantera { m_waterSS->setState_TP(m_tlast, m_plast); } + + void VPSSMgr_Water_ConstVol::_updateStandardStateThermo() { + // _updateRefStateThermo(); doublereal RT = GasConstant * m_tlast; - doublereal del_pRT = (m_plast - m_p0) / (RT); + doublereal del_pRT = (m_plast - OneAtm) / (RT); for (int k = 1; k < m_kk; k++) { m_hss_RT[k] = m_h0_RT[k] + del_pRT * m_Vss[k]; @@ -183,6 +196,8 @@ namespace Cantera { m_sss_R[k] = m_s0_R[k]; m_gss_RT[k] = m_hss_RT[k] - m_sss_R[k]; // m_Vss[k] constant + PDSS *kPDSS = m_vptp_ptr->providePDSS(k); + kPDSS->setState_TP(m_tlast, m_plast); } // Do the water m_waterSS->setState_TP(m_tlast, m_plast); diff --git a/Cantera/src/thermo/VPSSMgr_Water_ConstVol.h b/Cantera/src/thermo/VPSSMgr_Water_ConstVol.h index d4ae22ab5..459cb55ee 100644 --- a/Cantera/src/thermo/VPSSMgr_Water_ConstVol.h +++ b/Cantera/src/thermo/VPSSMgr_Water_ConstVol.h @@ -127,9 +127,17 @@ namespace Cantera { */ virtual void _updateStandardStateThermo(); + //! Updates the reference state thermodynamic functions at the + //! current T of the solution and the reference pressure + /*! + * Underscore updates never check for the state of the system + * They just do the calculation. + */ + virtual void _updateRefStateThermo () const; + public: - //@} + //@} /// @name Thermodynamic Values for the Species Reference States (VPStandardStateTP) /*! * There are also temporary @@ -207,9 +215,7 @@ namespace Cantera { */ virtual void getStandardVolumes_ref(doublereal *vol) const ; - - void updateRefStateThermo() const; - + //! @name Initialization Methods - For Internal use (VPStandardState) /*! * The following methods are used in the process of constructing diff --git a/Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp b/Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp index d5038c10d..b15628dab 100644 --- a/Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp +++ b/Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp @@ -38,7 +38,8 @@ namespace Cantera { VPSSMgr_Water_HKFT::VPSSMgr_Water_HKFT(VPStandardStateTP *vp_ptr, SpeciesThermo *spth) : VPSSMgr(vp_ptr, spth), - m_waterSS(0) + m_waterSS(0), + m_tlastRef(-1.0) { m_useTmpRefStateStorage = true; m_useTmpStandardStateStorage = true; @@ -52,7 +53,8 @@ namespace Cantera { VPSSMgr_Water_HKFT::VPSSMgr_Water_HKFT(const VPSSMgr_Water_HKFT &right) : VPSSMgr(right.m_vptp_ptr, right.m_spthermo), - m_waterSS(0) + m_waterSS(0), + m_tlastRef(-1.0) { m_useTmpRefStateStorage = true; m_useTmpStandardStateStorage = true; @@ -66,6 +68,7 @@ namespace Cantera { if (&b == this) return *this; VPSSMgr::operator=(b); m_waterSS = (PDSS_Water *) m_vptp_ptr->providePDSS(0); + m_tlastRef = -1.0; return *this; } @@ -77,6 +80,7 @@ namespace Cantera { void VPSSMgr_Water_HKFT::getEnthalpy_RT_ref(doublereal *hrt) const{ + updateRefStateThermo(); // Everything should be OK except for the water SS if (m_p0 != m_plast) { doublereal RT = GasConstant * m_tlast; @@ -91,6 +95,7 @@ namespace Cantera { void VPSSMgr_Water_HKFT::getGibbs_RT_ref(doublereal *grt) const{ + updateRefStateThermo(); // Everything should be OK except for the water SS if (m_p0 != m_plast) { doublereal RT = GasConstant * m_tlast; @@ -105,6 +110,7 @@ namespace Cantera { void VPSSMgr_Water_HKFT::getGibbs_ref(doublereal *g) const{ + getGibbs_RT_ref(g); doublereal RT = GasConstant * m_tlast; for (int k = 0; k < m_kk; k++) { g[k] *= RT; @@ -113,6 +119,7 @@ namespace Cantera { void VPSSMgr_Water_HKFT::getEntropy_R_ref(doublereal *sr) const{ + updateRefStateThermo(); // Everything should be OK except for the water SS if (m_p0 != m_plast) { m_waterSS->setState_TP(m_tlast, m_p0); @@ -126,6 +133,7 @@ namespace Cantera { void VPSSMgr_Water_HKFT::getCp_R_ref(doublereal *cpr) const{ + updateRefStateThermo(); // Everything should be OK except for the water SS if (m_p0 != m_plast) { m_waterSS->setState_TP(m_tlast, m_p0); @@ -139,6 +147,7 @@ namespace Cantera { void VPSSMgr_Water_HKFT::getStandardVolumes_ref(doublereal *vol) const{ + updateRefStateThermo(); // Everything should be OK except for the water SS if (m_p0 != m_plast) { m_waterSS->setState_TP(m_tlast, m_p0); @@ -150,26 +159,65 @@ namespace Cantera { copy(m_V0.begin(), m_V0.end(), vol); } + void VPSSMgr_Water_HKFT::setState_P(doublereal pres) { + if (m_plast != pres) { + m_plast = pres; + _updateStandardStateThermo(); + } + } + + void VPSSMgr_Water_HKFT::setState_T(doublereal temp) { + if (m_tlast != temp) { + m_tlast = temp; + _updateStandardStateThermo(); + } + } + + void VPSSMgr_Water_HKFT::setState_TP(doublereal temp, doublereal pres) { + if (m_tlast != temp) { + m_tlast = temp; + m_plast = pres; + _updateStandardStateThermo(); + } else if (m_plast != pres) { + m_plast = pres; + _updateStandardStateThermo(); + } + } void VPSSMgr_Water_HKFT::updateRefStateThermo() const { + if (m_tlastRef != m_tlast) { + m_tlastRef = m_tlast; + _updateRefStateThermo(); + } + } + + void VPSSMgr_Water_HKFT::_updateRefStateThermo() const { // Fix up the water + m_p0 = m_waterSS->pref_safe(m_tlast); doublereal RT = GasConstant * m_tlast; m_waterSS->setState_TP(m_tlast, m_p0); m_h0_RT[0] = (m_waterSS->enthalpy_mole())/ RT; 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_waterSS->density()) / m_vptp_ptr->molecularWeight(0); + m_V0[0] = (m_waterSS->density()) / m_vptp_ptr->molecularWeight(0); m_waterSS->setState_TP(m_tlast, m_plast); for (int k = 1; k < m_kk; k++) { PDSS_HKFT *ps = (PDSS_HKFT *) m_vptp_ptr->providePDSS(k); ps->setState_TP(m_tlast, m_p0); - m_cpss_R[k] = ps->cp_R(); - m_sss_R[k] = ps->entropy_mole(); - m_gss_RT[k] = ps->gibbs_RT();; - m_hss_RT[k] = m_gss_RT[k] + m_sss_R[k]; - m_Vss[k] = ps->molarVolume(); + m_cp0_R[k] = ps->cp_R(); + m_s0_R[k] = ps->entropy_mole() / GasConstant; + m_g0_RT[k] = ps->gibbs_RT(); + + m_h0_RT[k] = m_g0_RT[k] + m_s0_R[k]; +#ifdef DEBUG_MODE + double h = ps->enthalpy_RT(); + if (fabs( m_h0_RT[k] - h) > 1.0E-4) { + printf("we are here\n"); + } +#endif + m_V0[k] = ps->molarVolume(); } } diff --git a/Cantera/src/thermo/VPSSMgr_Water_HKFT.h b/Cantera/src/thermo/VPSSMgr_Water_HKFT.h index f8c437135..c4f059aa4 100644 --- a/Cantera/src/thermo/VPSSMgr_Water_HKFT.h +++ b/Cantera/src/thermo/VPSSMgr_Water_HKFT.h @@ -179,7 +179,28 @@ namespace Cantera { */ virtual void getStandardVolumes_ref(doublereal *vol) const ; - + //! Set the temperature (K) and pressure (Pa) + /*! + * This sets the temperature and pressure and triggers + * calculation of underlying quantities + * + * @param T Temperature (K) + * @param P Pressure (Pa) + */ + virtual void setState_TP(doublereal T, doublereal P); + + //! Set the temperature (K) + /*! + * @param T Temperature (K) + */ + virtual void setState_T(doublereal T); + + //! Set the pressure (Pa) + /*! + * @param P Pressure (Pa) + */ + virtual void setState_P(doublereal P); + //@} /// @name Setting the Internal State of the System /*! @@ -210,8 +231,8 @@ namespace Cantera { */ virtual void updateRefStateThermo() const; - protected: - + protected: + virtual void _updateRefStateThermo() const; //! Updates the standard state thermodynamic functions at the current T and P of the solution. /*! * @internal @@ -239,6 +260,8 @@ namespace Cantera { */ virtual void _updateStandardStateThermo(); + + public: //@} @@ -327,6 +350,8 @@ namespace Cantera { //! Shallow pointer to the water object PDSS_Water *m_waterSS; + + mutable doublereal m_tlastRef; }; //@} } diff --git a/Cantera/src/thermo/VPStandardStateTP.cpp b/Cantera/src/thermo/VPStandardStateTP.cpp index 6518b9c19..5404615e9 100644 --- a/Cantera/src/thermo/VPStandardStateTP.cpp +++ b/Cantera/src/thermo/VPStandardStateTP.cpp @@ -411,7 +411,7 @@ namespace Cantera { void VPStandardStateTP::initThermoXML(XML_Node& phaseNode, std::string id) { VPStandardStateTP::initLengths(); - m_VPSS_ptr->initThermo(); + //m_VPSS_ptr->initThermo(); for (int k = 0; k < m_kk; k++) { PDSS *kPDSS = m_PDSS_storage[k]; AssertTrace(kPDSS != 0);