23 int PDSS_HKFT::s_InputInconsistencyErrorExit = 1;
25 PDSS_HKFT::PDSS_HKFT()
28 , m_deltaG_formation_tr_pr(NAN)
29 , m_deltaH_formation_tr_pr(NAN)
42 , m_domega_jdT_prtr(0.0)
60 return deltaH() + enthTRPR;
80 doublereal pbar =
m_pres * 1.0E-5;
81 doublereal c1term =
m_c1;
88 doublereal domega_jdT;
89 doublereal d2omega_jdT2;
95 doublereal nu = 166027;
101 doublereal r_e_j = r_e_j_pr_tr + fabs(
m_charge_j) * gval;
102 doublereal dr_e_jdT = fabs(
m_charge_j) * dgvaldT;
103 doublereal d2r_e_jdT2 = fabs(
m_charge_j) * d2gvaldT2;
104 doublereal r_e_j2 = r_e_j * r_e_j;
107 doublereal r_e_H = 3.082 + gval;
108 doublereal r_e_H2 = r_e_H * r_e_H;
109 omega_j = nu * (charge2 / r_e_j -
m_charge_j / r_e_H);
110 domega_jdT = nu * (-(charge2 / r_e_j2 * dr_e_jdT)
112 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
118 doublereal Y = drelepsilondT / (relepsilon * relepsilon);
121 doublereal X = d2relepsilondT2 / (relepsilon* relepsilon) - 2.0 * relepsilon * Y * Y;
122 doublereal Z = -1.0 / relepsilon;
124 doublereal yterm = 2.0 *
m_temp * Y * domega_jdT;
125 doublereal xterm = omega_j *
m_temp * X;
126 doublereal otterm =
m_temp * d2omega_jdT2 * (Z + 1.0);
129 doublereal Cp_calgmol = c1term + c2term + a3term + a4term + yterm + xterm + otterm + rterm;
132 doublereal Cp = Cp_calgmol *
toSI(
"cal/gmol");
141 doublereal a1term =
m_a1 * 1.0E-5;
142 doublereal a2term =
m_a2 / (2600.E5 +
m_pres);
143 doublereal a3term =
m_a3 * 1.0E-5/ (
m_temp - 228.);
147 doublereal domega_jdP;
152 doublereal nu = 166027.;
159 doublereal r_e_j = r_e_j_pr_tr + fabs(
m_charge_j) * gval;
160 doublereal r_e_H = 3.082 + gval;
162 omega_j = nu * (charge2 / r_e_j -
m_charge_j / r_e_H);
163 doublereal dr_e_jdP = fabs(
m_charge_j) * dgvaldP;
164 domega_jdP = - nu * (charge2 / (r_e_j * r_e_j) * dr_e_jdP)
165 + nu *
m_charge_j / (r_e_H * r_e_H) * dgvaldP;
170 doublereal Q = drelepsilondP / (relepsilon * relepsilon);
171 doublereal Z = -1.0 / relepsilon;
172 doublereal wterm = - domega_jdP * (Z + 1.0);
173 doublereal qterm = - omega_j * Q;
174 doublereal molVol_calgmolPascal = a1term + a2term + a3term + a4term + wterm + qterm;
177 return molVol_calgmolPascal *
toSI(
"cal/gmol");
187 doublereal m_psave =
m_pres;
196 doublereal m_psave =
m_pres;
205 doublereal m_psave =
m_pres;
214 doublereal m_psave =
m_pres;
216 doublereal ee =
cp_R();
223 doublereal m_psave =
m_pres;
272 m_Y_pr_tr = drelepsilondT / (relepsilon * relepsilon);
285 if (fabs(Hcalc -DHjmol) > 100.*
toSI(
"cal/gmol")) {
288 throw CanteraError(
"PDSS_HKFT::initThermo()",
"For {}, DHjmol is" 289 " not consistent with G and S: {} vs {} cal gmol-1",
292 writelog(
"PDSS_HKFT::initThermo() WARNING: DHjmol for {} is not" 293 " consistent with G and S: calculated {} vs input {} cal gmol-1",
295 writelog(
" : continuing with consistent DHjmol = {}", Hcalc/
toSI(
"cal/gmol"));
299 doublereal nu = 166027;
300 doublereal r_e_j_pr_tr;
312 doublereal r_e_j = r_e_j_pr_tr + fabs(
m_charge_j) * gval;
313 doublereal dr_e_jdT = fabs(
m_charge_j) * dgvaldT;
315 + nu *
m_charge_j / (3.082 + gval) / (3.082 + gval) * dgvaldT;
334 m_a3 = a[2] /
toSI(
"cal-K/gmol/bar");
357 "no thermo Node for species " + speciesNode.
name());
361 "thermo model for species isn't hkft: " 362 + speciesNode.
name());
367 "no Thermo::HKFT Node for species " + speciesNode.
name());
372 std::string p0string = hh->
attrib(
"Pref");
373 if (p0string !=
"") {
377 std::string minTstring = hh->
attrib(
"Tmin");
378 if (minTstring !=
"") {
382 std::string maxTstring = hh->
attrib(
"Tmax");
383 if (maxTstring !=
"") {
405 "no standardState Node for species " + speciesNode.
name());
409 "standardState model for species isn't hkft: " 410 + speciesNode.
name());
421 int isum = hasDGO + hasDHO + hasSO;
424 "Missing 2 or more of DG0_f_Pr_Tr, DH0_f_Pr_Tr, or S0_f_Pr_Tr fields. " 425 "Need to supply at least two of these fields");
431 doublereal pbar =
m_pres * 1.0E-5;
435 doublereal c2term = -
m_c2 * (1.0/(
m_temp - 228.) - 1.0/(298.15 - 228.));
438 doublereal a4term =
m_a4 * a3tmp * log((2600. + pbar)/(2600. +
m_presR_bar));
440 doublereal domega_jdT;
445 doublereal nu = 166027;
448 doublereal r_e_j = r_e_j_pr_tr + fabs(
m_charge_j) * gval;
450 doublereal dr_e_jdT = fabs(
m_charge_j) * dgvaldT;
453 + nu *
m_charge_j / (3.082 + gval) / (3.082 + gval) * dgvaldT;
459 doublereal Y = drelepsilondT / (relepsilon * relepsilon);
460 doublereal Z = -1.0 / relepsilon;
462 doublereal yterm =
m_temp * omega_j * Y;
465 doublereal wterm = - omega_j * (Z + 1.0);
468 doublereal otterm =
m_temp * domega_jdT * (Z + 1.0);
471 doublereal deltaH_calgmol = c1term + a1term + a2term + c2term + a3term + a4term
472 + yterm + yrterm + wterm + wrterm + otterm + otrterm;
475 return deltaH_calgmol *
toSI(
"cal/gmol");
480 doublereal pbar =
m_pres * 1.0E-5;
485 doublereal c2term = -
m_c2 * ((1.0/(
m_temp - 228.) - 1.0/(298.15 - 228.)) * (228. -
m_temp)/228.
494 doublereal nu = 166027;
497 doublereal r_e_j = r_e_j_pr_tr + fabs(
m_charge_j) * gval;
502 doublereal Z = -1.0 / relepsilon;
503 doublereal wterm = - omega_j * (Z + 1.0);
506 doublereal deltaG_calgmol = sterm + c1term + a1term + a2term + c2term + a3term + a4term + wterm + wrterm + yterm;
509 return deltaG_calgmol *
toSI(
"cal/gmol");
514 doublereal pbar =
m_pres * 1.0E-5;
516 doublereal c1term =
m_c1 * log(
m_temp/298.15);
517 doublereal c2term = -
m_c2 / 228. * ((1.0/(
m_temp - 228.) - 1.0/(298.15 - 228.))
518 + 1.0 / 228. * log((298.15*(
m_temp-228.)) / (
m_temp*(298.15-228.))));
523 doublereal domega_jdT;
528 doublereal nu = 166027;
532 doublereal r_e_j = r_e_j_pr_tr + fabs(
m_charge_j) * gval;
533 doublereal dr_e_jdT = fabs(
m_charge_j) * dgvaldT;
536 + nu *
m_charge_j / (3.082 + gval) / (3.082 + gval) * dgvaldT;
541 doublereal Y = drelepsilondT / (relepsilon * relepsilon);
542 doublereal Z = -1.0 / relepsilon;
543 doublereal wterm = omega_j * Y;
545 doublereal otterm = domega_jdT * (Z + 1.0);
547 doublereal deltaS_calgmol = c1term + c2term + a3term + a4term + wterm + wrterm + otterm + otrterm;
550 return deltaS_calgmol *
toSI(
"cal/gmol");
555 static doublereal ag_coeff[3] = { -2.037662, 5.747000E-3, -6.557892E-6};
557 return ag_coeff[0] + ag_coeff[1] * temp + ag_coeff[2] * temp * temp;
558 }
else if (ifunc == 1) {
559 return ag_coeff[1] + ag_coeff[2] * 2.0 * temp;
564 return ag_coeff[2] * 2.0;
569 static doublereal bg_coeff[3] = { 6.107361, -1.074377E-2, 1.268348E-5};
571 return bg_coeff[0] + bg_coeff[1] * temp + bg_coeff[2] * temp * temp;
572 }
else if (ifunc == 1) {
573 return bg_coeff[1] + bg_coeff[2] * 2.0 * temp;
578 return bg_coeff[2] * 2.0;
581 doublereal
PDSS_HKFT::f(
const doublereal temp,
const doublereal pres,
const int ifunc)
const 583 static doublereal af_coeff[3] = { 3.666666E1, -0.1504956E-9, 0.5107997E-13};
584 doublereal TC = temp - 273.15;
585 doublereal presBar = pres / 1.0E5;
590 TC = std::min(TC, 355.0);
591 if (presBar > 1000.) {
595 doublereal T1 = (TC-155.0)/300.;
596 doublereal p2 = (1000. - presBar) * (1000. - presBar);
597 doublereal p3 = (1000. - presBar) * p2;
598 doublereal p4 = p2 * p2;
599 doublereal fac2 = af_coeff[1] * p3 + af_coeff[2] * p4;
601 return pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0) * fac2;
602 }
else if (ifunc == 1) {
603 return (4.8 * pow(T1,3.8) + 16.0 * af_coeff[0] * pow(T1, 15.0)) / 300. * fac2;
604 }
else if (ifunc == 2) {
605 return (4.8 * 3.8 * pow(T1,2.8) + 16.0 * 15.0 * af_coeff[0] * pow(T1, 14.0)) / (300. * 300.) * fac2;
606 }
else if (ifunc == 3) {
607 double fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
608 fac2 = - (3.0 * af_coeff[1] * p2 + 4.0 * af_coeff[2] * p3)/ 1.0E5;
615 doublereal
PDSS_HKFT::g(
const doublereal temp,
const doublereal pres,
const int ifunc)
const 617 doublereal afunc =
ag(temp, 0);
618 doublereal bfunc =
bg(temp, 0);
623 doublereal gval = afunc * pow((1.0-dens), bfunc);
629 }
else if (ifunc == 1 || ifunc == 2) {
630 doublereal afuncdT =
ag(temp, 1);
631 doublereal bfuncdT =
bg(temp, 1);
634 doublereal fac1 = afuncdT * gval / afunc;
635 doublereal fac2 = bfuncdT * gval * log(1.0 - dens);
636 doublereal fac3 = gval * alpha * bfunc * dens / (1.0 - dens);
638 doublereal dgdt = fac1 + fac2 + fac3;
643 doublereal afuncdT2 =
ag(temp, 2);
644 doublereal bfuncdT2 =
bg(temp, 2);
645 doublereal dfac1dT = dgdt * afuncdT / afunc + afuncdT2 * gval / afunc
646 - afuncdT * afuncdT * gval / (afunc * afunc);
647 doublereal ddensdT = - alpha * dens;
648 doublereal dfac2dT = bfuncdT2 * gval * log(1.0 - dens)
649 + bfuncdT * dgdt * log(1.0 - dens)
650 - bfuncdT * gval /(1.0 - dens) * ddensdT;
652 doublereal dfac3dT = dgdt * alpha * bfunc * dens / (1.0 - dens)
653 + gval * dalphadT * bfunc * dens / (1.0 - dens)
654 + gval * alpha * bfuncdT * dens / (1.0 - dens)
655 + gval * alpha * bfunc * ddensdT / (1.0 - dens)
656 + gval * alpha * bfunc * dens / ((1.0 - dens) * (1.0 - dens)) * ddensdT;
658 return dfac1dT + dfac2dT + dfac3dT;
659 }
else if (ifunc == 3) {
661 return - bfunc * gval * dens * beta / (1.0 - dens);
668 doublereal
PDSS_HKFT::gstar(
const doublereal temp,
const doublereal pres,
const int ifunc)
const 670 doublereal gval =
g(temp, pres, ifunc);
671 doublereal fval =
f(temp, pres, ifunc);
672 double res = gval - fval;
680 throw CanteraError(
"PDSS_HKFT::LookupGe",
"element " + elemName +
" not found");
685 "element " + elemName +
" does not have a supplied entropy298");
687 return geValue * -298.15;
693 doublereal totalSum = 0.0;
712 doublereal& minTemp_,
713 doublereal& maxTemp_,
714 doublereal& refPressure_)
const doublereal m_a4
Input a4 coefficient (cal K gmol-1)
size_t nElements() const
Number of elements.
doublereal m_densWaterSS
density of standard-state water. internal temporary variable
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
doublereal deltaH() const
Routine that actually calculates the enthalpy difference between the reference state at Tr...
virtual doublereal gibbs_RT() const
Return the molar Gibbs free energy divided by RT.
std::string name() const
Returns the name of the XML node.
virtual doublereal enthalpy_mole() const
Return the molar enthalpy in units of J kmol-1.
virtual void setState_TP(doublereal temp, doublereal pres)
Set the internal temperature and pressure.
void set_c(double *c)
Set "c" coefficients (array of 2 elements).
const XML_Node * findByName(const std::string &nm, int depth=100000) const
This routine carries out a recursive search for an XML node based on the name of the node...
doublereal m_deltaG_formation_tr_pr
Input value of deltaG of Formation at Tr and Pr (cal gmol-1)
const doublereal OneAtm
One atmosphere [Pa].
virtual doublereal gibbs_RT_ref() const
Return the molar Gibbs free energy divided by RT at reference pressure.
virtual void initThermo()
Initialization routine.
doublereal bg(const doublereal temp, const int ifunc=0) const
Internal formula for the calculation of b_g()
doublereal deltaS() const
Main routine that actually calculates the entropy difference between the reference state at Tr...
virtual doublereal entropy_mole() const
Return the molar entropy in units of J kmol-1 K-1.
doublereal toSI(const std::string &unit)
Return the conversion factor to convert unit std::string 'unit' to SI units.
virtual doublereal intEnergy_mole() const
Return the molar internal Energy in units of J kmol-1.
virtual void setPressure(doublereal pres)
Sets the pressure in the object.
const size_t npos
index returned by functions to indicate "no position"
virtual void setParametersFromXML(const XML_Node &speciesNode)
Initialization routine for the PDSS object based on the speciesNode.
size_t elementIndex(const std::string &name) const
Return the index of element named 'name'.
doublereal g(const doublereal temp, const doublereal pres, const int ifunc=0) const
function g appearing in the formulation
doublereal m_Mu0_tr_pr
Value of the Absolute Gibbs Free Energy NIST scale at T_r and P_r.
virtual void reportParams(size_t &kindex, int &type, doublereal *const c, doublereal &minTemp, doublereal &maxTemp, doublereal &refPressure) const
This utility function reports back the type of parameterization and all of the parameters for the spe...
doublereal m_pres
State of the system - pressure.
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
doublereal m_c2
Input c2 coefficient (cal K gmol-1)
Class XML_Node is a tree-based representation of the contents of an XML file.
doublereal gstar(const doublereal temp, const doublereal pres, const int ifunc=0) const
Evaluate the Gstar value appearing in the HKFT formulation.
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
virtual doublereal cp_R_ref() const
Return the molar heat capacity divided by R at reference pressure.
doublereal m_a3
Input a3 coefficient (cal K gmol-1 bar-1)
virtual doublereal cp_R() const
Return the molar const pressure heat capacity divided by RT.
virtual doublereal molarVolume() const
Return the molar volume at standard state.
doublereal m_omega_pr_tr
Input omega_pr_tr coefficient(cal gmol-1)
doublereal m_Y_pr_tr
y = dZdT = 1/(esp*esp) desp/dT at 298.15 and 1 bar
virtual doublereal entropy_R() const
Return the standard state entropy divided by RT.
virtual doublereal isothermalCompressibility() const
Returns the isothermal compressibility. Units: 1/Pa.
static int s_InputInconsistencyErrorExit
Static variable determining error exiting.
virtual void initThermo()
Initialization routine.
doublereal enthalpy_mole2() const
Return the molar enthalpy in units of J kmol-1.
doublereal m_a2
Input a2 coefficient (cal gmol-1)
doublereal m_charge_j
Charge of the ion.
std::unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
virtual void setTemperature(doublereal temp)
Set the internal temperature.
virtual doublereal thermalExpansionCoeff() const
Return the volumetric thermal expansion coefficient. Units: 1/K.
The WaterProps class is used to house several approximation routines for properties of water...
virtual void setState_TP(doublereal temp, doublereal pres)
Set the internal temperature and pressure.
doublereal deltaG() const
Main routine that actually calculates the Gibbs free energy difference between the reference state at...
std::string speciesName(size_t k) const
Name of the species with index k.
Class for the liquid water pressure dependent standard state.
doublereal m_deltaH_formation_tr_pr
Input value of deltaH of Formation at Tr and Pr (cal gmol-1)
doublereal entropyElement298(size_t m) const
Entropy of the element in its standard state at 298 K and 1 bar.
doublereal m_Entrop_tr_pr
Input value of S_j at Tr and Pr (cal gmol-1 K-1)
bool caseInsensitiveEquals(const std::string &input, const std::string &test)
Case insensitive equality predicate.
virtual doublereal density() const
Return the standard state density at standard state.
virtual doublereal entropy_R_ref() const
Return the molar entropy divided by R at reference pressure.
Base class for exceptions thrown by Cantera classes.
doublereal f(const doublereal temp, const doublereal pres, const int ifunc=0) const
Difference function f appearing in the formulation.
void convertDGFormation()
Translate a Gibbs free energy of formation value to a NIST-based Chemical potential.
doublereal pref_safe(doublereal temp) const
Returns a reference pressure value that can be safely calculated by the underlying real equation of s...
doublereal m_presR_bar
Reference pressure is 1 atm in units of bar= 1.0132.
void setDeltaG0(double dg0)
Set Gibbs free energy of formation at Pr, Tr [J/kmol].
virtual doublereal molarVolume_ref() const
Return the molar volume at reference pressure.
virtual doublereal dthermalExpansionCoeffdT() const
Return the derivative of the volumetric thermal expansion coefficient.
doublereal m_maxTemp
Maximum temperature.
doublereal m_domega_jdT_prtr
small value that is not quite zero
doublereal m_minTemp
Minimum temperature.
#define ENTROPY298_UNKNOWN
Number indicating we don't know the entropy of the element in its most stable state at 298...
void set_a(double *a)
Set "a" coefficients (array of 4 elements).
void setOmega(double omega)
Set omega [J/kmol].
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
PDSS_Water * m_waterSS
Water standard state calculator.
doublereal LookupGe(const std::string &elemName)
Function to look up Element Free Energies.
virtual void reportParams(size_t &kindex, int &type, doublereal *const c, doublereal &minTemp, doublereal &maxTemp, doublereal &refPressure) const
This utility function reports back the type of parameterization and all of the parameters for the spe...
Header file for a derived class of ThermoPhase that handles variable pressure standard state methods ...
void setDeltaH0(double dh0)
Set enthalpy of formation at Pr, Tr [J/kmol].
doublereal m_Z_pr_tr
Z = -1 / relEpsilon at 298.15 and 1 bar.
std::string attrib(const std::string &attr) const
Function returns the value of an attribute.
doublereal m_temp
Current temperature used by the PDSS object.
VPStandardStateTP * m_tp
Parent VPStandardStateTP (ThermoPhase) object.
Contains declarations for string manipulation functions within Cantera.
doublereal getFloat(const XML_Node &parent, const std::string &name, const std::string &type)
Get a floating-point value from a child element.
size_t m_spindex
Index of this species within the parent phase.
virtual doublereal gibbs_mole() const
Return the molar Gibbs free energy in units of J kmol-1.
void setParametersFromXML(const XML_Node &speciesNode)
Initialization routine for the PDSS object based on the speciesNode.
virtual doublereal density() const
Return the standard state density at standard state.
doublereal ag(const doublereal temp, const int ifunc=0) const
Internal formula for the calculation of a_g()
doublereal m_c1
Input c1 coefficient (cal gmol-1 K-1)
Namespace for the Cantera kernel.
doublereal nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
doublereal strSItoDbl(const std::string &strSI)
Interpret one or two token string as a single double.
virtual doublereal enthalpy_RT_ref() const
Return the molar enthalpy divided by RT at reference pressure.
virtual doublereal enthalpy_RT() const
Return the standard state molar enthalpy divided by RT.
doublereal fpValueCheck(const std::string &val)
Translate a string into one doublereal value, with error checking.
virtual doublereal cp_mole() const
Return the molar const pressure heat capacity in units of J kmol-1 K-1.
void setS0(double s0)
Set entropy of formation at Pr, Tr [J/kmol/K].
doublereal charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
std::string elementName(size_t m) const
Name of the element with index m.
Declarations for the class PDSS_HKFT (pressure dependent standard state) which handles calculations f...
doublereal m_p0
Reference state pressure of the species.
doublereal m_a1
Input a1 coefficient (cal gmol-1 bar-1)
doublereal m_mw
Molecular Weight of the species.