Cantera  2.4.0
PDSS_HKFT.cpp
Go to the documentation of this file.
1 /**
2  * @file PDSS_HKFT.cpp
3  * Definitions for the class PDSS_HKFT (pressure dependent standard state)
4  * which handles calculations for a single species in a phase using the
5  * HKFT standard state
6  * (see \ref pdssthermo and class \link Cantera::PDSS_HKFT PDSS_HKFT\endlink).
7  */
8 
9 // This file is part of Cantera. See License.txt in the top-level directory or
10 // at http://www.cantera.org/license.txt for license and copyright information.
11 
12 #include "cantera/base/ctml.h"
17 
18 using namespace std;
19 
20 namespace Cantera
21 {
22 // Set the default to error exit if there is an input file inconsistency
23 int PDSS_HKFT::s_InputInconsistencyErrorExit = 1;
24 
25 PDSS_HKFT::PDSS_HKFT()
26  : m_waterSS(0)
27  , m_densWaterSS(-1.0)
28  , m_deltaG_formation_tr_pr(NAN)
29  , m_deltaH_formation_tr_pr(NAN)
30  , m_Mu0_tr_pr(0.0)
31  , m_Entrop_tr_pr(NAN)
32  , m_a1(0.0)
33  , m_a2(0.0)
34  , m_a3(0.0)
35  , m_a4(0.0)
36  , m_c1(0.0)
37  , m_c2(0.0)
38  , m_omega_pr_tr(0.0)
39  , m_Y_pr_tr(0.0)
40  , m_Z_pr_tr(0.0)
41  , m_presR_bar(0.0)
42  , m_domega_jdT_prtr(0.0)
43  , m_charge_j(0.0)
44 {
45  m_pres = OneAtm;
46  m_presR_bar = OneAtm * 1.0E-5;
47  m_presR_bar = 1.0;
48 }
49 
50 doublereal PDSS_HKFT::enthalpy_mole() const
51 {
52  // Ok we may change this evaluation method in the future.
53  doublereal h = gibbs_mole() + m_temp * entropy_mole();
54  return h;
55 }
56 
57 doublereal PDSS_HKFT::enthalpy_mole2() const
58 {
59  double enthTRPR = m_Mu0_tr_pr + 298.15 * m_Entrop_tr_pr * toSI("cal/gmol");
60  return deltaH() + enthTRPR;
61 }
62 
63 doublereal PDSS_HKFT::intEnergy_mole() const
64 {
65  return enthalpy_RT() - molarVolume() * m_pres;
66 }
67 
68 doublereal PDSS_HKFT::entropy_mole() const
69 {
70  return m_Entrop_tr_pr * toSI("cal/gmol") + deltaS();
71 }
72 
73 doublereal PDSS_HKFT::gibbs_mole() const
74 {
75  return m_Mu0_tr_pr + deltaG();
76 }
77 
78 doublereal PDSS_HKFT::cp_mole() const
79 {
80  doublereal pbar = m_pres * 1.0E-5;
81  doublereal c1term = m_c1;
82  doublereal c2term = m_c2 / (m_temp - 228.) / (m_temp - 228.);
83  doublereal a3term = -m_a3 / (m_temp - 228.) / (m_temp - 228.) / (m_temp - 228.) * 2.0 * m_temp * (pbar - m_presR_bar);
84  doublereal a4term = -m_a4 / (m_temp - 228.) / (m_temp - 228.) / (m_temp - 228.) * 2.0 * m_temp
85  * log((2600. + pbar)/(2600. + m_presR_bar));
86 
87  doublereal omega_j;
88  doublereal domega_jdT;
89  doublereal d2omega_jdT2;
90  if (m_charge_j == 0.0) {
91  omega_j = m_omega_pr_tr;
92  domega_jdT = 0.0;
93  d2omega_jdT2 = 0.0;
94  } else {
95  doublereal nu = 166027;
96  doublereal r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
97  doublereal gval = gstar(m_temp, m_pres, 0);
98  doublereal dgvaldT = gstar(m_temp, m_pres, 1);
99  doublereal d2gvaldT2 = gstar(m_temp, m_pres, 2);
100 
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;
105 
106  doublereal charge2 = m_charge_j * m_charge_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)
111  +(m_charge_j / r_e_H2 * dgvaldT));
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
113  -2.0*m_charge_j*dgvaldT*dgvaldT/(r_e_H2*r_e_H) + m_charge_j*d2gvaldT2 /r_e_H2);
114  }
115 
116  doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
117  doublereal drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
118  doublereal Y = drelepsilondT / (relepsilon * relepsilon);
119  doublereal d2relepsilondT2 = m_waterProps->relEpsilon(m_temp, m_pres, 2);
120 
121  doublereal X = d2relepsilondT2 / (relepsilon* relepsilon) - 2.0 * relepsilon * Y * Y;
122  doublereal Z = -1.0 / relepsilon;
123 
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);
127  doublereal rterm = - m_domega_jdT_prtr * (m_Z_pr_tr + 1.0);
128 
129  doublereal Cp_calgmol = c1term + c2term + a3term + a4term + yterm + xterm + otterm + rterm;
130 
131  // Convert to Joules / kmol
132  doublereal Cp = Cp_calgmol * toSI("cal/gmol");
133 
134  return Cp;
135 }
136 
137 doublereal PDSS_HKFT::molarVolume() const
138 {
139  // Initially do all calculations in (cal/gmol/Pa)
140 
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.);
144  doublereal a4term = m_a4 / (m_temp - 228.) / (2600.E5 + m_pres);
145 
146  doublereal omega_j;
147  doublereal domega_jdP;
148  if (m_charge_j == 0.0) {
149  omega_j = m_omega_pr_tr;
150  domega_jdP = 0.0;
151  } else {
152  doublereal nu = 166027.;
153  doublereal charge2 = m_charge_j * m_charge_j;
154  doublereal r_e_j_pr_tr = charge2 / (m_omega_pr_tr/nu + m_charge_j/3.082);
155 
156  doublereal gval = gstar(m_temp, m_pres, 0);
157  doublereal dgvaldP = gstar(m_temp, m_pres, 3);
158 
159  doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
160  doublereal r_e_H = 3.082 + gval;
161 
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;
166  }
167 
168  doublereal drelepsilondP = m_waterProps->relEpsilon(m_temp, m_pres, 3);
169  doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
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;
175 
176  // Convert to m**3 / kmol from (cal/gmol/Pa)
177  return molVol_calgmolPascal * toSI("cal/gmol");
178 }
179 
180 doublereal PDSS_HKFT::density() const
181 {
182  return m_mw / molarVolume();
183 }
184 
185 doublereal PDSS_HKFT::gibbs_RT_ref() const
186 {
187  doublereal m_psave = m_pres;
189  doublereal ee = gibbs_RT();
190  m_pres = m_psave;
191  return ee;
192 }
193 
194 doublereal PDSS_HKFT::enthalpy_RT_ref() const
195 {
196  doublereal m_psave = m_pres;
198  doublereal hh = enthalpy_RT();
199  m_pres = m_psave;
200  return hh;
201 }
202 
203 doublereal PDSS_HKFT::entropy_R_ref() const
204 {
205  doublereal m_psave = m_pres;
207  doublereal ee = entropy_R();
208  m_pres = m_psave;
209  return ee;
210 }
211 
212 doublereal PDSS_HKFT::cp_R_ref() const
213 {
214  doublereal m_psave = m_pres;
216  doublereal ee = cp_R();
217  m_pres = m_psave;
218  return ee;
219 }
220 
221 doublereal PDSS_HKFT::molarVolume_ref() const
222 {
223  doublereal m_psave = m_pres;
225  doublereal ee = molarVolume();
226  m_pres = m_psave;
227  return ee;
228 }
229 
230 void PDSS_HKFT::setState_TP(doublereal temp, doublereal pres)
231 {
232  setTemperature(temp);
233  setPressure(pres);
234 }
235 
237 {
239 
240  // Ok, if we are missing one, then we construct its value from the other two.
241  // This code has been internally verified.
243  if (std::isnan(m_deltaH_formation_tr_pr)) {
245  doublereal Hcalc = m_Mu0_tr_pr + 298.15 * (m_Entrop_tr_pr * toSI("cal/gmol"));
246  m_deltaH_formation_tr_pr = Hcalc / toSI("cal/gmol");
247  } else if (std::isnan(m_deltaG_formation_tr_pr)) {
248  doublereal DHjmol = m_deltaH_formation_tr_pr * toSI("cal/gmol");
249  m_Mu0_tr_pr = DHjmol - 298.15 * (m_Entrop_tr_pr * toSI("cal/gmol"));
250  m_deltaG_formation_tr_pr = m_Mu0_tr_pr / toSI("cal/gmol");
251  double tmp = m_Mu0_tr_pr;
253  double totalSum = m_Mu0_tr_pr - tmp;
254  m_Mu0_tr_pr = tmp;
255  m_deltaG_formation_tr_pr = (m_Mu0_tr_pr - totalSum)/ toSI("cal/gmol");
256  } else if (std::isnan(m_Entrop_tr_pr)) {
258  doublereal DHjmol = m_deltaH_formation_tr_pr * toSI("cal/gmol");
259  m_Entrop_tr_pr = (DHjmol - m_Mu0_tr_pr) / (298.15 * toSI("cal/gmol"));
260  }
261 
262  m_waterSS = &dynamic_cast<PDSS_Water&>(*m_tp->providePDSS(0));
263 
264  // Section to initialize m_Z_pr_tr and m_Y_pr_tr
265  m_temp = 273.15 + 25.;
266  m_pres = OneAtm;
267  doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
270  m_Z_pr_tr = -1.0 / relepsilon;
271  doublereal drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
272  m_Y_pr_tr = drelepsilondT / (relepsilon * relepsilon);
273  m_waterProps.reset(new WaterProps(m_waterSS));
274  m_presR_bar = OneAtm / 1.0E5;
275  m_presR_bar = 1.0;
277 
278  // Ok, we have mu. Let's check it against the input value
279  // of DH_F to see that we have some internal consistency
280  doublereal Hcalc = m_Mu0_tr_pr + 298.15 * (m_Entrop_tr_pr * toSI("cal/gmol"));
281  doublereal DHjmol = m_deltaH_formation_tr_pr * toSI("cal/gmol");
282 
283  // If the discrepancy is greater than 100 cal gmol-1, print
284  // an error and exit.
285  if (fabs(Hcalc -DHjmol) > 100.* toSI("cal/gmol")) {
286  std::string sname = m_tp->speciesName(m_spindex);
288  throw CanteraError("PDSS_HKFT::initThermo()", "For {}, DHjmol is"
289  " not consistent with G and S: {} vs {} cal gmol-1",
290  sname, Hcalc/toSI("cal/gmol"), m_deltaH_formation_tr_pr);
291  } else {
292  writelog("PDSS_HKFT::initThermo() WARNING: DHjmol for {} is not"
293  " consistent with G and S: calculated {} vs input {} cal gmol-1",
294  sname, Hcalc/toSI("cal/gmol"), m_deltaH_formation_tr_pr);
295  writelog(" : continuing with consistent DHjmol = {}", Hcalc/toSI("cal/gmol"));
296  m_deltaH_formation_tr_pr = Hcalc / toSI("cal/gmol");
297  }
298  }
299  doublereal nu = 166027;
300  doublereal r_e_j_pr_tr;
301  if (m_charge_j != 0.0) {
302  r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
303  } else {
304  r_e_j_pr_tr = 0.0;
305  }
306 
307  if (m_charge_j == 0.0) {
308  m_domega_jdT_prtr = 0.0;
309  } else {
310  doublereal gval = gstar(m_temp, m_pres, 0);
311  doublereal dgvaldT = gstar(m_temp, m_pres, 1);
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;
314  m_domega_jdT_prtr = - nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdT)
315  + nu * m_charge_j / (3.082 + gval) / (3.082 + gval) * dgvaldT;
316  }
317 }
318 
319 void PDSS_HKFT::setDeltaH0(double dh0) {
320  m_deltaH_formation_tr_pr = dh0 / toSI("cal/gmol");
321 }
322 
323 void PDSS_HKFT::setDeltaG0(double dg0) {
324  m_deltaG_formation_tr_pr = dg0 / toSI("cal/gmol");
325 }
326 
327 void PDSS_HKFT::setS0(double s0) {
328  m_Entrop_tr_pr = s0 / toSI("cal/gmol/K");
329 }
330 
331 void PDSS_HKFT::set_a(double* a) {
332  m_a1 = a[0] / toSI("cal/gmol/bar");
333  m_a2 = a[1] / toSI("cal/gmol");
334  m_a3 = a[2] / toSI("cal-K/gmol/bar");
335  m_a4 = a[3] / toSI("cal-K/gmol");
336 }
337 
338 void PDSS_HKFT::set_c(double* c) {
339  m_c1 = c[0] / toSI("cal/gmol/K");
340  m_c2 = c[1] / toSI("cal-K/gmol");
341 }
342 
343 void PDSS_HKFT::setOmega(double omega) {
344  m_omega_pr_tr = omega / toSI("cal/gmol");
345 }
346 
348 {
349  PDSS::setParametersFromXML(speciesNode);
350  int hasDGO = 0;
351  int hasSO = 0;
352  int hasDHO = 0;
353 
354  const XML_Node* tn = speciesNode.findByName("thermo");
355  if (!tn) {
356  throw CanteraError("PDSS_HKFT::constructPDSSXML",
357  "no thermo Node for species " + speciesNode.name());
358  }
359  if (!caseInsensitiveEquals(tn->attrib("model"), "hkft")) {
360  throw CanteraError("PDSS_HKFT::initThermoXML",
361  "thermo model for species isn't hkft: "
362  + speciesNode.name());
363  }
364  const XML_Node* hh = tn->findByName("HKFT");
365  if (!hh) {
366  throw CanteraError("PDSS_HKFT::constructPDSSXML",
367  "no Thermo::HKFT Node for species " + speciesNode.name());
368  }
369 
370  // go get the attributes
371  m_p0 = OneAtm;
372  std::string p0string = hh->attrib("Pref");
373  if (p0string != "") {
374  m_p0 = strSItoDbl(p0string);
375  }
376 
377  std::string minTstring = hh->attrib("Tmin");
378  if (minTstring != "") {
379  m_minTemp = fpValueCheck(minTstring);
380  }
381 
382  std::string maxTstring = hh->attrib("Tmax");
383  if (maxTstring != "") {
384  m_maxTemp = fpValueCheck(maxTstring);
385  }
386 
387  if (hh->hasChild("DG0_f_Pr_Tr")) {
388  setDeltaG0(getFloat(*hh, "DG0_f_Pr_Tr", "toSI"));
389  hasDGO = 1;
390  }
391 
392  if (hh->hasChild("DH0_f_Pr_Tr")) {
393  setDeltaH0(getFloat(*hh, "DH0_f_Pr_Tr", "toSI"));
394  hasDHO = 1;
395  }
396 
397  if (hh->hasChild("S0_Pr_Tr")) {
398  setS0(getFloat(*hh, "S0_Pr_Tr", "toSI"));
399  hasSO = 1;
400  }
401 
402  const XML_Node* ss = speciesNode.findByName("standardState");
403  if (!ss) {
404  throw CanteraError("PDSS_HKFT::constructPDSSXML",
405  "no standardState Node for species " + speciesNode.name());
406  }
407  if (!caseInsensitiveEquals(ss->attrib("model"), "hkft")) {
408  throw CanteraError("PDSS_HKFT::initThermoXML",
409  "standardState model for species isn't hkft: "
410  + speciesNode.name());
411  }
412  double a[4] = {getFloat(*ss, "a1", "toSI"), getFloat(*ss, "a2", "toSI"),
413  getFloat(*ss, "a3", "toSI"), getFloat(*ss, "a4", "toSI")};
414  set_a(a);
415 
416  double c[2] = {getFloat(*ss, "c1", "toSI"), getFloat(*ss, "c2", "toSI")};
417  set_c(c);
418 
419  setOmega(getFloat(*ss, "omega_Pr_Tr", "toSI"));
420 
421  int isum = hasDGO + hasDHO + hasSO;
422  if (isum < 2) {
423  throw CanteraError("PDSS_HKFT::constructPDSSXML",
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");
426  }
427 }
428 
429 doublereal PDSS_HKFT::deltaH() const
430 {
431  doublereal pbar = m_pres * 1.0E-5;
432  doublereal c1term = m_c1 * (m_temp - 298.15);
433  doublereal a1term = m_a1 * (pbar - m_presR_bar);
434  doublereal a2term = m_a2 * log((2600. + pbar)/(2600. + m_presR_bar));
435  doublereal c2term = -m_c2 * (1.0/(m_temp - 228.) - 1.0/(298.15 - 228.));
436  double a3tmp = (2.0 * m_temp - 228.)/ (m_temp - 228.) /(m_temp - 228.);
437  doublereal a3term = m_a3 * a3tmp * (pbar - m_presR_bar);
438  doublereal a4term = m_a4 * a3tmp * log((2600. + pbar)/(2600. + m_presR_bar));
439  doublereal omega_j;
440  doublereal domega_jdT;
441  if (m_charge_j == 0.0) {
442  omega_j = m_omega_pr_tr;
443  domega_jdT = 0.0;
444  } else {
445  doublereal nu = 166027;
446  doublereal r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
447  doublereal gval = gstar(m_temp, m_pres, 0);
448  doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
449  doublereal dgvaldT = gstar(m_temp, m_pres, 1);
450  doublereal dr_e_jdT = fabs(m_charge_j) * dgvaldT;
451  omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval));
452  domega_jdT = - nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdT)
453  + nu * m_charge_j / (3.082 + gval) / (3.082 + gval) * dgvaldT;
454  }
455 
456  doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
457  doublereal drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
458 
459  doublereal Y = drelepsilondT / (relepsilon * relepsilon);
460  doublereal Z = -1.0 / relepsilon;
461 
462  doublereal yterm = m_temp * omega_j * Y;
463  doublereal yrterm = - 298.15 * m_omega_pr_tr * m_Y_pr_tr;
464 
465  doublereal wterm = - omega_j * (Z + 1.0);
466  doublereal wrterm = + m_omega_pr_tr * (m_Z_pr_tr + 1.0);
467 
468  doublereal otterm = m_temp * domega_jdT * (Z + 1.0);
469  doublereal otrterm = - m_temp * m_domega_jdT_prtr * (m_Z_pr_tr + 1.0);
470 
471  doublereal deltaH_calgmol = c1term + a1term + a2term + c2term + a3term + a4term
472  + yterm + yrterm + wterm + wrterm + otterm + otrterm;
473 
474  // Convert to Joules / kmol
475  return deltaH_calgmol * toSI("cal/gmol");
476 }
477 
478 doublereal PDSS_HKFT::deltaG() const
479 {
480  doublereal pbar = m_pres * 1.0E-5;
481  doublereal sterm = - m_Entrop_tr_pr * (m_temp - 298.15);
482  doublereal c1term = -m_c1 * (m_temp * log(m_temp/298.15) - (m_temp - 298.15));
483  doublereal a1term = m_a1 * (pbar - m_presR_bar);
484  doublereal a2term = m_a2 * log((2600. + pbar)/(2600. + m_presR_bar));
485  doublereal c2term = -m_c2 * ((1.0/(m_temp - 228.) - 1.0/(298.15 - 228.)) * (228. - m_temp)/228.
486  - m_temp / (228.*228.) * log((298.15*(m_temp-228.)) / (m_temp*(298.15-228.))));
487  doublereal a3term = m_a3 / (m_temp - 228.) * (pbar - m_presR_bar);
488  doublereal a4term = m_a4 / (m_temp - 228.) * log((2600. + pbar)/(2600. + m_presR_bar));
489 
490  doublereal omega_j;
491  if (m_charge_j == 0.0) {
492  omega_j = m_omega_pr_tr;
493  } else {
494  doublereal nu = 166027;
495  doublereal r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
496  doublereal gval = gstar(m_temp, m_pres, 0);
497  doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
498  omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval));
499  }
500 
501  doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
502  doublereal Z = -1.0 / relepsilon;
503  doublereal wterm = - omega_j * (Z + 1.0);
504  doublereal wrterm = m_omega_pr_tr * (m_Z_pr_tr + 1.0);
505  doublereal yterm = m_omega_pr_tr * m_Y_pr_tr * (m_temp - 298.15);
506  doublereal deltaG_calgmol = sterm + c1term + a1term + a2term + c2term + a3term + a4term + wterm + wrterm + yterm;
507 
508  // Convert to Joules / kmol
509  return deltaG_calgmol * toSI("cal/gmol");
510 }
511 
512 doublereal PDSS_HKFT::deltaS() const
513 {
514  doublereal pbar = m_pres * 1.0E-5;
515 
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.))));
519  doublereal a3term = m_a3 / (m_temp - 228.) / (m_temp - 228.) * (pbar - m_presR_bar);
520  doublereal a4term = m_a4 / (m_temp - 228.) / (m_temp - 228.) * log((2600. + pbar)/(2600. + m_presR_bar));
521 
522  doublereal omega_j;
523  doublereal domega_jdT;
524  if (m_charge_j == 0.0) {
525  omega_j = m_omega_pr_tr;
526  domega_jdT = 0.0;
527  } else {
528  doublereal nu = 166027;
529  doublereal r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
530  doublereal gval = gstar(m_temp, m_pres, 0);
531  doublereal dgvaldT = gstar(m_temp, m_pres, 1);
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;
534  omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval));
535  domega_jdT = - nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdT)
536  + nu * m_charge_j / (3.082 + gval) / (3.082 + gval) * dgvaldT;
537  }
538 
539  doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
540  doublereal drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
541  doublereal Y = drelepsilondT / (relepsilon * relepsilon);
542  doublereal Z = -1.0 / relepsilon;
543  doublereal wterm = omega_j * Y;
544  doublereal wrterm = - m_omega_pr_tr * m_Y_pr_tr;
545  doublereal otterm = domega_jdT * (Z + 1.0);
546  doublereal otrterm = - m_domega_jdT_prtr * (m_Z_pr_tr + 1.0);
547  doublereal deltaS_calgmol = c1term + c2term + a3term + a4term + wterm + wrterm + otterm + otrterm;
548 
549  // Convert to Joules / kmol
550  return deltaS_calgmol * toSI("cal/gmol");
551 }
552 
553 doublereal PDSS_HKFT::ag(const doublereal temp, const int ifunc) const
554 {
555  static doublereal ag_coeff[3] = { -2.037662, 5.747000E-3, -6.557892E-6};
556  if (ifunc == 0) {
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;
560  }
561  if (ifunc != 2) {
562  return 0.0;
563  }
564  return ag_coeff[2] * 2.0;
565 }
566 
567 doublereal PDSS_HKFT::bg(const doublereal temp, const int ifunc) const
568 {
569  static doublereal bg_coeff[3] = { 6.107361, -1.074377E-2, 1.268348E-5};
570  if (ifunc == 0) {
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;
574  }
575  if (ifunc != 2) {
576  return 0.0;
577  }
578  return bg_coeff[2] * 2.0;
579 }
580 
581 doublereal PDSS_HKFT::f(const doublereal temp, const doublereal pres, const int ifunc) const
582 {
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;
586 
587  if (TC < 155.0) {
588  return 0.0;
589  }
590  TC = std::min(TC, 355.0);
591  if (presBar > 1000.) {
592  return 0.0;
593  }
594 
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;
600  if (ifunc == 0) {
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;
609  return fac1 * fac2;
610  } else {
611  throw CanteraError("HKFT_PDSS::gg", "unimplemented");
612  }
613 }
614 
615 doublereal PDSS_HKFT::g(const doublereal temp, const doublereal pres, const int ifunc) const
616 {
617  doublereal afunc = ag(temp, 0);
618  doublereal bfunc = bg(temp, 0);
619  m_waterSS->setState_TP(temp, pres);
621  // density in gm cm-3
622  doublereal dens = m_densWaterSS * 1.0E-3;
623  doublereal gval = afunc * pow((1.0-dens), bfunc);
624  if (dens >= 1.0) {
625  return 0.0;
626  }
627  if (ifunc == 0) {
628  return gval;
629  } else if (ifunc == 1 || ifunc == 2) {
630  doublereal afuncdT = ag(temp, 1);
631  doublereal bfuncdT = bg(temp, 1);
632  doublereal alpha = m_waterSS->thermalExpansionCoeff();
633 
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);
637 
638  doublereal dgdt = fac1 + fac2 + fac3;
639  if (ifunc == 1) {
640  return dgdt;
641  }
642 
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;
651  doublereal dalphadT = m_waterSS->dthermalExpansionCoeffdT();
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;
657 
658  return dfac1dT + dfac2dT + dfac3dT;
659  } else if (ifunc == 3) {
660  doublereal beta = m_waterSS->isothermalCompressibility();
661  return - bfunc * gval * dens * beta / (1.0 - dens);
662  } else {
663  throw CanteraError("HKFT_PDSS::g", "unimplemented");
664  }
665  return 0.0;
666 }
667 
668 doublereal PDSS_HKFT::gstar(const doublereal temp, const doublereal pres, const int ifunc) const
669 {
670  doublereal gval = g(temp, pres, ifunc);
671  doublereal fval = f(temp, pres, ifunc);
672  double res = gval - fval;
673  return res;
674 }
675 
676 doublereal PDSS_HKFT::LookupGe(const std::string& elemName)
677 {
678  size_t iE = m_tp->elementIndex(elemName);
679  if (iE == npos) {
680  throw CanteraError("PDSS_HKFT::LookupGe", "element " + elemName + " not found");
681  }
682  doublereal geValue = m_tp->entropyElement298(iE);
683  if (geValue == ENTROPY298_UNKNOWN) {
684  throw CanteraError("PDSS_HKFT::LookupGe",
685  "element " + elemName + " does not have a supplied entropy298");
686  }
687  return geValue * -298.15;
688 }
689 
691 {
692  // Ok let's get the element compositions and conversion factors.
693  doublereal totalSum = 0.0;
694  for (size_t m = 0; m < m_tp->nElements(); m++) {
695  double na = m_tp->nAtoms(m_spindex, m);
696  if (na > 0.0) {
697  totalSum += na * LookupGe(m_tp->elementName(m));
698  }
699  }
700  // Add in the charge
701  if (m_charge_j != 0.0) {
702  totalSum -= m_charge_j * LookupGe("H");
703  }
704  // Ok, now do the calculation. Convert to joules kmol-1
705  doublereal dg = m_deltaG_formation_tr_pr * toSI("cal/gmol");
706  //! Store the result into an internal variable.
707  m_Mu0_tr_pr = dg + totalSum;
708 }
709 
710 void PDSS_HKFT::reportParams(size_t& kindex, int& type,
711  doublereal* const c,
712  doublereal& minTemp_,
713  doublereal& maxTemp_,
714  doublereal& refPressure_) const
715 {
716  // Fill in the first part
717  PDSS::reportParams(kindex, type, c, minTemp_, maxTemp_,
718  refPressure_);
719 
722  c[2] = m_Mu0_tr_pr;
723  c[3] = m_Entrop_tr_pr;
724  c[4] = m_a1;
725  c[5] = m_a2;
726  c[6] = m_a3;
727  c[7] = m_a4;
728  c[8] = m_c1;
729  c[9] = m_c2;
730  c[10] = m_omega_pr_tr;
731 }
732 
733 }
doublereal m_a4
Input a4 coefficient (cal K gmol-1)
Definition: PDSS_HKFT.h:313
size_t nElements() const
Number of elements.
Definition: Phase.cpp:88
doublereal m_densWaterSS
density of standard-state water. internal temporary variable
Definition: PDSS_HKFT.h:265
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...
Definition: PDSS_HKFT.cpp:429
virtual doublereal gibbs_RT() const
Return the molar Gibbs free energy divided by RT.
Definition: PDSS.cpp:221
std::string name() const
Returns the name of the XML node.
Definition: xml.h:370
virtual doublereal enthalpy_mole() const
Return the molar enthalpy in units of J kmol-1.
Definition: PDSS_HKFT.cpp:50
virtual void setState_TP(doublereal temp, doublereal pres)
Set the internal temperature and pressure.
Definition: PDSS_HKFT.cpp:230
void set_c(double *c)
Set "c" coefficients (array of 2 elements).
Definition: PDSS_HKFT.cpp:338
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...
Definition: xml.cpp:695
doublereal m_deltaG_formation_tr_pr
Input value of deltaG of Formation at Tr and Pr (cal gmol-1)
Definition: PDSS_HKFT.h:277
const doublereal OneAtm
One atmosphere [Pa].
Definition: ct_defs.h:69
virtual doublereal gibbs_RT_ref() const
Return the molar Gibbs free energy divided by RT at reference pressure.
Definition: PDSS_HKFT.cpp:185
virtual void initThermo()
Initialization routine.
Definition: PDSS.h:452
doublereal bg(const doublereal temp, const int ifunc=0) const
Internal formula for the calculation of b_g()
Definition: PDSS_HKFT.cpp:567
doublereal deltaS() const
Main routine that actually calculates the entropy difference between the reference state at Tr...
Definition: PDSS_HKFT.cpp:512
virtual doublereal entropy_mole() const
Return the molar entropy in units of J kmol-1 K-1.
Definition: PDSS_HKFT.cpp:68
doublereal toSI(const std::string &unit)
Return the conversion factor to convert unit std::string &#39;unit&#39; to SI units.
Definition: global.cpp:135
virtual doublereal intEnergy_mole() const
Return the molar internal Energy in units of J kmol-1.
Definition: PDSS_HKFT.cpp:63
virtual void setPressure(doublereal pres)
Sets the pressure in the object.
Definition: PDSS.cpp:157
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
virtual void setParametersFromXML(const XML_Node &speciesNode)
Initialization routine for the PDSS object based on the speciesNode.
Definition: PDSS.h:459
size_t elementIndex(const std::string &name) const
Return the index of element named &#39;name&#39;.
Definition: Phase.cpp:113
doublereal g(const doublereal temp, const doublereal pres, const int ifunc=0) const
function g appearing in the formulation
Definition: PDSS_HKFT.cpp:615
doublereal m_Mu0_tr_pr
Value of the Absolute Gibbs Free Energy NIST scale at T_r and P_r.
Definition: PDSS_HKFT.h:295
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...
Definition: PDSS_HKFT.cpp:710
doublereal m_pres
State of the system - pressure.
Definition: PDSS.h:483
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:153
doublereal m_c2
Input c2 coefficient (cal K gmol-1)
Definition: PDSS_HKFT.h:319
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:97
doublereal gstar(const doublereal temp, const doublereal pres, const int ifunc=0) const
Evaluate the Gstar value appearing in the HKFT formulation.
Definition: PDSS_HKFT.cpp:668
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
STL namespace.
virtual doublereal cp_R_ref() const
Return the molar heat capacity divided by R at reference pressure.
Definition: PDSS_HKFT.cpp:212
doublereal m_a3
Input a3 coefficient (cal K gmol-1 bar-1)
Definition: PDSS_HKFT.h:310
virtual doublereal cp_R() const
Return the molar const pressure heat capacity divided by RT.
Definition: PDSS.cpp:226
virtual doublereal molarVolume() const
Return the molar volume at standard state.
Definition: PDSS_HKFT.cpp:137
doublereal m_omega_pr_tr
Input omega_pr_tr coefficient(cal gmol-1)
Definition: PDSS_HKFT.h:322
doublereal m_Y_pr_tr
y = dZdT = 1/(esp*esp) desp/dT at 298.15 and 1 bar
Definition: PDSS_HKFT.h:325
virtual doublereal entropy_R() const
Return the standard state entropy divided by RT.
Definition: PDSS.cpp:216
virtual doublereal isothermalCompressibility() const
Returns the isothermal compressibility. Units: 1/Pa.
Definition: PDSS_Water.cpp:191
static int s_InputInconsistencyErrorExit
Static variable determining error exiting.
Definition: PDSS_HKFT.h:344
virtual void initThermo()
Initialization routine.
Definition: PDSS_HKFT.cpp:236
doublereal enthalpy_mole2() const
Return the molar enthalpy in units of J kmol-1.
Definition: PDSS_HKFT.cpp:57
doublereal m_a2
Input a2 coefficient (cal gmol-1)
Definition: PDSS_HKFT.h:307
doublereal m_charge_j
Charge of the ion.
Definition: PDSS_HKFT.h:337
std::unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
Definition: PDSS_HKFT.h:268
virtual void setTemperature(doublereal temp)
Set the internal temperature.
Definition: PDSS.cpp:167
virtual doublereal thermalExpansionCoeff() const
Return the volumetric thermal expansion coefficient. Units: 1/K.
Definition: PDSS_Water.cpp:170
The WaterProps class is used to house several approximation routines for properties of water...
Definition: WaterProps.h:94
virtual void setState_TP(doublereal temp, doublereal pres)
Set the internal temperature and pressure.
Definition: PDSS_Water.cpp:228
doublereal deltaG() const
Main routine that actually calculates the Gibbs free energy difference between the reference state at...
Definition: PDSS_HKFT.cpp:478
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:191
Class for the liquid water pressure dependent standard state.
Definition: PDSS_Water.h:49
doublereal m_deltaH_formation_tr_pr
Input value of deltaH of Formation at Tr and Pr (cal gmol-1)
Definition: PDSS_HKFT.h:286
doublereal entropyElement298(size_t m) const
Entropy of the element in its standard state at 298 K and 1 bar.
Definition: Phase.cpp:133
doublereal m_Entrop_tr_pr
Input value of S_j at Tr and Pr (cal gmol-1 K-1)
Definition: PDSS_HKFT.h:301
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.
Definition: PDSS_HKFT.cpp:180
virtual doublereal entropy_R_ref() const
Return the molar entropy divided by R at reference pressure.
Definition: PDSS_HKFT.cpp:203
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
doublereal f(const doublereal temp, const doublereal pres, const int ifunc=0) const
Difference function f appearing in the formulation.
Definition: PDSS_HKFT.cpp:581
void convertDGFormation()
Translate a Gibbs free energy of formation value to a NIST-based Chemical potential.
Definition: PDSS_HKFT.cpp:690
doublereal pref_safe(doublereal temp) const
Returns a reference pressure value that can be safely calculated by the underlying real equation of s...
Definition: PDSS_Water.cpp:241
doublereal m_presR_bar
Reference pressure is 1 atm in units of bar= 1.0132.
Definition: PDSS_HKFT.h:331
void setDeltaG0(double dg0)
Set Gibbs free energy of formation at Pr, Tr [J/kmol].
Definition: PDSS_HKFT.cpp:323
virtual doublereal molarVolume_ref() const
Return the molar volume at reference pressure.
Definition: PDSS_HKFT.cpp:221
virtual doublereal dthermalExpansionCoeffdT() const
Return the derivative of the volumetric thermal expansion coefficient.
Definition: PDSS_Water.cpp:175
doublereal m_maxTemp
Maximum temperature.
Definition: PDSS.h:492
doublereal m_domega_jdT_prtr
small value that is not quite zero
Definition: PDSS_HKFT.h:334
doublereal m_minTemp
Minimum temperature.
Definition: PDSS.h:489
#define ENTROPY298_UNKNOWN
Number indicating we don&#39;t know the entropy of the element in its most stable state at 298...
Definition: Elements.h:87
void set_a(double *a)
Set "a" coefficients (array of 4 elements).
Definition: PDSS_HKFT.cpp:331
void setOmega(double omega)
Set omega [J/kmol].
Definition: PDSS_HKFT.cpp:343
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
Definition: xml.cpp:536
PDSS_Water * m_waterSS
Water standard state calculator.
Definition: PDSS_HKFT.h:262
doublereal LookupGe(const std::string &elemName)
Function to look up Element Free Energies.
Definition: PDSS_HKFT.cpp:676
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...
Definition: PDSS.cpp:196
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].
Definition: PDSS_HKFT.cpp:319
doublereal m_Z_pr_tr
Z = -1 / relEpsilon at 298.15 and 1 bar.
Definition: PDSS_HKFT.h:328
std::string attrib(const std::string &attr) const
Function returns the value of an attribute.
Definition: xml.cpp:500
doublereal m_temp
Current temperature used by the PDSS object.
Definition: PDSS.h:480
VPStandardStateTP * m_tp
Parent VPStandardStateTP (ThermoPhase) object.
Definition: PDSS_HKFT.h:142
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.
Definition: ctml.cpp:164
size_t m_spindex
Index of this species within the parent phase.
Definition: PDSS_HKFT.h:143
virtual doublereal gibbs_mole() const
Return the molar Gibbs free energy in units of J kmol-1.
Definition: PDSS_HKFT.cpp:73
void setParametersFromXML(const XML_Node &speciesNode)
Initialization routine for the PDSS object based on the speciesNode.
Definition: PDSS_HKFT.cpp:347
virtual doublereal density() const
Return the standard state density at standard state.
Definition: PDSS_Water.cpp:217
doublereal ag(const doublereal temp, const int ifunc=0) const
Internal formula for the calculation of a_g()
Definition: PDSS_HKFT.cpp:553
doublereal m_c1
Input c1 coefficient (cal gmol-1 K-1)
Definition: PDSS_HKFT.h:316
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
doublereal nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
Definition: Phase.cpp:161
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.
Definition: PDSS_HKFT.cpp:194
virtual doublereal enthalpy_RT() const
Return the standard state molar enthalpy divided by RT.
Definition: PDSS.cpp:211
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.
Definition: PDSS_HKFT.cpp:78
void setS0(double s0)
Set entropy of formation at Pr, Tr [J/kmol/K].
Definition: PDSS_HKFT.cpp:327
doublereal charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition: Phase.h:577
std::string elementName(size_t m) const
Name of the element with index m.
Definition: Phase.cpp:107
Declarations for the class PDSS_HKFT (pressure dependent standard state) which handles calculations f...
doublereal m_p0
Reference state pressure of the species.
Definition: PDSS.h:486
doublereal m_a1
Input a1 coefficient (cal gmol-1 bar-1)
Definition: PDSS_HKFT.h:304
doublereal m_mw
Molecular Weight of the species.
Definition: PDSS.h:495