Use ThermoPhase::RT() where convenient

This commit is contained in:
Ray Speth 2015-11-12 15:36:17 -05:00
parent 0d0c96257b
commit 5e69a9949a
22 changed files with 61 additions and 70 deletions

View file

@ -346,7 +346,7 @@ public:
* \see SpeciesThermo
*/
virtual doublereal enthalpy_mole() const {
return GasConstant * temperature() * mean_X(enthalpy_RT_ref());
return RT() * mean_X(enthalpy_RT_ref());
}
/**
@ -404,7 +404,7 @@ public:
* @param p Pressure (Pa)
*/
virtual void setPressure(doublereal p) {
setDensity(p * meanMolecularWeight() / (GasConstant * temperature()));
setDensity(p * meanMolecularWeight() / RT());
}
//! Set the density and pressure at constant composition.

View file

@ -57,7 +57,7 @@ void AqueousKinetics::updateKc()
// compute Delta G^0 for all reversible reactions
getRevReactionDelta(m_grt.data(), m_rkcn.data());
doublereal rrt = 1.0/(GasConstant * thermo().temperature());
doublereal rrt = 1.0 / thermo().RT();
for (size_t i = 0; i < m_revindex.size(); i++) {
size_t irxn = m_revindex[i];
m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt);
@ -82,7 +82,7 @@ void AqueousKinetics::getEquilibriumConstants(doublereal* kc)
// compute Delta G^0 for all reactions
getReactionDelta(m_grt.data(), m_rkcn.data());
doublereal rrt = 1.0/(GasConstant * thermo().temperature());
doublereal rrt = 1.0 / thermo().RT();
for (size_t i = 0; i < nReactions(); i++) {
kc[i] = exp(-m_rkcn[i]*rrt);
}

View file

@ -103,7 +103,7 @@ void GasKinetics::updateKc()
// compute Delta G^0 for all reversible reactions
getRevReactionDelta(m_grt.data(), m_rkcn.data());
doublereal rrt = 1.0/(GasConstant * thermo().temperature());
doublereal rrt = 1.0 / thermo().RT();
for (size_t i = 0; i < m_revindex.size(); i++) {
size_t irxn = m_revindex[i];
m_rkcn[irxn] = std::min(exp(m_rkcn[irxn]*rrt - m_dn[irxn]*m_logStandConc),
@ -124,7 +124,7 @@ void GasKinetics::getEquilibriumConstants(doublereal* kc)
// compute Delta G^0 for all reactions
getReactionDelta(m_grt.data(), m_rkcn.data());
doublereal rrt = 1.0/(GasConstant * thermo().temperature());
doublereal rrt = 1.0 / thermo().RT();
for (size_t i = 0; i < nReactions(); i++) {
kc[i] = exp(-m_rkcn[i]*rrt + m_dn[i]*m_logStandConc);
}

View file

@ -205,7 +205,7 @@ void InterfaceKinetics::updateKc()
* and m_mu0_Kc[]
*/
updateMu0();
doublereal rrt = 1.0 / (GasConstant * thermo(0).temperature());
doublereal rrt = 1.0 / thermo(0).RT();
// compute Delta mu^0 for all reversible reactions
getRevReactionDelta(m_mu0_Kc.data(), m_rkcn.data());
@ -281,7 +281,7 @@ void InterfaceKinetics::checkPartialEquil()
void InterfaceKinetics::getEquilibriumConstants(doublereal* kc)
{
updateMu0();
doublereal rrt = 1.0 / (GasConstant * thermo(0).temperature());
doublereal rrt = 1.0 / thermo(0).RT();
std::fill(kc, kc + nReactions(), 0.0);
getReactionDelta(m_mu0_Kc.data(), kc);
for (size_t i = 0; i < nReactions(); i++) {

View file

@ -47,8 +47,7 @@ int ConstDensityThermo::eosType() const
doublereal ConstDensityThermo::enthalpy_mole() const
{
doublereal p0 = m_spthermo->refPressure();
return GasConstant * temperature() *
mean_X(enthalpy_RT()) + (pressure() - p0)/molarDensity();
return RT() * mean_X(enthalpy_RT()) + (pressure() - p0)/molarDensity();
}
doublereal ConstDensityThermo::entropy_mole() const

View file

@ -240,7 +240,7 @@ void FixedChemPotSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_
chemPot_ = val;
} else {
_updateThermo();
chemPot_ = (m_h0_RT[0] - m_s0_R[0]) * GasConstant * temperature();
chemPot_ = (m_h0_RT[0] - m_s0_R[0]) * RT();
}
}

View file

@ -83,7 +83,7 @@ doublereal IdealGasPhase::cv_mole() const
doublereal IdealGasPhase::standardConcentration(size_t k) const
{
return pressure() / (GasConstant * temperature());
return pressure() / RT();
}
void IdealGasPhase::getActivityCoefficients(doublereal* ac) const
@ -97,8 +97,7 @@ void IdealGasPhase::getStandardChemPotentials(doublereal* muStar) const
{
const vector_fp& gibbsrt = gibbs_RT_ref();
scale(gibbsrt.begin(), gibbsrt.end(), muStar, RT());
double tmp = log(pressure() / m_spthermo->refPressure());
tmp *= GasConstant * temperature();
double tmp = log(pressure() / m_spthermo->refPressure()) * RT();
for (size_t k = 0; k < m_kk; k++) {
muStar[k] += tmp; // add RT*ln(P/P_0)
}
@ -317,7 +316,7 @@ void IdealGasPhase::_updateThermo() const
for (size_t k = 0; k < m_kk; k++) {
m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
}
m_logc0 = log(m_p0 / (GasConstant * tnow));
m_logc0 = log(m_p0 / RT());
}
}
}

View file

@ -112,7 +112,7 @@ int IdealSolidSolnPhase::eosType() const
doublereal IdealSolidSolnPhase::enthalpy_mole() const
{
doublereal htp = GasConstant * temperature() * mean_X(enthalpy_RT_ref());
doublereal htp = RT() * mean_X(enthalpy_RT_ref());
return htp + (pressure() - m_Pref)/molarDensity();
}
@ -123,7 +123,7 @@ doublereal IdealSolidSolnPhase::entropy_mole() const
doublereal IdealSolidSolnPhase::gibbs_mole() const
{
return GasConstant * temperature() * (mean_X(gibbs_RT_ref()) + sum_xlogx());
return RT() * (mean_X(gibbs_RT_ref()) + sum_xlogx());
}
doublereal IdealSolidSolnPhase::cp_mole() const
@ -304,7 +304,7 @@ void IdealSolidSolnPhase::getChemPotentials_RT(doublereal* mu) const
void IdealSolidSolnPhase::getPartialMolarEnthalpies(doublereal* hbar) const
{
const vector_fp& _h = enthalpy_RT_ref();
scale(_h.begin(), _h.end(), hbar, GasConstant * temperature());
scale(_h.begin(), _h.end(), hbar, RT());
}
void IdealSolidSolnPhase::getPartialMolarEntropies(doublereal* sbar) const
@ -367,7 +367,7 @@ void IdealSolidSolnPhase::getEntropy_R(doublereal* sr) const
void IdealSolidSolnPhase::getIntEnergy_RT(doublereal* urt) const
{
const vector_fp& _h = enthalpy_RT_ref();
doublereal prefrt = m_Pref / (GasConstant * temperature());
doublereal prefrt = m_Pref / RT();
for (size_t k = 0; k < m_kk; k++) {
urt[k] = _h[k] - prefrt * m_speciesMolarVolume[k];
}
@ -405,7 +405,7 @@ void IdealSolidSolnPhase::getGibbs_RT_ref(doublereal* grt) const
void IdealSolidSolnPhase::getGibbs_ref(doublereal* g) const
{
_updateThermo();
double tmp = GasConstant * temperature();
double tmp = RT();
for (size_t k = 0; k != m_kk; k++) {
g[k] = tmp * m_g0_RT[k];
}
@ -414,7 +414,7 @@ void IdealSolidSolnPhase::getGibbs_ref(doublereal* g) const
void IdealSolidSolnPhase::getIntEnergy_RT_ref(doublereal* urt) const
{
const vector_fp& _h = enthalpy_RT_ref();
doublereal prefrt = m_Pref / (GasConstant * temperature());
doublereal prefrt = m_Pref / RT();
for (size_t k = 0; k < m_kk; k++) {
urt[k] = _h[k] - prefrt * m_speciesMolarVolume[k];
}
@ -566,7 +566,7 @@ void IdealSolidSolnPhase::_updateThermo() const
// Update the thermodynamic functions of the reference state.
m_spthermo->update(tnow, m_cp0_R.data(), m_h0_RT.data(), m_s0_R.data());
m_tlast = tnow;
doublereal rrt = 1.0 / (GasConstant * tnow);
doublereal rrt = 1.0 / RT();
for (size_t k = 0; k < m_kk; k++) {
double deltaE = rrt * m_pe[k];
m_h0_RT[k] += deltaE;

View file

@ -84,7 +84,7 @@ int IdealSolnGasVPSS::eosType() const
doublereal IdealSolnGasVPSS::enthalpy_mole() const
{
updateStandardStateThermo();
return GasConstant * temperature() * mean_X(m_VPSS_ptr->enthalpy_RT());
return RT() * mean_X(m_VPSS_ptr->enthalpy_RT());
}
doublereal IdealSolnGasVPSS::entropy_mole() const
@ -115,8 +115,7 @@ void IdealSolnGasVPSS::calcDensity()
{
// Calculate the molarVolume of the solution (m**3 kmol-1)
if (m_idealGas) {
double dens = (m_Pcurrent * meanMolecularWeight()
/(GasConstant * temperature()));
double dens = m_Pcurrent * meanMolecularWeight() / RT();
Phase::setDensity(dens);
} else {
const doublereal* const dtmp = moleFractdivMMW();
@ -168,8 +167,7 @@ void IdealSolnGasVPSS::getActivityConcentrations(doublereal* c) const
doublereal IdealSolnGasVPSS::standardConcentration(size_t k) const
{
if (m_idealGas) {
double p = pressure();
return p/(GasConstant * temperature());
return pressure() / RT();
} else {
const vector_fp& vss = m_VPSS_ptr->getStandardVolumes();
switch (m_formGC) {
@ -205,7 +203,7 @@ void IdealSolnGasVPSS::getChemPotentials(doublereal* mu) const
void IdealSolnGasVPSS::getPartialMolarEnthalpies(doublereal* hbar) const
{
getEnthalpy_RT(hbar);
scale(hbar, hbar+m_kk, hbar, GasConstant * temperature());
scale(hbar, hbar+m_kk, hbar, RT());
}
void IdealSolnGasVPSS::getPartialMolarEntropies(doublereal* sbar) const
@ -221,7 +219,7 @@ void IdealSolnGasVPSS::getPartialMolarEntropies(doublereal* sbar) const
void IdealSolnGasVPSS::getPartialMolarIntEnergies(doublereal* ubar) const
{
getIntEnergy_RT(ubar);
scale(ubar, ubar+m_kk, ubar, GasConstant * temperature());
scale(ubar, ubar+m_kk, ubar, RT());
}
void IdealSolnGasVPSS::getPartialMolarCp(doublereal* cpbar) const

View file

@ -288,15 +288,13 @@ void IonsFromNeutralVPSSTP::getChemPotentials(doublereal* mu) const
// Get the standard chemical potentials of netural molecules
neutralMoleculePhase_->getStandardChemPotentials(muNeutralMolecule_.data());
doublereal RT_ = GasConstant * temperature();
switch (ionSolnType_) {
case cIonSolnType_PASSTHROUGH:
neutralMoleculePhase_->getChemPotentials(mu);
break;
case cIonSolnType_SINGLEANION:
neutralMoleculePhase_->getLnActivityCoefficients(lnActCoeff_NeutralMolecule_.data());
fact2 = 2.0 * RT_ * log(2.0);
fact2 = 2.0 * RT() * log(2.0);
// Do the cation list
for (size_t k = 0; k < cationList_.size(); k++) {
@ -304,21 +302,21 @@ void IonsFromNeutralVPSSTP::getChemPotentials(doublereal* mu) const
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
xx = std::max(SmallNumber, moleFractions_[icat]);
mu[icat] = muNeutralMolecule_[jNeut] + fact2 + RT_ * (lnActCoeff_NeutralMolecule_[jNeut] + log(xx));
mu[icat] = muNeutralMolecule_[jNeut] + fact2 + RT() * (lnActCoeff_NeutralMolecule_[jNeut] + log(xx));
}
// Do the anion list
icat = anionList_[0];
jNeut = fm_invert_ionForNeutral[icat];
xx = std::max(SmallNumber, moleFractions_[icat]);
mu[icat] = RT_ * log(xx);
mu[icat] = RT() * log(xx);
// Do the list of neutral molecules
for (size_t k = 0; k < passThroughList_.size(); k++) {
icat = passThroughList_[k];
jNeut = fm_invert_ionForNeutral[icat];
xx = std::max(SmallNumber, moleFractions_[icat]);
mu[icat] = muNeutralMolecule_[jNeut] + RT_ * (lnActCoeff_NeutralMolecule_[jNeut] + log(xx));
mu[icat] = muNeutralMolecule_[jNeut] + RT() * (lnActCoeff_NeutralMolecule_[jNeut] + log(xx));
}
break;

View file

@ -65,7 +65,7 @@ ThermoPhase* LatticePhase::duplMyselfAsThermoPhase() const
doublereal LatticePhase::enthalpy_mole() const
{
return GasConstant * temperature() * mean_X(enthalpy_RT_ref()) +
return RT() * mean_X(enthalpy_RT_ref()) +
(pressure() - m_Pref)/molarDensity();
}
@ -162,7 +162,7 @@ void LatticePhase::getChemPotentials(doublereal* mu) const
void LatticePhase::getPartialMolarEnthalpies(doublereal* hbar) const
{
const vector_fp& _h = enthalpy_RT_ref();
scale(_h.begin(), _h.end(), hbar, GasConstant * temperature());
scale(_h.begin(), _h.end(), hbar, RT());
}
void LatticePhase::getPartialMolarEntropies(doublereal* sbar) const
@ -230,7 +230,7 @@ void LatticePhase::getGibbs_ref(doublereal* g) const
{
getGibbs_RT_ref(g);
for (size_t k = 0; k < m_kk; k++) {
g[k] *= GasConstant * temperature();
g[k] *= RT();
}
}

View file

@ -320,7 +320,7 @@ void LatticeSolidPhase::getGibbs_ref(doublereal* g) const
{
getGibbs_RT_ref(g);
for (size_t k = 0; k < m_kk; k++) {
g[k] *= GasConstant * temperature();
g[k] *= RT();
}
}

View file

@ -71,7 +71,7 @@ void MaskellSolidSolnPhase::getActivityConcentrations(doublereal* c) const
doublereal MaskellSolidSolnPhase::enthalpy_mole() const
{
_updateThermo();
const doublereal h0 = GasConstant * temperature() * mean_X(m_h0_RT);
const doublereal h0 = RT() * mean_X(m_h0_RT);
const doublereal r = moleFraction(product_species_index);
const doublereal fmval = fm(r);
return h0 + r * fmval * h_mixing;
@ -145,7 +145,7 @@ void MaskellSolidSolnPhase::getActivityCoefficients(doublereal* ac) const
const doublereal rfm = r * fm(r);
const doublereal A = (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval)) /
(std::pow(1 - r - rfm, 1 + pval) * (1 - r));
const doublereal B = pval * h_mixing / (GasConstant * temperature());
const doublereal B = pval * h_mixing / RT();
cached.value[product_species_index] = A * std::exp(B);
cached.value[reactant_species_index] = 1 / (A * r * (1-r) ) * std::exp(-B);
}
@ -159,7 +159,7 @@ void MaskellSolidSolnPhase::getChemPotentials(doublereal* mu) const
const doublereal pval = p(r);
const doublereal rfm = r * fm(r);
const doublereal DgbarDr = pval * h_mixing +
GasConstant * temperature() *
RT() *
std::log( (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval) * r) /
(std::pow(1 - r - rfm, 1 + pval) * (1 - r)) );
mu[product_species_index] = RT() * m_g0_RT[product_species_index] + DgbarDr;
@ -199,9 +199,8 @@ void MaskellSolidSolnPhase::getPartialMolarVolumes(doublereal* vbar) const
void MaskellSolidSolnPhase::getPureGibbs(doublereal* gpure) const
{
_updateThermo();
const doublereal RT = GasConstant * temperature();
for (size_t sp=0; sp < m_kk; ++sp) {
gpure[sp] = RT * m_g0_RT[sp];
gpure[sp] = RT() * m_g0_RT[sp];
}
}
@ -285,7 +284,7 @@ void MaskellSolidSolnPhase::_updateThermo() const
doublereal MaskellSolidSolnPhase::s() const
{
return 1 + std::exp(h_mixing / (GasConstant * temperature()));
return 1 + std::exp(h_mixing / RT());
}
doublereal MaskellSolidSolnPhase::fm(const doublereal r) const

View file

@ -106,7 +106,7 @@ doublereal MetalSHEelectrons::logStandardConc(size_t k) const
void MetalSHEelectrons::getStandardChemPotentials(doublereal* mu0) const
{
getGibbs_RT(mu0);
mu0[0] *= GasConstant * temperature();
mu0[0] *= RT();
}
void MetalSHEelectrons::getEnthalpy_RT(doublereal* hrt) const

View file

@ -115,7 +115,7 @@ doublereal MineralEQ3::logStandardConc(size_t k) const
void MineralEQ3::getStandardChemPotentials(doublereal* mu0) const
{
getGibbs_RT(mu0);
mu0[0] *= GasConstant * temperature();
mu0[0] *= RT();
}
void MineralEQ3::getEnthalpy_RT(doublereal* hrt) const

View file

@ -249,7 +249,7 @@ void PureFluidPhase::getGibbs_RT_ref(doublereal* grt) const
void PureFluidPhase::getGibbs_ref(doublereal* g) const
{
getGibbs_RT_ref(g);
g[0] *= (GasConstant * temperature());
g[0] *= RT();
}
void PureFluidPhase::getEntropy_R_ref(doublereal* er) const

View file

@ -16,8 +16,8 @@ SemiconductorPhase::SemiconductorPhase(std::string infile,
void SemiconductorPhase::getChemPotentials(doublereal* mu) const
{
getActivityConcentrations(m_work.data());
mu[0] = ec() + GasConstant*temperature()*(JoyceDixon(m_work[0]/nc()));
mu[1] = ev() + GasConstant*temperature()*(log(m_work[1]/nv()));
mu[0] = ec() + RT()*(JoyceDixon(m_work[0]/nc()));
mu[1] = ev() + RT()*(log(m_work[1]/nv()));
}
// units: kmol/m^3

View file

@ -128,7 +128,7 @@ void SingleSpeciesTP::getChemPotentials(doublereal* mu) const
void SingleSpeciesTP::getChemPotentials_RT(doublereal* murt) const
{
getStandardChemPotentials(murt);
murt[0] /= GasConstant * temperature();
murt[0] /= RT();
}
void SingleSpeciesTP::getElectrochemPotentials(doublereal* mu) const
@ -139,13 +139,13 @@ void SingleSpeciesTP::getElectrochemPotentials(doublereal* mu) const
void SingleSpeciesTP::getPartialMolarEnthalpies(doublereal* hbar) const
{
getEnthalpy_RT(hbar);
hbar[0] *= GasConstant * temperature();
hbar[0] *= RT();
}
void SingleSpeciesTP::getPartialMolarIntEnergies(doublereal* ubar) const
{
getIntEnergy_RT(ubar);
ubar[0] *= GasConstant * temperature();
ubar[0] *= RT();
}
void SingleSpeciesTP::getPartialMolarEntropies(doublereal* sbar) const
@ -170,7 +170,7 @@ void SingleSpeciesTP::getPartialMolarVolumes(doublereal* vbar) const
void SingleSpeciesTP::getPureGibbs(doublereal* gpure) const
{
getGibbs_RT(gpure);
gpure[0] *= GasConstant * temperature();
gpure[0] *= RT();
}
void SingleSpeciesTP::getStandardVolumes(doublereal* vbar) const
@ -195,7 +195,7 @@ void SingleSpeciesTP::getGibbs_RT_ref(doublereal* grt) const
void SingleSpeciesTP::getGibbs_ref(doublereal* g) const
{
getGibbs_RT_ref(g);
g[0] *= GasConstant * temperature();
g[0] *= RT();
}
void SingleSpeciesTP::getEntropy_R_ref(doublereal* er) const

View file

@ -102,14 +102,14 @@ doublereal StoichSubstance::logStandardConc(size_t k) const
void StoichSubstance::getStandardChemPotentials(doublereal* mu0) const
{
getGibbs_RT(mu0);
mu0[0] *= GasConstant * temperature();
mu0[0] *= RT();
}
void StoichSubstance::getEnthalpy_RT(doublereal* hrt) const
{
getEnthalpy_RT_ref(hrt);
doublereal presCorrect = (m_press - m_p0) / molarDensity();
hrt[0] += presCorrect / (GasConstant * temperature());
hrt[0] += presCorrect / RT();
}
void StoichSubstance::getEntropy_R(doublereal* sr) const
@ -132,7 +132,7 @@ void StoichSubstance::getCp_R(doublereal* cpr) const
void StoichSubstance::getIntEnergy_RT(doublereal* urt) const
{
_updateThermo();
urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / (GasConstant * temperature());
urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / RT();
}
// ---- Thermodynamic Values for the Species Reference States ----
@ -140,7 +140,7 @@ void StoichSubstance::getIntEnergy_RT(doublereal* urt) const
void StoichSubstance::getIntEnergy_RT_ref(doublereal* urt) const
{
_updateThermo();
urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / (GasConstant * temperature());
urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / RT();
}
// ---- Initialization and Internal functions

View file

@ -146,8 +146,7 @@ void SurfPhase::getChemPotentials(doublereal* mu) const
copy(m_mu0.begin(), m_mu0.end(), mu);
getActivityConcentrations(m_work.data());
for (size_t k = 0; k < m_kk; k++) {
mu[k] += GasConstant * temperature() *
(log(m_work[k]) - logStandardConc(k));
mu[k] += RT() * (log(m_work[k]) - logStandardConc(k));
}
}
@ -184,13 +183,13 @@ void SurfPhase::getPureGibbs(doublereal* g) const
void SurfPhase::getGibbs_RT(doublereal* grt) const
{
_updateThermo();
scale(m_mu0.begin(), m_mu0.end(), grt, 1.0/(GasConstant*temperature()));
scale(m_mu0.begin(), m_mu0.end(), grt, 1.0/RT());
}
void SurfPhase::getEnthalpy_RT(doublereal* hrt) const
{
_updateThermo();
scale(m_h0.begin(), m_h0.end(), hrt, 1.0/(GasConstant*temperature()));
scale(m_h0.begin(), m_h0.end(), hrt, 1.0/RT());
}
void SurfPhase::getEntropy_R(doublereal* sr) const

View file

@ -804,14 +804,14 @@ void ThermoPhase::setElementPotentials(const vector_fp& lambda)
if (!m_hasElementPotentials) {
m_lambdaRRT.resize(mm);
}
scale(lambda.begin(), lambda.end(), m_lambdaRRT.begin(), 1.0/(GasConstant* temperature()));
scale(lambda.begin(), lambda.end(), m_lambdaRRT.begin(), 1.0/RT());
m_hasElementPotentials = true;
}
bool ThermoPhase::getElementPotentials(doublereal* lambda) const
{
if (m_hasElementPotentials) {
scale(m_lambdaRRT.begin(), m_lambdaRRT.end(), lambda, GasConstant* temperature());
scale(m_lambdaRRT.begin(), m_lambdaRRT.end(), lambda, RT());
}
return m_hasElementPotentials;
}

View file

@ -154,7 +154,7 @@ void WaterSSTP::setParametersFromXML(const XML_Node& eosdata)
void WaterSSTP::getEnthalpy_RT(doublereal* hrt) const
{
*hrt = (m_sub.enthalpy() + EW_Offset)/(GasConstant*temperature());
*hrt = (m_sub.enthalpy() + EW_Offset) / RT();
}
void WaterSSTP::getIntEnergy_RT(doublereal* ubar) const
@ -169,8 +169,7 @@ void WaterSSTP::getEntropy_R(doublereal* sr) const
void WaterSSTP::getGibbs_RT(doublereal* grt) const
{
double T = temperature();
*grt = (m_sub.Gibbs() + EW_Offset - SW_Offset*T) / (GasConstant * T);
*grt = (m_sub.Gibbs() + EW_Offset) / RT() - SW_Offset / GasConstant;
if (!m_ready) {
throw CanteraError("waterSSTP::", "Phase not ready");
}
@ -209,7 +208,7 @@ void WaterSSTP::getEnthalpy_RT_ref(doublereal* hrt) const
throw CanteraError("setPressure", "error");
}
doublereal h = m_sub.enthalpy();
*hrt = (h + EW_Offset) / (GasConstant * T);
*hrt = (h + EW_Offset) / RT();
dd = m_sub.density(T, p, waterState, dens);
}
@ -229,7 +228,7 @@ void WaterSSTP::getGibbs_RT_ref(doublereal* grt) const
}
m_sub.setState_TR(T, dd);
doublereal g = m_sub.Gibbs();
*grt = (g + EW_Offset - SW_Offset*T)/ (GasConstant * T);
*grt = (g + EW_Offset - SW_Offset*T)/ RT();
dd = m_sub.density(T, p, waterState, dens);
}