research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Fixed various errors.

Browse source 
Fixed PDSS_HKFT's Cp calculation. It now is correct
   took out out of place initThermo and initThermoXML calls.
This commit is contained in:
Harry Moffat 2009-01-04 06:34:19 +00:00
parent a879666984
commit 76ebcf3726
13 changed files with 304 additions and 54 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
Cantera/src/thermo/PDSS.cpp
View file

@ -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,

2
Cantera/src/thermo/PDSS_ConstVol.cpp
View file

@ -122,7 +122,7 @@ namespace Cantera {
m_constMolarVolume = getFloat(*ss, "molarVolume", "toSI");
std::string id = "";
initThermoXML(phaseNode, id);
// initThermoXML(phaseNode, id);
}

151
Cantera/src/thermo/PDSS_HKFT.cpp
View file

@ -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;
}

11
Cantera/src/thermo/PDSS_HKFT.h
View file

@ -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

11
Cantera/src/thermo/PDSS_IdealGas.cpp
View file

@ -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);
}
/**

3
Cantera/src/thermo/SpeciesThermoInterpType.cpp
View file

@ -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();

37
Cantera/src/thermo/VPSSMgr.cpp
View file

@ -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;
}

4
Cantera/src/thermo/VPSSMgr.h
View file

@ -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;

23
Cantera/src/thermo/VPSSMgr_Water_ConstVol.cpp
View file

@ -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);

14
Cantera/src/thermo/VPSSMgr_Water_ConstVol.h
View file

@ -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

64
Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp
View file

@ -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();
}
}

31
Cantera/src/thermo/VPSSMgr_Water_HKFT.h
View file

@ -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;
};
//@}
}

2
Cantera/src/thermo/VPStandardStateTP.cpp
View file

@ -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);