Doxygen update -> no real changes.

Working on doxygen docs and SAND report and validation for
  EQ3 import.
This commit is contained in:
Harry Moffat 2008-09-29 16:05:40 +00:00
parent 0d846144f0
commit 54988b2b5f
6 changed files with 339 additions and 220 deletions

View file

@ -179,13 +179,19 @@ namespace Cantera {
*/
void addElement(const XML_Node& e);
//! Adde an element, checking for uniqueness
//! Add an element, checking for uniqueness
/*!
* The uniqueness is checked by comparing the string symbol. If
* not unique, nothing is done.
*
* @param symbol String symbol of the element
* @param weight Atomic weight of the element (kg kmol-1).
* @param atomicNumber Atomic number of the element (unitless)
* @param entropy298 Entropy of the element at 298 K and 1 bar
* in its most stable form. The default is
* the value ENTROPY298_UNKNOWN, which is
* interpreted as an unknown, and if used
* will cause Cantera to throw an error.
*/
void addUniqueElement(const std::string& symbol, doublereal weight,
int atomicNumber = 0,

View file

@ -154,16 +154,16 @@ namespace Cantera {
doublereal
PDSS_HKFT::enthalpy_mole() const {
// Ok we may change this evaluation method in the future.
double GG = gibbs_mole();
double SS = entropy_mole();
double h = GG + m_temp * SS;
doublereal GG = gibbs_mole();
doublereal SS = entropy_mole();
doublereal h = GG + m_temp * SS;
return h;
}
doublereal
PDSS_HKFT::enthalpy_RT() const {
double hh = enthalpy_mole();
double RT = GasConstant * m_temp;
doublereal hh = enthalpy_mole();
doublereal RT = GasConstant * m_temp;
return hh / RT;
}
@ -173,8 +173,8 @@ namespace Cantera {
*/
doublereal
PDSS_HKFT::intEnergy_mole() const {
double hh = enthalpy_RT();
double mv = molarVolume();
doublereal hh = enthalpy_RT();
doublereal mv = molarVolume();
return (hh - mv * m_pres);
}
@ -184,7 +184,7 @@ namespace Cantera {
*/
doublereal
PDSS_HKFT::entropy_mole() const {
double delS = deltaS();
doublereal delS = deltaS();
return (m_Entrop_tr_pr * 1.0E3 * 4.184 + delS);
}
@ -193,7 +193,7 @@ namespace Cantera {
* J kmol-1
*/
doublereal PDSS_HKFT::gibbs_mole() const {
double delG = deltaG();
doublereal delG = deltaG();
return (m_Mu0_tr_pr + delG);
}
@ -203,56 +203,56 @@ namespace Cantera {
*/
doublereal PDSS_HKFT::cp_mole() const {
double pbar = m_pres * 1.0E-5;
double c1term = m_c1;
doublereal pbar = m_pres * 1.0E-5;
doublereal c1term = m_c1;
double c2term = m_c2 / (m_temp - 228.) / (m_temp - 228.);
doublereal c2term = m_c2 / (m_temp - 228.) / (m_temp - 228.);
double 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 * (m_pres - OneAtm);
double 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));
double nu = 166027;
double r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
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);
double gval = gstar(m_temp, m_pres, 0);
doublereal gval = gstar(m_temp, m_pres, 0);
double dgvaldT = gstar(m_temp, m_pres, 1);
double d2gvaldT2 = gstar(m_temp, m_pres, 2);
doublereal dgvaldT = gstar(m_temp, m_pres, 1);
doublereal d2gvaldT2 = gstar(m_temp, m_pres, 2);
double r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
double dr_e_jdT = fabs(m_charge_j) * dgvaldT;
doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
doublereal dr_e_jdT = fabs(m_charge_j) * dgvaldT;
double omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval) );
doublereal omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval) );
double domega_jdT = - 2.0 * nu * (m_charge_j * m_charge_j * m_charge_j * m_charge_j / (r_e_j * r_e_j* r_e_j)
doublereal domega_jdT = - 2.0 * nu * (m_charge_j * m_charge_j * m_charge_j * m_charge_j / (r_e_j * r_e_j* r_e_j)
- m_charge_j / (3.082 + gval) / (3.082 + gval) / (3.082 + gval)) * dgvaldT * dgvaldT
- nu * (m_charge_j * m_charge_j * fabs(m_charge_j) / (r_e_j * r_e_j)
- m_charge_j / (3.082 + gval) / (3.082 + gval)) * d2gvaldT2;
double d2omega_jdT2 = nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdT)
doublereal d2omega_jdT2 = 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;
double relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
double drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
doublereal drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
double Y = drelepsilondT / (relepsilon * relepsilon);
doublereal Y = drelepsilondT / (relepsilon * relepsilon);
double d2relepsilondT2 = m_waterProps->relEpsilon(m_temp, m_pres, 2);
doublereal d2relepsilondT2 = m_waterProps->relEpsilon(m_temp, m_pres, 2);
double X = d2relepsilondT2 / (relepsilon* relepsilon) - 2.0 * relepsilon * Y * Y;
doublereal X = d2relepsilondT2 / (relepsilon* relepsilon) - 2.0 * relepsilon * Y * Y;
double Z = -1.0 / relepsilon;
doublereal Z = -1.0 / relepsilon;
double yterm = 2.0 * m_temp * Y * domega_jdT;
doublereal yterm = 2.0 * m_temp * Y * domega_jdT;
double xterm = omega_j * m_temp * X;
doublereal xterm = omega_j * m_temp * X;
double otterm = m_temp * d2omega_jdT2 * (Z + 1.0);
doublereal otterm = m_temp * d2omega_jdT2 * (Z + 1.0);
double Cp_calgmol = c1term + c2term + a3term + a4term + yterm + xterm + otterm;
doublereal Cp_calgmol = c1term + c2term + a3term + a4term + yterm + xterm + otterm;
// Convert to Joules / kmol
doublereal Cp = Cp_calgmol * 1.0E3 * 4.184;
@ -272,97 +272,97 @@ namespace Cantera {
doublereal
PDSS_HKFT::molarVolume() const {
// double pbar = m_pres * 1.0E-5;
// doublereal pbar = m_pres * 1.0E-5;
double a1term = m_a1 * 1.0E-5;
doublereal a1term = m_a1 * 1.0E-5;
double a2term = m_a2 / (2600.E5 + m_pres);
doublereal a2term = m_a2 / (2600.E5 + m_pres);
double a3term = m_a3 * 1.0E-5/ (m_temp - 228.);
doublereal a3term = m_a3 * 1.0E-5/ (m_temp - 228.);
double a4term = m_a4 / (m_temp - 228.) / (2600.E5 + m_pres);
doublereal a4term = m_a4 / (m_temp - 228.) / (2600.E5 + m_pres);
double nu = 166027.;
double r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
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);
double gval = gstar(m_temp, m_pres, 0);
double dgvaldP = gstar(m_temp, m_pres, 3);
doublereal gval = gstar(m_temp, m_pres, 0);
doublereal dgvaldP = gstar(m_temp, m_pres, 3);
double r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
double omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval) );
doublereal omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval) );
double relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
double dr_e_jdP = fabs(m_charge_j) * dgvaldP;
doublereal dr_e_jdP = fabs(m_charge_j) * dgvaldP;
double domega_jdP = - nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdP)
doublereal domega_jdP = - nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdP)
+ nu * m_charge_j / (3.082 + gval) / (3.082 + gval) * dgvaldP;
double drelepsilondP = m_waterProps->relEpsilon(m_temp, m_pres, 3);
doublereal drelepsilondP = m_waterProps->relEpsilon(m_temp, m_pres, 3);
double Q = drelepsilondP / (relepsilon * relepsilon);
doublereal Q = drelepsilondP / (relepsilon * relepsilon);
double Z = -1.0 / relepsilon;
doublereal Z = -1.0 / relepsilon;
double wterm = - domega_jdP * (Z + 1.0);
doublereal wterm = - domega_jdP * (Z + 1.0);
double qterm = - omega_j * Q;
doublereal qterm = - omega_j * Q;
double molVol_calgmolPascal = a1term + a2term + a3term + a4term + wterm + qterm;
doublereal molVol_calgmolPascal = a1term + a2term + a3term + a4term + wterm + qterm;
// Convert to m**3 / kmol
double molVol = molVol_calgmolPascal * 4.184 * 1.0E3;
doublereal molVol = molVol_calgmolPascal * 4.184 * 1.0E3;
return molVol;
}
doublereal
PDSS_HKFT::density() const {
double val = molarVolume();
doublereal val = molarVolume();
return (m_mw/val);
}
doublereal
PDSS_HKFT::gibbs_RT_ref() const {
double m_psave = m_pres;
doublereal m_psave = m_pres;
m_pres = OneAtm;
double ee = gibbs_RT();
doublereal ee = gibbs_RT();
m_pres = m_psave;
return ee;
}
doublereal
PDSS_HKFT::enthalpy_RT_ref() const {
double m_psave = m_pres;
doublereal m_psave = m_pres;
m_pres = OneAtm;
double hh = enthalpy_RT();
doublereal hh = enthalpy_RT();
m_pres = m_psave;
return hh;
}
doublereal
PDSS_HKFT::entropy_R_ref() const {
double m_psave = m_pres;
doublereal m_psave = m_pres;
m_pres = OneAtm;
double ee = entropy_R();
doublereal ee = entropy_R();
m_pres = m_psave;
return ee;
}
doublereal
PDSS_HKFT::cp_R_ref() const {
double m_psave = m_pres;
doublereal m_psave = m_pres;
m_pres = OneAtm;
double ee = cp_R();
doublereal ee = cp_R();
m_pres = m_psave;
return ee;
}
doublereal
PDSS_HKFT::molarVolume_ref() const {
double m_psave = m_pres;
doublereal m_psave = m_pres;
m_pres = OneAtm;
double ee = molarVolume();
doublereal ee = molarVolume();
m_pres = m_psave;
return ee;
}
@ -426,15 +426,15 @@ namespace Cantera {
*/
m_temp = 273.15 + 25.;
m_pres = OneAtm;
double relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
m_waterSS->setState_TP(m_temp, m_pres);
m_densWaterSS = m_waterSS->density();
m_Z_pr_tr = -1.0 / relepsilon;
//double m_Z_pr_tr = -0.0127803;
//doublereal m_Z_pr_tr = -0.0127803;
//printf("m_Z_pr_tr = %20.10g\n", m_Z_pr_tr );
double drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
//double m_Y_pr_tr = -5.799E-5;
doublereal drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
//doublereal m_Y_pr_tr = -5.799E-5;
m_Y_pr_tr = drelepsilondT / (relepsilon * relepsilon);
//printf("m_Y_pr_tr = %20.10g\n", m_Y_pr_tr );
@ -447,9 +447,9 @@ namespace Cantera {
//! Ok, we have mu. Let's check it against the input value
// of DH_F to see that we have some internal consistency
double Hcalc = m_Mu0_tr_pr + 298.15 * (m_Entrop_tr_pr * 1.0E3 * 4.184);
doublereal Hcalc = m_Mu0_tr_pr + 298.15 * (m_Entrop_tr_pr * 1.0E3 * 4.184);
double DHjmol = m_deltaH_formation_tr_pr * 1.0E3 * 4.184;
doublereal DHjmol = m_deltaH_formation_tr_pr * 1.0E3 * 4.184;
// If the discrepency is greater than 100 cal gmol-1, print
// an error and exit.
@ -458,15 +458,15 @@ namespace Cantera {
"DHjmol is not consistent with G and S" +
fp2str(Hcalc) + " vs " + fp2str(DHjmol));
}
double nu = 166027;
double r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
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);
double gval = gstar(m_temp, m_pres, 0);
doublereal gval = gstar(m_temp, m_pres, 0);
double dgvaldT = gstar(m_temp, m_pres, 1);
doublereal dgvaldT = gstar(m_temp, m_pres, 1);
double r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
double dr_e_jdT = fabs(m_charge_j) * dgvaldT;
doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
doublereal dr_e_jdT = fabs(m_charge_j) * dgvaldT;
m_domega_jdT_prtr = - nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdT)
@ -526,21 +526,21 @@ namespace Cantera {
}
if (hh->hasChild("DG0_f_Pr_Tr")) {
double val = getFloat(*hh, "DG0_f_Pr_Tr");
doublereal val = getFloat(*hh, "DG0_f_Pr_Tr");
m_deltaG_formation_tr_pr = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing DG0_f_Pr_Tr field");
}
if (hh->hasChild("DH0_f_Pr_Tr")) {
double val = getFloat(*hh, "DH0_f_Pr_Tr");
doublereal val = getFloat(*hh, "DH0_f_Pr_Tr");
m_deltaH_formation_tr_pr = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing DH0_f_Pr_Tr field");
}
if (hh->hasChild("S0_Pr_Tr")) {
double val = getFloat(*hh, "S0_Pr_Tr");
doublereal val = getFloat(*hh, "S0_Pr_Tr");
m_Entrop_tr_pr= val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing S0_Pr_Tr field");
@ -558,44 +558,44 @@ namespace Cantera {
+ speciesNode.name());
}
if (ss->hasChild("a1")) {
double val = getFloat(*ss, "a1");
doublereal val = getFloat(*ss, "a1");
m_a1 = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a1 field");
}
if (ss->hasChild("a2")) {
double val = getFloat(*ss, "a2");
doublereal val = getFloat(*ss, "a2");
m_a2 = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a2 field");
}
if (ss->hasChild("a3")) {
double val = getFloat(*ss, "a3");
doublereal val = getFloat(*ss, "a3");
m_a3 = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a3 field");
}
if (ss->hasChild("a4")) {
double val = getFloat(*ss, "a4");
doublereal val = getFloat(*ss, "a4");
m_a4 = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a4 field");
}
if (ss->hasChild("c1")) {
double val = getFloat(*ss, "c1");
doublereal val = getFloat(*ss, "c1");
m_c1 = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing c1 field");
}
if (ss->hasChild("c2")) {
double val = getFloat(*ss, "c2");
doublereal val = getFloat(*ss, "c2");
m_c2 = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing c2 field");
}
if (ss->hasChild("omega_Pr_Tr")) {
double val = getFloat(*ss, "omega_Pr_Tr");
doublereal val = getFloat(*ss, "omega_Pr_Tr");
m_omega_pr_tr = val;
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing omega_Pr_Tr field");
@ -645,117 +645,112 @@ namespace Cantera {
double PDSS_HKFT::deltaG() const {
doublereal PDSS_HKFT::deltaG() const {
double pbar = m_pres * 1.0E-5;
//double m_presR_bar = OneAtm * 1.0E-5;
doublereal pbar = m_pres * 1.0E-5;
//doublereal m_presR_bar = OneAtm * 1.0E-5;
double sterm = - m_Entrop_tr_pr * (m_temp - 298.15);
doublereal sterm = - m_Entrop_tr_pr * (m_temp - 298.15);
double c1term = -m_c1 * (m_temp * log(m_temp/298.15) - (m_temp - 298.15));
double a1term = m_a1 * (pbar - m_presR_bar);
doublereal c1term = -m_c1 * (m_temp * log(m_temp/298.15) - (m_temp - 298.15));
doublereal a1term = m_a1 * (pbar - m_presR_bar);
double a2term = m_a2 * log((2600. + pbar)/(2600. + m_presR_bar));
doublereal a2term = m_a2 * log((2600. + pbar)/(2600. + m_presR_bar));
double c2term = -m_c2 * (( 1.0/(m_temp - 228.) - 1.0/(298.15 - 228.) ) * (228. - m_temp)/228.
doublereal c2term = -m_c2 * (( 1.0/(m_temp - 228.) - 1.0/(298.15 - 228.) ) * (228. - m_temp)/228.
- m_temp / (228.*228.) * log( (298.15*(m_temp-228.)) / (m_temp*(298.15-228.)) ));
double a3term = m_a3 / (m_temp - 228.) * (pbar - m_presR_bar);
doublereal a3term = m_a3 / (m_temp - 228.) * (pbar - m_presR_bar);
double a4term = m_a4 / (m_temp - 228.) * log((2600. + pbar)/(2600. + m_presR_bar));
doublereal a4term = m_a4 / (m_temp - 228.) * log((2600. + pbar)/(2600. + m_presR_bar));
double nu = 166027;
double r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
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);
double gval = gstar(m_temp, m_pres, 0);
doublereal gval = gstar(m_temp, m_pres, 0);
double r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
double omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval) );
doublereal omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval) );
double relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
double Z = -1.0 / relepsilon;
doublereal Z = -1.0 / relepsilon;
double wterm = - omega_j * (Z + 1.0);
doublereal wterm = - omega_j * (Z + 1.0);
double wrterm = m_omega_pr_tr * (m_Z_pr_tr + 1.0);
doublereal wrterm = m_omega_pr_tr * (m_Z_pr_tr + 1.0);
double yterm = m_omega_pr_tr * m_Y_pr_tr * (m_temp - 298.15);
doublereal yterm = m_omega_pr_tr * m_Y_pr_tr * (m_temp - 298.15);
double deltaG_calgmol = sterm + c1term + a1term + a2term + c2term + a3term + a4term + wterm + wrterm + yterm;
doublereal deltaG_calgmol = sterm + c1term + a1term + a2term + c2term + a3term + a4term + wterm + wrterm + yterm;
// Convert to Joules / kmol
double deltaG = deltaG_calgmol * 1.0E3 * 4.184;
doublereal deltaG = deltaG_calgmol * 1.0E3 * 4.184;
return deltaG;
}
double PDSS_HKFT::deltaS() const {
doublereal PDSS_HKFT::deltaS() const {
double pbar = m_pres * 1.0E-5;
doublereal pbar = m_pres * 1.0E-5;
double c1term = m_c1 * log(m_temp/298.15);
doublereal c1term = m_c1 * log(m_temp/298.15);
double c2term = -m_c2 / 228. * (( 1.0/(m_temp - 228.) - 1.0/(298.15 - 228.) )
doublereal c2term = -m_c2 / 228. * (( 1.0/(m_temp - 228.) - 1.0/(298.15 - 228.) )
+ 1.0 / 228. * log( (298.15*(m_temp-228.)) / (m_temp*(298.15-228.)) ));
double a3term = m_a3 / (m_temp - 228.) / (m_temp - 228.) * (pbar - m_presR_bar);
doublereal a3term = m_a3 / (m_temp - 228.) / (m_temp - 228.) * (pbar - m_presR_bar);
double a4term = m_a4 / (m_temp - 228.) / (m_temp - 228.) * log((2600. + pbar)/(2600. + m_presR_bar));
doublereal a4term = m_a4 / (m_temp - 228.) / (m_temp - 228.) * log((2600. + pbar)/(2600. + m_presR_bar));
double nu = 166027;
double r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
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);
double gval = gstar(m_temp, m_pres, 0);
doublereal gval = gstar(m_temp, m_pres, 0);
double dgvaldT = gstar(m_temp, m_pres, 1);
doublereal dgvaldT = gstar(m_temp, m_pres, 1);
double r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
double dr_e_jdT = fabs(m_charge_j) * dgvaldT;
doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
doublereal dr_e_jdT = fabs(m_charge_j) * dgvaldT;
double omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval) );
doublereal omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval) );
double domega_jdT = - nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdT)
doublereal 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;
double relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
double drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
doublereal drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
double Y = drelepsilondT / (relepsilon * relepsilon);
doublereal Y = drelepsilondT / (relepsilon * relepsilon);
double Z = -1.0 / relepsilon;
doublereal Z = -1.0 / relepsilon;
double wterm = omega_j * Y;
doublereal wterm = omega_j * Y;
double wrterm = - m_omega_pr_tr * m_Y_pr_tr;
doublereal wrterm = - m_omega_pr_tr * m_Y_pr_tr;
double otterm = domega_jdT * (Z + 1.0);
doublereal otterm = domega_jdT * (Z + 1.0);
double otrterm = - m_domega_jdT_prtr * (m_Z_pr_tr + 1.0);
doublereal otrterm = - m_domega_jdT_prtr * (m_Z_pr_tr + 1.0);
double deltaS_calgmol = c1term + c2term + a3term + a4term + wterm + wrterm + otterm + otrterm;
doublereal deltaS_calgmol = c1term + c2term + a3term + a4term + wterm + wrterm + otterm + otrterm;
// Convert to Joules / kmol
double deltaS = deltaS_calgmol * 1.0E3 * 4.184;
doublereal deltaS = deltaS_calgmol * 1.0E3 * 4.184;
return deltaS;
}
double PDSS_HKFT::electrostatic_radii_calc() {
return 0.0;
}
//! Internal formula for the calculation of a_g()
// Internal formula for the calculation of a_g()
/*
* The output of this is in units of Angstroms
*/
double PDSS_HKFT::ag(const double temp, const int ifunc) const {
static double ag_coeff[3] = { -2.037662, 5.747000E-3, -6.557892E-6};
doublereal PDSS_HKFT::ag(const doublereal temp, const int ifunc) const {
static doublereal ag_coeff[3] = { -2.037662, 5.747000E-3, -6.557892E-6};
if (ifunc == 0) {
double t2 = temp * temp;
double val = ag_coeff[0] + ag_coeff[1] * temp + ag_coeff[2] * t2;
doublereal t2 = temp * temp;
doublereal val = ag_coeff[0] + ag_coeff[1] * temp + ag_coeff[2] * t2;
return val;
} else if (ifunc == 1) {
return ag_coeff[1] + ag_coeff[2] * 2.0 * temp;
@ -767,15 +762,15 @@ namespace Cantera {
}
//! Internal formula for the calculation of b_g()
// Internal formula for the calculation of b_g()
/*
* the output of this is unitless
*/
double PDSS_HKFT::bg(const double temp, const int ifunc) const {
static double bg_coeff[3] = { 6.107361, -1.074377E-2, 1.268348E-5};
doublereal PDSS_HKFT::bg(const doublereal temp, const int ifunc) const {
static doublereal bg_coeff[3] = { 6.107361, -1.074377E-2, 1.268348E-5};
if (ifunc == 0) {
double t2 = temp * temp;
double val = bg_coeff[0] + bg_coeff[1] * temp + bg_coeff[2] * t2;
doublereal t2 = temp * temp;
doublereal val = bg_coeff[0] + bg_coeff[1] * temp + bg_coeff[2] * t2;
return val;
} else if (ifunc == 1) {
return bg_coeff[1] + bg_coeff[2] * 2.0 * temp;
@ -787,24 +782,24 @@ namespace Cantera {
}
double PDSS_HKFT::f(const double temp, const double pres, const int ifunc) const {
doublereal PDSS_HKFT::f(const doublereal temp, const doublereal pres, const int ifunc) const {
static double af_coeff[3] = { 3.666666E1, -0.1504956E-9, 0.5107997E-13};
double TC = temp - 273.15;
double presBar = pres / 1.0E5;
static doublereal af_coeff[3] = { 3.666666E1, -0.1504956E-9, 0.5107997E-13};
doublereal TC = temp - 273.15;
doublereal presBar = pres / 1.0E5;
if (TC < 155.0) return 0.0;
if (TC > 355.0) TC = 355.0;
if (presBar > 1000.) return 0.0;
double T1 = (TC-155.0)/300.;
double fac1;
doublereal T1 = (TC-155.0)/300.;
doublereal fac1;
double p2 = (1000. - presBar) * (1000. - presBar);
double p3 = (1000. - presBar) * p2;
double p4 = p2 * p2;
double fac2 = af_coeff[1] * p3 + af_coeff[2] * p4;
doublereal p2 = (1000. - presBar) * (1000. - presBar);
doublereal p3 = (1000. - presBar) * p2;
doublereal p4 = p2 * p2;
doublereal fac2 = af_coeff[1] * p3 + af_coeff[2] * p4;
if (ifunc == 0) {
fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
return fac1 * fac2;
@ -825,14 +820,14 @@ namespace Cantera {
}
double PDSS_HKFT::g(const double temp, const double pres, const int ifunc) const {
double afunc = ag(temp, 0);
double bfunc = bg(temp, 0);
doublereal PDSS_HKFT::g(const doublereal temp, const doublereal pres, const int ifunc) const {
doublereal afunc = ag(temp, 0);
doublereal bfunc = bg(temp, 0);
m_waterSS->setState_TP(temp, pres);
m_densWaterSS = m_waterSS->density();
// density in gm cm-3
double dens = m_densWaterSS * 1.0E-3;
double gval = afunc * pow((1.0-dens), bfunc);
doublereal dens = m_densWaterSS * 1.0E-3;
doublereal gval = afunc * pow((1.0-dens), bfunc);
if (dens >= 1.0) {
return 0.0;
}
@ -840,33 +835,33 @@ namespace Cantera {
return gval;
} else if (ifunc == 1 || ifunc == 2) {
double afuncdT = ag(temp, 1);
double bfuncdT = bg(temp, 1);
double alpha = m_waterSS->thermalExpansionCoeff();
doublereal afuncdT = ag(temp, 1);
doublereal bfuncdT = bg(temp, 1);
doublereal alpha = m_waterSS->thermalExpansionCoeff();
double fac1 = afuncdT * gval / afunc;
double fac2 = bfuncdT * gval * log(1.0 - dens);
double fac3 = gval * alpha * bfunc * dens / (1.0 - dens);
doublereal fac1 = afuncdT * gval / afunc;
doublereal fac2 = bfuncdT * gval * log(1.0 - dens);
doublereal fac3 = gval * alpha * bfunc * dens / (1.0 - dens);
double dgdt = fac1 + fac2 + fac3;
doublereal dgdt = fac1 + fac2 + fac3;
if (ifunc == 1) {
return dgdt;
}
double afuncdT2 = ag(temp, 2);
double bfuncdT2 = bg(temp, 2);
doublereal afuncdT2 = ag(temp, 2);
doublereal bfuncdT2 = bg(temp, 2);
double dfac1dT = dgdt * afuncdT / afunc + afuncdT2 * gval / afunc
doublereal dfac1dT = dgdt * afuncdT / afunc + afuncdT2 * gval / afunc
- afuncdT * afuncdT * gval / (afunc * afunc);
double ddensdT = - alpha * dens;
double dfac2dT = bfuncdT2 * gval * log(1.0 - dens)
doublereal ddensdT = - alpha * dens;
doublereal dfac2dT = bfuncdT2 * gval * log(1.0 - dens)
+ bfuncdT * dgdt * log(1.0 - dens)
- bfuncdT * gval /(1.0 - dens) * ddensdT;
double dalphadT = m_waterSS->dthermalExpansionCoeffdT();
doublereal dalphadT = m_waterSS->dthermalExpansionCoeffdT();
double dfac3dT = dgdt * alpha * bfunc * dens / (1.0 - dens)
doublereal dfac3dT = dgdt * alpha * bfunc * dens / (1.0 - dens)
+ gval * dalphadT * bfunc * dens / (1.0 - dens)
+ gval * alpha * bfuncdT * dens / (1.0 - dens)
+ gval * alpha * bfunc * ddensdT / (1.0 - dens)
@ -875,9 +870,9 @@ namespace Cantera {
return dfac1dT + dfac2dT + dfac3dT;
} else if (ifunc == 3) {
double beta = m_waterSS->isothermalCompressibility();
doublereal beta = m_waterSS->isothermalCompressibility();
double dgdp = - bfunc * gval * dens * beta / (1.0 - dens);
doublereal dgdp = - bfunc * gval * dens * beta / (1.0 - dens);
return dgdp;
} else {
@ -887,9 +882,9 @@ namespace Cantera {
}
double PDSS_HKFT::gstar(const double temp, const double pres, const int ifunc) const {
double gval = g(temp, pres, ifunc);
double fval = f(temp, pres, ifunc);
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;
}
@ -914,7 +909,7 @@ namespace Cantera {
*/
struct GeData {
char name[4]; ///< Null Terminated name, First letter capitalized
double GeValue; /// < Gibbs free energies of elements J kmol-1
doublereal GeValue; /// < Gibbs free energies of elements J kmol-1
};
//! Values of G_elements(T=298.15,1atm)
@ -951,12 +946,12 @@ namespace Cantera {
* @exception CanteraError
* If a match is not found, a CanteraError is thrown as well
*/
double PDSS_HKFT::LookupGe(const std::string& s) {
doublereal PDSS_HKFT::LookupGe(const std::string& elemName) {
#ifdef OLDWAY
int num = sizeof(geDataTable) / sizeof(struct GeData);
string s3 = s.substr(0,3);
string s3 = elemName.substr(0,3);
for (int i = 0; i < num; i++) {
//if (!std::strncmp(s.c_str(), aWTable[i].name, 3)) {
//if (!std::strncmp(elemName.c_str(), aWTable[i].name, 3)) {
if (s3 == geDataTable[i].name) {
return (geDataTable[i].GeValue);
}
@ -964,14 +959,14 @@ namespace Cantera {
throw CanteraError("LookupGe", "element " + s + " not found");
return -1.0;
#else
int iE = m_tp->elementIndex(s);
int iE = m_tp->elementIndex(elemName);
if (iE < 0) {
throw CanteraError("PDSS_HKFT::LookupGe", "element " + s + " not found");
throw CanteraError("PDSS_HKFT::LookupGe", "element " + elemName + " not found");
}
doublereal geValue = m_tp->entropyElement298(iE);
if (geValue == ENTROPY298_UNKNOWN) {
throw CanteraError("PDSS_HKFT::LookupGe",
"element " + s + " doesn not have a supplied entropy298");
"element " + elemName + " doesn not have a supplied entropy298");
}
geValue *= (-298.15);
return geValue;
@ -983,11 +978,11 @@ namespace Cantera {
* Ok let's get the element compositions and conversion factors.
*/
int ne = m_tp->nElements();
double na;
double ge;
doublereal na;
doublereal ge;
string ename;
double totalSum = 0.0;
doublereal totalSum = 0.0;
for (int m = 0; m < ne; m++) {
na = m_tp->nAtoms(m_spindex, m);
if (na > 0.0) {
@ -1003,7 +998,7 @@ namespace Cantera {
totalSum -= m_charge_j * ge;
}
// Ok, now do the calculation. Convert to joules kmol-1
double dg = m_deltaG_formation_tr_pr * 4.184 * 1.0E3;
doublereal dg = m_deltaG_formation_tr_pr * 4.184 * 1.0E3;
//! Store the result into an internal variable.
m_Mu0_tr_pr = dg + totalSum;
}

View file

@ -424,7 +424,7 @@ namespace Cantera {
* This is eEqn. 59 in Johnson et al. (1992).
*
*/
double deltaG() const;
doublereal deltaG() const;
//! Main routine that actually calculates the entropy difference
//! between the reference state at Tr, Pr and T,P
@ -432,15 +432,49 @@ namespace Cantera {
* This is eEqn. 61 in Johnson et al. (1992). Actually, there appears to
* be an error in the latter. This is a correction.
*/
double deltaS() const;
doublereal deltaS() const;
//! Internal formula for the calculation of a_g()
/*!
* The output of this is in units of Angstroms
*
* @param temp Temperature (K)
*
* @param ifunc parameters specifying the desired information
* - 0 function value
* - 1 derivative wrt temperature
* - 2 2nd derivative wrt temperature
* - 3 derivative wrt pressure
*/
doublereal ag(const doublereal temp, const int ifunc = 0) const;
double electrostatic_radii_calc();
//! Internal formula for the calculation of b_g()
/*!
* the output of this is unitless
*
* @param temp Temperature (K)
*
* @param ifunc parameters specifying the desired information
* - 0 function value
* - 1 derivative wrt temperature
* - 2 2nd derivative wrt temperature
* - 3 derivative wrt pressure
*/
doublereal bg(const doublereal temp, const int ifunc = 0) const;
double ag(const double temp, const int ifunc = 0) const;
double bg(const double temp, const int ifunc = 0) const;
double g(const double temp, const double pres, const int ifunc = 0) const;
//! function g appearing in the formulation
/*!
* Function g appearing in the Johnson et al formulation
*
* @param temp Temperature kelvin
* @param pres Pressure (pascal)
* @param ifunc parameters specifying the desired information
* - 0 function value
* - 1 derivative wrt temperature
* - 2 2nd derivative wrt temperature
* - 3 derivative wrt pressure
*/
doublereal g(const doublereal temp, const doublereal pres, const int ifunc = 0) const;
//! Difference function f appearing in the formulation
/*!
@ -455,10 +489,44 @@ namespace Cantera {
* - 2 2nd derivative wrt temperature
* - 3 derivative wrt pressure
*/
double f(const double temp, const double pres, const int ifunc = 0) const;
double gstar(const double temp, const double pres, const int ifunc = 0) const;
doublereal f(const doublereal temp, const doublereal pres, const int ifunc = 0) const;
double LookupGe(const std::string& s);
//! Evaluate the Gstar value appearing in the HKFT formulation
/*!
*
* @param temp Temperature kelvin
* @param pres Pressure (pascal)
* @param ifunc parameters specifying the desired information
* - 0 function value
* - 1 derivative wrt temperature
* - 2 2nd derivative wrt temperature
* - 3 derivative wrt pressure
*/
doublereal gstar(const doublereal temp, const doublereal pres, const int ifunc = 0) const;
//! Function to look up Element Free Energies
/*!
*
* This static function looks up the argument string in the
* element database and returns the associated 298 K Gibbs Free energy
* of the element in its stable state
*
* @param elemName String. Only the first 3 characters are significant
*
* @return
* Return value contains the Gibbs free energy for that element
*
* @exception CanteraError
* If a match is not found, a CanteraError is thrown as well
*/
doublereal LookupGe(const std::string& elemName);
//! Translate a Gibbs free energy of formation value to a NIST-based Chemical potential
/*!
* Internally, this function is used to translate the input value, m_deltaG_formation_tr_pr,
* to the internally storred value, m_Mu0_tr_pr.
*/
void convertDGFormation();
private:
@ -475,7 +543,7 @@ namespace Cantera {
/*!
* internal temporary variable
*/
mutable double m_densWaterSS;
mutable doublereal m_densWaterSS;
/**
* Pointer to the water property calculator
@ -545,11 +613,11 @@ namespace Cantera {
//! omega_pr_tr coefficient(cal gmol-1)
doublereal m_omega_pr_tr;
//! y = dZdT = 1/(esp*esp) desp/dT
double m_Y_pr_tr;
//! y = dZdT = 1/(esp*esp) desp/dT at 298.15 and 1 bar
doublereal m_Y_pr_tr;
double m_Z_pr_tr;
//! Z = -1 / relEpsilon at 298.15 and 1 bar
doublereal m_Z_pr_tr;
//! Reference pressure is 1 atm in units of bar= 1.0132
doublereal m_presR_bar;

View file

@ -63,8 +63,8 @@ namespace Cantera {
}
void STITbyPDSS::initAllPtrs(int k, VPSSMgr *vpssmgr_ptr, PDSS *PDSS_ptr) {
AssertThrow(k == m_speciesIndex, "STITbyPDSS::initAllPtrs internal confusion");
void STITbyPDSS::initAllPtrs(int speciesIndex, VPSSMgr *vpssmgr_ptr, PDSS *PDSS_ptr) {
AssertThrow(speciesIndex == m_speciesIndex, "STITbyPDSS::initAllPtrs internal confusion");
m_vpssmgr_ptr = vpssmgr_ptr;
m_PDSS_ptr = PDSS_ptr;
}

View file

@ -151,7 +151,20 @@ namespace Cantera {
};
//! Class for the thermoydnamic manager for an individual species' reference state
//! which usess the PDSS base class to satisfy the requests.
/*!
*
* This class is a pass-through class for handling thermodynamics calls
* for reference state thermo to an pressure dependent standard state (PDSS)
* class. For some situations, it makes no sense to have a reference state
* at all. One example of this is the real water standard state.
*
* What this class does is just to pass through the calls for thermo at (T , p0)
* to the PDSS class, which evaluates the calls at (T, p0).
*
* @ingroup spthermo
*/
class STITbyPDSS : public SpeciesThermoInterpType {
public:
@ -159,11 +172,18 @@ namespace Cantera {
//! Constructor
STITbyPDSS();
//! Constructor
STITbyPDSS(int k, VPSSMgr *vpssmgr_ptr, PDSS *PDSS_ptr);
//! Main Constructor
/*!
*
* @param speciesIndex species index for this object. Note, this must
* agree with what was internally set before.
*
* @param vpssmgr_ptr Pointer to the Variable pressure standard state manager
* that owns the PDSS object that will handle calls for this object
*
* @param PDSS_ptr Pointer to the PDSS object that handles calls for this object
*/
STITbyPDSS(int speciesIndex, VPSSMgr *vpssmgr_ptr, PDSS *PDSS_ptr);
//! copy constructor
/*!
@ -176,8 +196,23 @@ namespace Cantera {
//! duplicator
virtual SpeciesThermoInterpType *duplMyselfAsSpeciesThermoInterpType() const;
void initAllPtrs(int k, VPSSMgr *vpssmgr_ptr, PDSS *PDSS_ptr);
//! Initialize and/or Reinitialize all the pointers for this object
/*!
* This routine is needed because the STITbyPDSS object doesn't own the
* underlying objects. Therefore, shallow copies during duplication operations
* may fail.
*
* @param speciesIndex species index for this object. Note, this must
* agree with what was internally set before.
*
* @param vpssmgr_ptr Pointer to the Variable pressure standard state manager
* that owns the PDSS object that will handle calls for this object
*
* @param PDSS_ptr Pointer to the PDSS object that handles calls for this object
*
*/
void initAllPtrs(int speciesIndex, VPSSMgr *vpssmgr_ptr, PDSS *PDSS_ptr);
//! Returns the minimum temperature that the thermo
//! parameterization is valid
@ -270,8 +305,17 @@ namespace Cantera {
private:
//! Pointer to the Variable pressure standard state manager
//! that owns the PDSS object that will handle calls for this object
VPSSMgr *m_vpssmgr_ptr;
//! Pointer to the PDSS object that handles calls for this object
/*!
* This object is not owned by the current one.
*/
PDSS *m_PDSS_ptr;
//! Species index within the phase
int m_speciesIndex;
};

View file

@ -125,6 +125,7 @@ namespace Cantera {
*/
class UnknownSpeciesThermo : public CanteraError {
public:
//! constructor
/*!
* @param proc name of the procecdure
@ -134,6 +135,11 @@ namespace Cantera {
CanteraError(proc, "Specified species parameterization type (" + int2str(type)
+ ") does not match any known type.") {}
//! Alternate constructor
/*!
* @param proc name of the procecdure
* @param stype String name for the unknown type
*/
UnknownSpeciesThermo(std::string proc, std::string stype) :
CanteraError(proc, "Specified species parameterization type (" + stype
+ ") does not match any known type.") {}