diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp index 88b729c3f..94e9f9e66 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp @@ -1390,12 +1390,11 @@ namespace VCSnonideal { * Now, calculate a sample naught gibbs free energy calculation * at the specified temperature. */ - VolPhase->G0_calc(vprob->T); double R = vcsUtil_gasConstant(vprob->m_VCS_UnitsFormat); for (k = 0; k < nSpPhase; k++) { vcs_SpeciesProperties *sProp = VolPhase->ListSpeciesPtr[k]; ts_ptr = sProp->SpeciesThermo; - ts_ptr->SS0_feSave = VolPhase->SS0ChemicalPotential[k] / R; + ts_ptr->SS0_feSave = VolPhase->G0_calc_one(k)/ R; ts_ptr->SS0_TSave = vprob->T; } diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index dff1edc6c..0859f671d 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -62,6 +62,7 @@ namespace VCSnonideal { m_UpToDate_VolStar(false), m_UpToDate_VolPM(false), m_UpToDate_GStar(false), + m_UpToDate_G0(false), Temp(273.15), Pres(1.01325E5), RefPres(1.01325E5) @@ -124,6 +125,7 @@ namespace VCSnonideal { m_UpToDate_VolStar(false), m_UpToDate_VolPM(false), m_UpToDate_GStar(false), + m_UpToDate_G0(false), Temp(b.Temp), Pres(b.Pres) { @@ -238,6 +240,7 @@ namespace VCSnonideal { m_UpToDate_VolStar = false; m_UpToDate_VolPM = false; m_UpToDate_GStar = false; + m_UpToDate_G0 = false; Temp = b.Temp; Pres = b.Pres; setState_TP(Temp, Pres); @@ -299,7 +302,7 @@ namespace VCSnonideal { IndSpecies.resize(nspecies,-1); - if ((int) ListSpeciesPtr.size() >= NVolSpecies) { + if ((int) ListSpeciesPtr.size() >= NVolSpecies) { for (int i = 0; i < NVolSpecies; i++) { if (ListSpeciesPtr[i]) { delete ListSpeciesPtr[i]; @@ -332,6 +335,7 @@ namespace VCSnonideal { m_UpToDate_VolStar = false; m_UpToDate_VolPM = false; m_UpToDate_GStar = false; + m_UpToDate_G0 = false; } /************************************************************************************/ @@ -384,31 +388,21 @@ namespace VCSnonideal { // Gibbs free energy calculation at a temperature for the reference state // of each species /* - * @param TKelvin temperature */ - void vcs_VolPhase::G0_calc(double tkelvin) { - bool lsame = false; - if (Temp == tkelvin) { - lsame = true; - } - - bool doit = !lsame; - setState_TP(tkelvin, Pres); - if (SS0ChemicalPotential[0] == -1) doit = true; - if (doit) { - if (m_useCanteraCalls) { - TP_ptr->getGibbs_ref(VCS_DATA_PTR(SS0ChemicalPotential)); - } else { - double R = vcsUtil_gasConstant(m_VCS_UnitsFormat); - for (int k = 0; k < NVolSpecies; k++) { - int kglob = IndSpecies[k]; - vcs_SpeciesProperties *sProp = ListSpeciesPtr[k]; - VCS_SPECIES_THERMO *sTherm = sProp->SpeciesThermo; - SS0ChemicalPotential[k] = - R * (sTherm->G0_R_calc(kglob, tkelvin)); - } + void vcs_VolPhase::_updateG0() const { + if (m_useCanteraCalls) { + TP_ptr->getGibbs_ref(VCS_DATA_PTR(SS0ChemicalPotential)); + } else { + double R = vcsUtil_gasConstant(m_VCS_UnitsFormat); + for (int k = 0; k < NVolSpecies; k++) { + int kglob = IndSpecies[k]; + vcs_SpeciesProperties *sProp = ListSpeciesPtr[k]; + VCS_SPECIES_THERMO *sTherm = sProp->SpeciesThermo; + SS0ChemicalPotential[k] = + R * (sTherm->G0_R_calc(kglob, Temp)); } } + m_UpToDate_G0 = true; } /*******************************************************************************/ @@ -420,8 +414,10 @@ namespace VCSnonideal { * * @return return value of the gibbs free energy */ - double vcs_VolPhase::G0_calc_one(int kspec, double tkelvin) { - G0_calc(tkelvin); + double vcs_VolPhase::G0_calc_one(int kspec) const { + if (!m_UpToDate_G0) { + _updateG0(); + } return SS0ChemicalPotential[kspec]; } /*******************************************************************************/ @@ -463,7 +459,7 @@ namespace VCSnonideal { * @return Gstar[kspec] returns the gibbs free energy for the * standard state of the kspec species. */ - double vcs_VolPhase::GStar_calc_one(int kspec) { + double vcs_VolPhase::GStar_calc_one(int kspec) const { if (!m_UpToDate_GStar) { _updateGStar(); } @@ -741,7 +737,7 @@ namespace VCSnonideal { * in all of the phases in a VCS problem. Only the * entries for the current phase are filled in. */ - void vcs_VolPhase::sendToVCS_GStar(double * const gstar){ + void vcs_VolPhase::sendToVCS_GStar(double * const gstar) const { if (!m_UpToDate_GStar) { _updateGStar(); } @@ -782,7 +778,7 @@ namespace VCSnonideal { * @param temperature_Kelvin (Kelvin) * @param pressure_PA Pressure (MKS units - Pascal) */ - void vcs_VolPhase::setState_TP(double temp, double pres) + void vcs_VolPhase::setState_TP(const double temp, const double pres) { if (Temp == temp) { if (Pres == pres) { @@ -799,6 +795,21 @@ namespace VCSnonideal { m_UpToDate_VolStar = false; m_UpToDate_VolPM = false; m_UpToDate_GStar = false; + m_UpToDate_G0 = false; + } + /****************************************************************************/ + + // Sets the temperature in this object and + // underlying objects + /* + * Sets the temperature and pressure in this object and + * underlying objects. The underlying objects refers to the + * Cantera's ThermoPhase object for this phase. + * + * @param temperature_Kelvin (Kelvin) + */ + void vcs_VolPhase::setState_T(const double temp) { + setState_TP(temp, Pres); } /**************************************************************************/ @@ -841,7 +852,7 @@ namespace VCSnonideal { * state */ double vcs_VolPhase::VolStar_calc_one(int kspec) const { - if (!m_UpToDate_VolStar) { + if (!m_UpToDate_VolStar) { _updateVolStar(); } return StarMolarVol[kspec]; diff --git a/Cantera/src/equil/vcs_VolPhase.h b/Cantera/src/equil/vcs_VolPhase.h index ffec88d12..a2acabdd4 100644 --- a/Cantera/src/equil/vcs_VolPhase.h +++ b/Cantera/src/equil/vcs_VolPhase.h @@ -245,14 +245,9 @@ namespace VCSnonideal { * @return Gstar[kspec] returns the gibbs free energy for the * standard state of the kth species. */ - double GStar_calc_one(int kspec); + double GStar_calc_one(int kspec) const; - //! Gibbs free energy calculation at a temperature for the reference state - //! of each species - /*! - * @param TKelvin temperature - */ - void G0_calc(double TKelvin); + //! Gibbs free energy calculation at a temperature for the reference state //! of a species, return a value for one species @@ -262,7 +257,7 @@ namespace VCSnonideal { * * @return return value of the gibbs free energy */ - double G0_calc_one(int kspec, double TKelvin); + double G0_calc_one(int kspec) const; //! Molar volume calculation for standard state of one species /*! @@ -301,7 +296,7 @@ namespace VCSnonideal { * in all of the phases in a VCS problem. Only the * entries for the current phase are filled in. */ - void sendToVCS_GStar(double * const gstar); + void sendToVCS_GStar(double * const gstar) const; //! Sets the temperature and pressure in this object and //! underlying objects @@ -313,7 +308,18 @@ namespace VCSnonideal { * @param temperature_Kelvin (Kelvin) * @param pressure_PA Pressure (MKS units - Pascal) */ - void setState_TP(double temperature_Kelvin, double pressure_PA); + void setState_TP(const double temperature_Kelvin, const double pressure_PA); + + //! Sets the temperature in this object and + //! underlying objects + /*! + * Sets the temperature and pressure in this object and + * underlying objects. The underlying objects refers to the + * Cantera's ThermoPhase object for this phase. + * + * @param temperature_Kelvin (Kelvin) + */ + void setState_T(const double temperature_Kelvin); // Downloads the ln ActCoeff jacobian into the VCS version of the // ln ActCoeff jacobian. @@ -417,6 +423,13 @@ namespace VCSnonideal { */ void _updateGStar() const; + //! Gibbs free energy calculation at a temperature for the reference state + //! of each species + /*! + * + */ + void _updateG0() const; + //! Molar volume calculation for standard states /*! * Calculate the molar volume for the standard states @@ -703,7 +716,8 @@ namespace VCSnonideal { * units are m**3 */ mutable double m_totalVol; - + + private: //! Vector of calculated SS0 chemical potentials for the //! current Temperature. /*! @@ -716,7 +730,6 @@ namespace VCSnonideal { */ mutable std::vector SS0ChemicalPotential; - private: //! Vector of calculated Star chemical potentials for the //! current Temperature and pressure. /*! @@ -806,6 +819,13 @@ namespace VCSnonideal { */ mutable bool m_UpToDate_GStar; + + //! Boolean indicating whether G0 is uptodate. + /*! + * G0 is sensitive to the temperature and the pressure, only + */ + mutable bool m_UpToDate_G0; + //! Current value of the temperature for this object, and underlying objects double Temp; diff --git a/Cantera/src/equil/vcs_prob.cpp b/Cantera/src/equil/vcs_prob.cpp index 3c0088e13..611ddb261 100644 --- a/Cantera/src/equil/vcs_prob.cpp +++ b/Cantera/src/equil/vcs_prob.cpp @@ -320,7 +320,6 @@ namespace VCSnonideal { " SS0ChemPot StarChemPot\n"); for (iphase = 0; iphase < NPhase; iphase++) { Vphase = VPhaseList[iphase]; - Vphase->G0_calc(T); Vphase->setState_TP(T, PresPA); for (int kindex = 0; kindex < Vphase->NVolSpecies; kindex++) { int kglob = Vphase->IndSpecies[kindex]; @@ -331,8 +330,7 @@ namespace VCSnonideal { plogf(" "); } - plogf("%16g %16g\n", - Vphase->SS0ChemicalPotential[kindex], + plogf("%16g %16g\n", Vphase->G0_calc_one(kindex), Vphase->GStar_calc_one(kindex)); } } diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 05a417920..758c71979 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -970,7 +970,7 @@ namespace VCSnonideal { pubPhase->setMoleFractions(VCS_DATA_PTR(vPhase->moleFractions())); for (int k = 0; k < pubPhase->NVolSpecies; k++) { kT = pubPhase->IndSpecies[k]; - pubPhase->SS0ChemicalPotential[k] = vPhase->SS0ChemicalPotential[k]; + //pubPhase->SS0ChemicalPotential[k] = vPhase->SS0ChemicalPotential[k]; pubPhase->StarMolarVol[k] = vPhase->StarMolarVol[k]; pubPhase->PartialMolarVol[k] = vPhase->PartialMolarVol[k]; pubPhase->ActCoeff[k] = vPhase->ActCoeff[k]; diff --git a/Cantera/src/equil/vcs_species_thermo.cpp b/Cantera/src/equil/vcs_species_thermo.cpp index 469567d50..d84fdd979 100644 --- a/Cantera/src/equil/vcs_species_thermo.cpp +++ b/Cantera/src/equil/vcs_species_thermo.cpp @@ -280,7 +280,8 @@ double VCS_SPECIES_THERMO::G0_R_calc(int kglob, double TKelvin) if (UseCanteraCalls) { AssertThrowVCS(m_VCS_UnitsFormat == VCS_UNITS_MKS, "Possible inconsistency"); int kspec = IndexSpeciesPhase; - fe = OwningPhase->G0_calc_one(kspec, TKelvin); + OwningPhase->setState_T(TKelvin); + fe = OwningPhase->G0_calc_one(kspec); double R = vcsUtil_gasConstant(m_VCS_UnitsFormat); fe /= R; } else {