Cantera  2.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  * Copyright (2006) Sandia Corporation. Under the terms of
10  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
11  * U.S. Government retains certain rights in this software.
12  */
13 
14 #include "cantera/base/ctml.h"
15 #include "cantera/thermo/PDSS_HKFT.h"
16 #include "cantera/thermo/WaterProps.h"
17 #include "cantera/thermo/PDSS_Water.h"
18 #include "cantera/base/stringUtils.h"
19 
20 #include <cstdlib>
21 #include <fstream>
22 
23 using namespace std;
24 using namespace ctml;
25 
26 namespace Cantera
27 {
28 
29 /*
30  * Basic list of constructors and duplicators
31  */
32 PDSS_HKFT::PDSS_HKFT(VPStandardStateTP* tp, size_t spindex) :
33  PDSS(tp, spindex),
34  m_waterSS(0),
35  m_densWaterSS(-1.0),
36  m_waterProps(0),
37  m_born_coeff_j(-1.0),
38  m_r_e_j(-1.0),
39  m_deltaG_formation_tr_pr(0.0),
40  m_deltaH_formation_tr_pr(0.0),
41  m_Mu0_tr_pr(0.0),
42  m_Entrop_tr_pr(0.0),
43  m_a1(0.0),
44  m_a2(0.0),
45  m_a3(0.0),
46  m_a4(0.0),
47  m_c1(0.0),
48  m_c2(0.0),
49  m_omega_pr_tr(0.0),
50  m_Y_pr_tr(0.0),
51  m_Z_pr_tr(0.0),
52  m_presR_bar(0.0),
53  m_domega_jdT_prtr(0.0),
54  m_charge_j(0.0)
55 {
56  m_pres = OneAtm;
57  m_pdssType = cPDSS_MOLAL_HKFT;
58  m_presR_bar = OneAtm * 1.0E-5;
59 }
60 
61 
62 PDSS_HKFT::PDSS_HKFT(VPStandardStateTP* tp, size_t spindex, std::string inputFile, std::string id) :
63  PDSS(tp, spindex),
64  m_waterSS(0),
65  m_densWaterSS(-1.0),
66  m_waterProps(0),
67  m_born_coeff_j(-1.0),
68  m_r_e_j(-1.0),
69  m_deltaG_formation_tr_pr(0.0),
70  m_deltaH_formation_tr_pr(0.0),
71  m_Mu0_tr_pr(0.0),
72  m_Entrop_tr_pr(0.0),
73  m_a1(0.0),
74  m_a2(0.0),
75  m_a3(0.0),
76  m_a4(0.0),
77  m_c1(0.0),
78  m_c2(0.0),
79  m_omega_pr_tr(0.0),
80  m_Y_pr_tr(0.0),
81  m_Z_pr_tr(0.0),
82  m_presR_bar(0.0),
83  m_domega_jdT_prtr(0.0),
84  m_charge_j(0.0)
85 {
86  m_pres = OneAtm;
87  m_pdssType = cPDSS_MOLAL_HKFT;
88  m_presR_bar = OneAtm * 1.0E-5;
89  constructPDSSFile(tp, spindex, inputFile, id);
90 }
91 
92 PDSS_HKFT::PDSS_HKFT(VPStandardStateTP* tp, size_t spindex, const XML_Node& speciesNode,
93  const XML_Node& phaseRoot, bool spInstalled) :
94  PDSS(tp, spindex),
95  m_waterSS(0),
96  m_densWaterSS(-1.0),
97  m_waterProps(0),
98  m_born_coeff_j(-1.0),
99  m_r_e_j(-1.0),
100  m_deltaG_formation_tr_pr(0.0),
101  m_deltaH_formation_tr_pr(0.0),
102  m_Mu0_tr_pr(0.0),
103  m_Entrop_tr_pr(0.0),
104  m_a1(0.0),
105  m_a2(0.0),
106  m_a3(0.0),
107  m_a4(0.0),
108  m_c1(0.0),
109  m_c2(0.0),
110  m_omega_pr_tr(0.0),
111  m_Y_pr_tr(0.0),
112  m_Z_pr_tr(0.0),
113  m_presR_bar(0.0),
114  m_domega_jdT_prtr(0.0),
115  m_charge_j(0.0)
116 {
117  m_pres = OneAtm;
118  m_pdssType = cPDSS_MOLAL_HKFT;
119  m_presR_bar = OneAtm * 1.0E-5;
120  // We have to read the info from here
121  constructPDSSXML(tp, spindex, speciesNode, phaseRoot, spInstalled);
122 }
123 
124 PDSS_HKFT::PDSS_HKFT(const PDSS_HKFT& b) :
125  PDSS(b),
126  m_waterSS(0),
127  m_densWaterSS(-1.0),
128  m_waterProps(0),
129  m_born_coeff_j(-1.0),
130  m_r_e_j(-1.0),
131  m_deltaG_formation_tr_pr(0.0),
132  m_deltaH_formation_tr_pr(0.0),
133  m_Mu0_tr_pr(0.0),
134  m_Entrop_tr_pr(0.0),
135  m_a1(0.0),
136  m_a2(0.0),
137  m_a3(0.0),
138  m_a4(0.0),
139  m_c1(0.0),
140  m_c2(0.0),
141  m_omega_pr_tr(0.0),
142  m_Y_pr_tr(0.0),
143  m_Z_pr_tr(0.0),
144  m_presR_bar(0.0),
145  m_domega_jdT_prtr(0.0),
146  m_charge_j(0.0)
147 {
148  m_pdssType = cPDSS_MOLAL_HKFT;
149  m_presR_bar = OneAtm * 1.0E-5;
150  /*
151  * Use the assignment operator to do the brunt
152  * of the work for the copy constructor.
153  */
154  *this = b;
155 }
156 
157 /*
158  * Assignment operator
159  */
160 PDSS_HKFT& PDSS_HKFT::operator=(const PDSS_HKFT& b)
161 {
162  if (&b == this) {
163  return *this;
164  }
165  /*
166  * Call the base class operator
167  */
168  PDSS::operator=(b);
169 
170  //! Need to call initAllPtrs AFTER, to get the correct m_waterSS
171  m_waterSS = 0;
172  m_densWaterSS = b.m_densWaterSS;
173  //! Need to call initAllPtrs AFTER, to get the correct m_waterProps
174  if (m_waterProps) {
175  delete m_waterProps;
176  }
177  m_waterProps = 0;
178  m_born_coeff_j = b.m_born_coeff_j;
179  m_r_e_j = b.m_r_e_j;
180  m_deltaG_formation_tr_pr = b.m_deltaG_formation_tr_pr;
181  m_deltaH_formation_tr_pr = b.m_deltaH_formation_tr_pr;
182  m_Mu0_tr_pr = b.m_Mu0_tr_pr;
183  m_Entrop_tr_pr = b.m_Entrop_tr_pr;
184  m_a1 = b.m_a1;
185  m_a2 = b.m_a2;
186  m_a3 = b.m_a3;
187  m_a4 = b.m_a4;
188  m_c1 = b.m_c1;
189  m_c2 = b.m_c2;
190  m_omega_pr_tr = b.m_omega_pr_tr;
191  m_Y_pr_tr = b.m_Y_pr_tr;
192  m_Z_pr_tr = b.m_Z_pr_tr;
193  m_presR_bar = b.m_presR_bar;
194  m_domega_jdT_prtr = b.m_domega_jdT_prtr;
195  m_charge_j = b.m_charge_j;
196 
197  // Here we just fill these in so that local copies within the VPSS object work.
198  m_waterSS = b.m_waterSS;
199  m_waterProps = new WaterProps(m_waterSS);
200 
201  return *this;
202 }
203 
204 /*
205  * Destructor for the PDSS_HKFT class
206  */
207 PDSS_HKFT::~PDSS_HKFT()
208 {
209  delete m_waterProps;
210 }
211 
212 // Duplicator
213 PDSS* PDSS_HKFT::duplMyselfAsPDSS() const
214 {
215  PDSS_HKFT* idg = new PDSS_HKFT(*this);
216  return (PDSS*) idg;
217 }
218 
219 /*
220  * Return the molar enthalpy in units of J kmol-1
221  */
222 doublereal PDSS_HKFT::enthalpy_mole() const
223 {
224  // Ok we may change this evaluation method in the future.
225  doublereal GG = gibbs_mole();
226  doublereal SS = entropy_mole();
227  doublereal h = GG + m_temp * SS;
228 
229 #ifdef DEBUG_MODE_NOT
230  doublereal h2 = enthalpy_mole2();
231  if (fabs(h - h2) > 1.0E-1) {
232  printf("we are here, h = %g, h2 = %g, k = %d, T = %g, P = %g p0 = %g\n",
233  h, h2, m_spindex, m_temp, m_pres,
234  m_p0);
235  }
236 #endif
237  return h;
238 }
239 
240 doublereal PDSS_HKFT::enthalpy_RT() const
241 {
242  doublereal hh = enthalpy_mole();
243  doublereal RT = GasConstant * m_temp;
244  return hh / RT;
245 }
246 
247 #ifdef DEBUG_MODE
248 doublereal PDSS_HKFT::enthalpy_mole2() const
249 {
250  doublereal delH = deltaH();
251  double enthTRPR = m_Mu0_tr_pr + 298.15 * m_Entrop_tr_pr * 1.0E3 * 4.184;
252  double res = delH + enthTRPR;
253  return res;
254 }
255 #endif
256 
257 /*
258  * Calculate the internal energy in mks units of
259  * J kmol-1
260  */
261 doublereal PDSS_HKFT::intEnergy_mole() const
262 {
263  doublereal hh = enthalpy_RT();
264  doublereal mv = molarVolume();
265  return (hh - mv * m_pres);
266 }
267 
268 /*
269  * Calculate the entropy in mks units of
270  * J kmol-1 K-1
271  */
272 doublereal PDSS_HKFT::entropy_mole() const
273 {
274  doublereal delS = deltaS();
275  return (m_Entrop_tr_pr * 1.0E3 * 4.184 + delS);
276 }
277 
278 /*
279  * Calculate the Gibbs free energy in mks units of
280  * J kmol-1
281  */
282 doublereal PDSS_HKFT::gibbs_mole() const
283 {
284  doublereal delG = deltaG();
285  return (m_Mu0_tr_pr + delG);
286 }
287 
288 /*
289  * Calculate the constant pressure heat capacity
290  * in mks units of J kmol-1 K-1
291  */
292 doublereal PDSS_HKFT::cp_mole() const
293 {
294 
295  doublereal pbar = m_pres * 1.0E-5;
296  doublereal c1term = m_c1;
297 
298  doublereal c2term = m_c2 / (m_temp - 228.) / (m_temp - 228.);
299 
300  doublereal a3term = -m_a3 / (m_temp - 228.) / (m_temp - 228.) / (m_temp - 228.) * 2.0 * m_temp * (pbar - m_presR_bar);
301 
302  doublereal a4term = -m_a4 / (m_temp - 228.) / (m_temp - 228.) / (m_temp - 228.) * 2.0 * m_temp
303  * log((2600. + pbar)/(2600. + m_presR_bar));
304 
305  doublereal omega_j;
306  doublereal domega_jdT;
307  doublereal d2omega_jdT2;
308  if (m_charge_j == 0.0) {
309  omega_j = m_omega_pr_tr;
310  domega_jdT = 0.0;
311  d2omega_jdT2 = 0.0;
312  } else {
313  doublereal nu = 166027;
314  doublereal r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
315 
316  doublereal gval = gstar(m_temp, m_pres, 0);
317 
318  doublereal dgvaldT = gstar(m_temp, m_pres, 1);
319  doublereal d2gvaldT2 = gstar(m_temp, m_pres, 2);
320 
321  doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
322  doublereal dr_e_jdT = fabs(m_charge_j) * dgvaldT;
323  doublereal d2r_e_jdT2 = fabs(m_charge_j) * d2gvaldT2;
324 
325  doublereal r_e_j2 = r_e_j * r_e_j;
326 
327  doublereal charge2 = m_charge_j * m_charge_j;
328 
329  doublereal r_e_H = 3.082 + gval;
330  doublereal r_e_H2 = r_e_H * r_e_H;
331 
332  omega_j = nu * (charge2 / r_e_j - m_charge_j / r_e_H);
333 
334  domega_jdT = nu * (-(charge2 / r_e_j2 * dr_e_jdT)
335  +(m_charge_j / r_e_H2 * dgvaldT));
336 
337  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
338  -2.0*m_charge_j*dgvaldT*dgvaldT/(r_e_H2*r_e_H) + m_charge_j*d2gvaldT2 /r_e_H2);
339  }
340 
341  doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
342  doublereal drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
343 
344  doublereal Y = drelepsilondT / (relepsilon * relepsilon);
345 
346  doublereal d2relepsilondT2 = m_waterProps->relEpsilon(m_temp, m_pres, 2);
347 
348 #ifdef DEBUG_MODE_NOT
349  doublereal d1 = m_waterProps->relEpsilon(m_temp, m_pres, 1);
350  doublereal d2 = m_waterProps->relEpsilon(m_temp + 0.0001, m_pres, 1);
351  doublereal d3 = (d2 - d1) / 0.0001;
352  if (fabs(d2relepsilondT2 - d3) > 1.0E-6) {
353  printf("we are here\n");
354  }
355 #endif
356 
357  doublereal X = d2relepsilondT2 / (relepsilon* relepsilon) - 2.0 * relepsilon * Y * Y;
358 
359  doublereal Z = -1.0 / relepsilon;
360 
361  doublereal yterm = 2.0 * m_temp * Y * domega_jdT;
362 
363  doublereal xterm = omega_j * m_temp * X;
364 
365  doublereal otterm = m_temp * d2omega_jdT2 * (Z + 1.0);
366 
367  doublereal rterm = - m_domega_jdT_prtr * (m_Z_pr_tr + 1.0);
368 
369  doublereal Cp_calgmol = c1term + c2term + a3term + a4term + yterm + xterm + otterm + rterm;
370 
371  // Convert to Joules / kmol
372  doublereal Cp = Cp_calgmol * 1.0E3 * 4.184;
373 
374 #ifdef DEBUG_MODE_NOT
375  double e1 = enthalpy_mole();
376  m_temp = m_temp - 0.001;
377  double e2 = enthalpy_mole();
378  m_temp = m_temp + 0.001;
379  double cpd = (e1 - e2) / 0.001;
380  if (fabs(Cp - cpd) > 10.0) {
381  printf("Cp difference : raw: %g, delta: %g, k = %d, T = %g, m_pres = %g\n",
382  Cp, cpd, m_spindex, m_temp, m_pres);
383  }
384 #endif
385  return Cp;
386 }
387 
388 /*
389  * Calculate the constant volume heat capacity
390  * in mks units of J kmol-1 K-1
391  */
392 doublereal
393 PDSS_HKFT::cv_mole() const
394 {
395  throw CanteraError("PDSS_HKFT::cv_mole()", "unimplemented");
396  return (0.0);
397 }
398 
399 doublereal PDSS_HKFT::molarVolume() const
400 {
401 
402  // Initially do all calculations in (cal/gmol/Pa)
403 
404  doublereal a1term = m_a1 * 1.0E-5;
405 
406  doublereal a2term = m_a2 / (2600.E5 + m_pres);
407 
408  doublereal a3term = m_a3 * 1.0E-5/ (m_temp - 228.);
409 
410  doublereal a4term = m_a4 / (m_temp - 228.) / (2600.E5 + m_pres);
411 
412  doublereal omega_j;
413  doublereal domega_jdP;
414  if (m_charge_j == 0.0) {
415  omega_j = m_omega_pr_tr;
416  domega_jdP = 0.0;
417  } else {
418  doublereal nu = 166027.;
419  doublereal charge2 = m_charge_j * m_charge_j;
420  doublereal r_e_j_pr_tr = charge2 / (m_omega_pr_tr/nu + m_charge_j/3.082);
421 
422  doublereal gval = gstar(m_temp, m_pres, 0);
423  doublereal dgvaldP = gstar(m_temp, m_pres, 3);
424 
425  doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
426  doublereal r_e_H = 3.082 + gval;
427 
428  omega_j = nu * (charge2 / r_e_j - m_charge_j / r_e_H);
429 
430  doublereal dr_e_jdP = fabs(m_charge_j) * dgvaldP;
431 
432  domega_jdP = - nu * (charge2 / (r_e_j * r_e_j) * dr_e_jdP)
433  + nu * m_charge_j / (r_e_H * r_e_H) * dgvaldP;
434  }
435 
436  doublereal drelepsilondP = m_waterProps->relEpsilon(m_temp, m_pres, 3);
437 
438  doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
439 
440  doublereal Q = drelepsilondP / (relepsilon * relepsilon);
441 
442  doublereal Z = -1.0 / relepsilon;
443 
444  doublereal wterm = - domega_jdP * (Z + 1.0);
445 
446  doublereal qterm = - omega_j * Q;
447 
448  doublereal molVol_calgmolPascal = a1term + a2term + a3term + a4term + wterm + qterm;
449 
450  // Convert to m**3 / kmol from (cal/gmol/Pa)
451  doublereal molVol = molVol_calgmolPascal * 4.184 * 1.0E3;
452  return molVol;
453 }
454 
455 doublereal
456 PDSS_HKFT::density() const
457 {
458  doublereal val = molarVolume();
459  return (m_mw/val);
460 }
461 
462 doublereal
463 PDSS_HKFT::gibbs_RT_ref() const
464 {
465  doublereal m_psave = m_pres;
466  m_pres = m_waterSS->pref_safe(m_temp);
467  doublereal ee = gibbs_RT();
468  m_pres = m_psave;
469  return ee;
470 }
471 
472 doublereal
473 PDSS_HKFT::enthalpy_RT_ref() const
474 {
475  doublereal m_psave = m_pres;
476  m_pres = m_waterSS->pref_safe(m_temp);
477  doublereal hh = enthalpy_RT();
478  m_pres = m_psave;
479  return hh;
480 }
481 
482 doublereal
483 PDSS_HKFT::entropy_R_ref() const
484 {
485  doublereal m_psave = m_pres;
486  m_pres = m_waterSS->pref_safe(m_temp);
487  doublereal ee = entropy_R();
488  m_pres = m_psave;
489  return ee;
490 }
491 
492 doublereal
493 PDSS_HKFT::cp_R_ref() const
494 {
495  doublereal m_psave = m_pres;
496  m_pres = m_waterSS->pref_safe(m_temp);
497  doublereal ee = cp_R();
498  m_pres = m_psave;
499  return ee;
500 }
501 
502 doublereal
503 PDSS_HKFT::molarVolume_ref() const
504 {
505  doublereal m_psave = m_pres;
506  m_pres = m_waterSS->pref_safe(m_temp);
507  doublereal ee = molarVolume();
508  m_pres = m_psave;
509  return ee;
510 }
511 
512 /*
513  * Calculate the pressure (Pascals), given the temperature and density
514  * Temperature: kelvin
515  * rho: density in kg m-3
516  */
517 doublereal
518 PDSS_HKFT::pressure() const
519 {
520  return m_pres;
521 }
522 
523 void
524 PDSS_HKFT::setPressure(doublereal p)
525 {
526  m_pres = p;
527 }
528 
529 void PDSS_HKFT::setTemperature(doublereal temp)
530 {
531  m_temp = temp;
532 }
533 
534 doublereal PDSS_HKFT::temperature() const
535 {
536  return m_temp;
537 }
538 
539 void PDSS_HKFT::setState_TP(doublereal temp, doublereal pres)
540 {
541  setTemperature(temp);
542  setPressure(pres);
543 }
544 
545 // critical temperature
546 doublereal
547 PDSS_HKFT::critTemperature() const
548 {
549  throw CanteraError("PDSS_HKFT::critTemperature()", "unimplemented");
550  return (0.0);
551 }
552 
553 // critical pressure
554 doublereal PDSS_HKFT::critPressure() const
555 {
556  throw CanteraError("PDSS_HKFT::critPressure()", "unimplemented");
557  return (0.0);
558 }
559 
560 // critical density
561 doublereal PDSS_HKFT::critDensity() const
562 {
563  throw CanteraError("PDSS_HKFT::critDensity()", "unimplemented");
564  return (0.0);
565 }
566 
567 
568 void PDSS_HKFT::initThermo()
569 {
570  PDSS::initThermo();
571 
572  m_waterSS = (PDSS_Water*) m_tp->providePDSS(0);
573  /*
574  * Section to initialize m_Z_pr_tr and m_Y_pr_tr
575  */
576  m_temp = 273.15 + 25.;
577  m_pres = OneAtm;
578  doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
579 
580  m_waterSS->setState_TP(m_temp, m_pres);
581  m_densWaterSS = m_waterSS->density();
582  m_Z_pr_tr = -1.0 / relepsilon;
583 
584  doublereal drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
585 
586  m_Y_pr_tr = drelepsilondT / (relepsilon * relepsilon);
587 
588  m_waterProps = new WaterProps(m_waterSS);
589 
590  m_presR_bar = OneAtm / 1.0E5;
591  m_charge_j = m_tp->charge(m_spindex);
592  convertDGFormation();
593 
594  //! Ok, we have mu. Let's check it against the input value
595  // of DH_F to see that we have some internal consistency
596 
597  doublereal Hcalc = m_Mu0_tr_pr + 298.15 * (m_Entrop_tr_pr * 1.0E3 * 4.184);
598 
599  doublereal DHjmol = m_deltaH_formation_tr_pr * 1.0E3 * 4.184;
600 
601  // If the discrepency is greater than 100 cal gmol-1, print
602  // an error and exit.
603  if (fabs(Hcalc -DHjmol) > 100.* 1.0E3 * 4.184) {
604  throw CanteraError(" PDSS_HKFT::initThermo()",
605  "DHjmol is not consistent with G and S: " +
606  fp2str(Hcalc/(4.184E3)) + " vs "
607  + fp2str(m_deltaH_formation_tr_pr) + "cal gmol-1");
608  }
609  doublereal nu = 166027;
610 
611  doublereal r_e_j_pr_tr;
612  if (m_charge_j != 0.0) {
613  r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
614  } else {
615  r_e_j_pr_tr = 0.0;
616  }
617 
618  if (m_charge_j == 0.0) {
619  m_domega_jdT_prtr = 0.0;
620  } else {
621  doublereal gval = gstar(m_temp, m_pres, 0);
622 
623  doublereal dgvaldT = gstar(m_temp, m_pres, 1);
624 
625  doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
626  doublereal dr_e_jdT = fabs(m_charge_j) * dgvaldT;
627 
628  m_domega_jdT_prtr = - nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdT)
629  + nu * m_charge_j / (3.082 + gval) / (3.082 + gval) * dgvaldT;
630  }
631 }
632 
633 
634 void PDSS_HKFT::initThermoXML(const XML_Node& phaseNode, std::string& id)
635 {
636  PDSS::initThermoXML(phaseNode, id);
637 }
638 
639 void PDSS_HKFT::initAllPtrs(VPStandardStateTP* vptp_ptr, VPSSMgr* vpssmgr_ptr,
640  SpeciesThermo* spthermo_ptr)
641 {
642 
643  PDSS::initAllPtrs(vptp_ptr, vpssmgr_ptr, spthermo_ptr);
644  m_waterSS = (PDSS_Water*) m_tp->providePDSS(0);
645  if (m_waterProps) {
646  delete m_waterProps;
647  }
648  m_waterProps = new WaterProps(m_waterSS);
649 }
650 
651 void PDSS_HKFT::constructPDSSXML(VPStandardStateTP* tp, size_t spindex,
652  const XML_Node& speciesNode,
653  const XML_Node& phaseNode, bool spInstalled)
654 {
655  int hasDGO = 0;
656  int hasSO = 0;
657  int hasDHO = 0;
658 
659  if (!spInstalled) {
660  throw CanteraError("PDSS_HKFT::constructPDSSXML", "spInstalled false not handled");
661  }
662 
663  const XML_Node* tn = speciesNode.findByName("thermo");
664  if (!tn) {
665  throw CanteraError("PDSS_HKFT::constructPDSSXML",
666  "no thermo Node for species " + speciesNode.name());
667  }
668  std::string model = lowercase((*tn)["model"]);
669  if (model != "hkft") {
670  throw CanteraError("PDSS_HKFT::initThermoXML",
671  "thermo model for species isn't hkft: "
672  + speciesNode.name());
673  }
674  const XML_Node* hh = tn->findByName("HKFT");
675  if (!hh) {
676  throw CanteraError("PDSS_HKFT::constructPDSSXML",
677  "no Thermo::HKFT Node for species " + speciesNode.name());
678  }
679 
680  // go get the attributes
681  m_p0 = OneAtm;
682  std::string p0string = (*hh)["Pref"];
683  if (p0string != "") {
684  m_p0 = strSItoDbl(p0string);
685  }
686 
687  std::string minTstring = (*hh)["Tmin"];
688  if (minTstring != "") {
689  m_minTemp = atofCheck(minTstring.c_str());
690  }
691 
692  std::string maxTstring = (*hh)["Tmax"];
693  if (maxTstring != "") {
694  m_maxTemp = atofCheck(maxTstring.c_str());
695  }
696 
697  if (hh->hasChild("DG0_f_Pr_Tr")) {
698  doublereal val = getFloat(*hh, "DG0_f_Pr_Tr");
699  m_deltaG_formation_tr_pr = val;
700  hasDGO = 1;
701  } else {
702  // throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing DG0_f_Pr_Tr field");
703  }
704 
705  if (hh->hasChild("DH0_f_Pr_Tr")) {
706  doublereal val = getFloat(*hh, "DH0_f_Pr_Tr");
707  m_deltaH_formation_tr_pr = val;
708  hasDHO = 1;
709  } else {
710  // throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing DH0_f_Pr_Tr field");
711  }
712 
713  if (hh->hasChild("S0_Pr_Tr")) {
714  doublereal val = getFloat(*hh, "S0_Pr_Tr");
715  m_Entrop_tr_pr= val;
716  hasSO = 1;
717  } else {
718  // throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing S0_Pr_Tr field");
719  }
720 
721  const XML_Node* ss = speciesNode.findByName("standardState");
722  if (!ss) {
723  throw CanteraError("PDSS_HKFT::constructPDSSXML",
724  "no standardState Node for species " + speciesNode.name());
725  }
726  model = lowercase((*ss)["model"]);
727  if (model != "hkft") {
728  throw CanteraError("PDSS_HKFT::initThermoXML",
729  "standardState model for species isn't hkft: "
730  + speciesNode.name());
731  }
732  if (ss->hasChild("a1")) {
733  doublereal val = getFloat(*ss, "a1");
734  m_a1 = val;
735  } else {
736  throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a1 field");
737  }
738  if (ss->hasChild("a2")) {
739  doublereal val = getFloat(*ss, "a2");
740  m_a2 = val;
741  } else {
742  throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a2 field");
743  }
744  if (ss->hasChild("a3")) {
745  doublereal val = getFloat(*ss, "a3");
746  m_a3 = val;
747  } else {
748  throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a3 field");
749  }
750  if (ss->hasChild("a4")) {
751  doublereal val = getFloat(*ss, "a4");
752  m_a4 = val;
753  } else {
754  throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a4 field");
755  }
756 
757  if (ss->hasChild("c1")) {
758  doublereal val = getFloat(*ss, "c1");
759  m_c1 = val;
760  } else {
761  throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing c1 field");
762  }
763  if (ss->hasChild("c2")) {
764  doublereal val = getFloat(*ss, "c2");
765  m_c2 = val;
766  } else {
767  throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing c2 field");
768  }
769  if (ss->hasChild("omega_Pr_Tr")) {
770  doublereal val = getFloat(*ss, "omega_Pr_Tr");
771  m_omega_pr_tr = val;
772  } else {
773  throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing omega_Pr_Tr field");
774  }
775 
776 
777  int isum = hasDGO + hasDHO + hasSO;
778  if (isum < 2) {
779  throw CanteraError("PDSS_HKFT::constructPDSSXML",
780  "Missing 2 or more of DG0_f_Pr_Tr, DH0_f_Pr_Tr, or S0_f_Pr_Tr fields. "
781  "Need to supply at least two of these fields");
782  }
783  // Ok, if we are missing one, then we construct its value from the other two.
784  // This code has been internally verified.
785  if (hasDHO == 0) {
786  m_charge_j = m_tp->charge(m_spindex);
787  convertDGFormation();
788  doublereal Hcalc = m_Mu0_tr_pr + 298.15 * (m_Entrop_tr_pr * 1.0E3 * 4.184);
789  m_deltaH_formation_tr_pr = Hcalc / (1.0E3 * 4.184);
790  }
791  if (hasDGO == 0) {
792  doublereal DHjmol = m_deltaH_formation_tr_pr * 1.0E3 * 4.184;
793  m_Mu0_tr_pr = DHjmol - 298.15 * (m_Entrop_tr_pr * 1.0E3 * 4.184);
794  m_deltaG_formation_tr_pr = m_Mu0_tr_pr / (1.0E3 * 4.184);
795  double tmp = m_Mu0_tr_pr;
796  m_charge_j = m_tp->charge(m_spindex);
797  convertDGFormation();
798  double totalSum = m_Mu0_tr_pr - tmp;
799  m_Mu0_tr_pr = tmp;
800  m_deltaG_formation_tr_pr = (m_Mu0_tr_pr - totalSum)/ (1.0E3 * 4.184);
801  }
802  if (hasSO == 0) {
803  m_charge_j = m_tp->charge(m_spindex);
804  convertDGFormation();
805  doublereal DHjmol = m_deltaH_formation_tr_pr * 1.0E3 * 4.184;
806  m_Entrop_tr_pr = (DHjmol - m_Mu0_tr_pr) / (298.15 * 1.0E3 * 4.184);
807  }
808 
809 }
810 
811 void PDSS_HKFT::constructPDSSFile(VPStandardStateTP* tp, size_t spindex,
812  std::string inputFile, std::string id)
813 {
814 
815  if (inputFile.size() == 0) {
816  throw CanteraError("PDSS_HKFT::initThermo",
817  "input file is null");
818  }
819  std::string path = findInputFile(inputFile);
820  ifstream fin(path.c_str());
821  if (!fin) {
822  throw CanteraError("PDSS_HKFT::initThermo","could not open "
823  +path+" for reading.");
824  }
825  /*
826  * The phase object automatically constructs an XML object.
827  * Use this object to store information.
828  */
829 
830  XML_Node* fxml = new XML_Node();
831  fxml->build(fin);
832  XML_Node* fxml_phase = findXMLPhase(fxml, id);
833  if (!fxml_phase) {
834  throw CanteraError("PDSS_HKFT::initThermo",
835  "ERROR: Can not find phase named " +
836  id + " in file named " + inputFile);
837  }
838 
839  XML_Node& speciesList = fxml_phase->child("speciesArray");
840  XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"],
841  &(fxml_phase->root()));
842  const vector<string>&sss = tp->speciesNames();
843  const XML_Node* s = speciesDB->findByAttr("name", sss[spindex]);
844 
845  constructPDSSXML(tp, spindex, *s, *fxml_phase, true);
846  delete fxml;
847 }
848 
849 #ifdef DEBUG_MODE
850 doublereal PDSS_HKFT::deltaH() const
851 {
852 
853  doublereal pbar = m_pres * 1.0E-5;
854 
855  doublereal c1term = m_c1 * (m_temp - 298.15);
856 
857  doublereal a1term = m_a1 * (pbar - m_presR_bar);
858 
859  doublereal a2term = m_a2 * log((2600. + pbar)/(2600. + m_presR_bar));
860 
861  doublereal c2term = -m_c2 * (1.0/(m_temp - 228.) - 1.0/(298.15 - 228.));
862 
863  double a3tmp = (2.0 * m_temp - 228.)/ (m_temp - 228.) /(m_temp - 228.);
864 
865  doublereal a3term = m_a3 * a3tmp * (pbar - m_presR_bar);
866 
867  doublereal a4term = m_a4 * a3tmp * log((2600. + pbar)/(2600. + m_presR_bar));
868 
869  doublereal omega_j;
870  doublereal domega_jdT;
871  if (m_charge_j == 0.0) {
872  omega_j = m_omega_pr_tr;
873  domega_jdT = 0.0;
874  } else {
875  doublereal nu = 166027;
876  doublereal r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
877  doublereal gval = gstar(m_temp, m_pres, 0);
878  doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
879 
880  doublereal dgvaldT = gstar(m_temp, m_pres, 1);
881  doublereal dr_e_jdT = fabs(m_charge_j) * dgvaldT;
882 
883  omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval));
884 
885  domega_jdT = - nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdT)
886  + nu * m_charge_j / (3.082 + gval) / (3.082 + gval) * dgvaldT;
887  }
888 
889  doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
890  doublereal drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
891 
892  doublereal Y = drelepsilondT / (relepsilon * relepsilon);
893 
894  doublereal Z = -1.0 / relepsilon;
895 
896  doublereal yterm = m_temp * omega_j * Y;
897  doublereal yrterm = - 298.15 * m_omega_pr_tr * m_Y_pr_tr;
898 
899  doublereal wterm = - omega_j * (Z + 1.0);
900  doublereal wrterm = + m_omega_pr_tr * (m_Z_pr_tr + 1.0);
901 
902  doublereal otterm = m_temp * domega_jdT * (Z + 1.0);
903  doublereal otrterm = - m_temp * m_domega_jdT_prtr * (m_Z_pr_tr + 1.0);
904 
905  doublereal deltaH_calgmol = c1term + a1term + a2term + c2term + a3term + a4term
906  + yterm + yrterm + wterm + wrterm + otterm + otrterm;
907 
908  // Convert to Joules / kmol
909  doublereal deltaH = deltaH_calgmol * 1.0E3 * 4.184;
910  return deltaH;
911 }
912 #endif
913 
914 doublereal PDSS_HKFT::deltaG() const
915 {
916 
917  doublereal pbar = m_pres * 1.0E-5;
918  //doublereal m_presR_bar = OneAtm * 1.0E-5;
919 
920  doublereal sterm = - m_Entrop_tr_pr * (m_temp - 298.15);
921 
922  doublereal c1term = -m_c1 * (m_temp * log(m_temp/298.15) - (m_temp - 298.15));
923  doublereal a1term = m_a1 * (pbar - m_presR_bar);
924 
925  doublereal a2term = m_a2 * log((2600. + pbar)/(2600. + m_presR_bar));
926 
927  doublereal c2term = -m_c2 * ((1.0/(m_temp - 228.) - 1.0/(298.15 - 228.)) * (228. - m_temp)/228.
928  - m_temp / (228.*228.) * log((298.15*(m_temp-228.)) / (m_temp*(298.15-228.))));
929 
930  doublereal a3term = m_a3 / (m_temp - 228.) * (pbar - m_presR_bar);
931 
932  doublereal a4term = m_a4 / (m_temp - 228.) * log((2600. + pbar)/(2600. + m_presR_bar));
933 
934  doublereal omega_j;
935  if (m_charge_j == 0.0) {
936  omega_j = m_omega_pr_tr;
937  } else {
938  doublereal nu = 166027;
939  doublereal r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
940  doublereal gval = gstar(m_temp, m_pres, 0);
941  doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
942  omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval));
943  }
944 
945  doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
946 
947  doublereal Z = -1.0 / relepsilon;
948 
949  doublereal wterm = - omega_j * (Z + 1.0);
950 
951  doublereal wrterm = m_omega_pr_tr * (m_Z_pr_tr + 1.0);
952 
953  doublereal yterm = m_omega_pr_tr * m_Y_pr_tr * (m_temp - 298.15);
954 
955  doublereal deltaG_calgmol = sterm + c1term + a1term + a2term + c2term + a3term + a4term + wterm + wrterm + yterm;
956 
957  // Convert to Joules / kmol
958  doublereal deltaG = deltaG_calgmol * 1.0E3 * 4.184;
959  return deltaG;
960 }
961 
962 
963 doublereal PDSS_HKFT::deltaS() const
964 {
965 
966  doublereal pbar = m_pres * 1.0E-5;
967 
968  doublereal c1term = m_c1 * log(m_temp/298.15);
969 
970  doublereal c2term = -m_c2 / 228. * ((1.0/(m_temp - 228.) - 1.0/(298.15 - 228.))
971  + 1.0 / 228. * log((298.15*(m_temp-228.)) / (m_temp*(298.15-228.))));
972 
973  doublereal a3term = m_a3 / (m_temp - 228.) / (m_temp - 228.) * (pbar - m_presR_bar);
974 
975  doublereal a4term = m_a4 / (m_temp - 228.) / (m_temp - 228.) * log((2600. + pbar)/(2600. + m_presR_bar));
976 
977  doublereal omega_j;
978  doublereal domega_jdT;
979  if (m_charge_j == 0.0) {
980  omega_j = m_omega_pr_tr;
981  domega_jdT = 0.0;
982  } else {
983 
984  doublereal nu = 166027;
985  doublereal r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
986 
987  doublereal gval = gstar(m_temp, m_pres, 0);
988 
989  doublereal dgvaldT = gstar(m_temp, m_pres, 1);
990 
991  doublereal r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
992  doublereal dr_e_jdT = fabs(m_charge_j) * dgvaldT;
993 
994  omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval));
995 
996  domega_jdT = - nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdT)
997  + nu * m_charge_j / (3.082 + gval) / (3.082 + gval) * dgvaldT;
998  }
999 
1000  doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
1001  doublereal drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
1002 
1003  doublereal Y = drelepsilondT / (relepsilon * relepsilon);
1004 
1005  doublereal Z = -1.0 / relepsilon;
1006 
1007  doublereal wterm = omega_j * Y;
1008 
1009  doublereal wrterm = - m_omega_pr_tr * m_Y_pr_tr;
1010 
1011  doublereal otterm = domega_jdT * (Z + 1.0);
1012 
1013  doublereal otrterm = - m_domega_jdT_prtr * (m_Z_pr_tr + 1.0);
1014 
1015  doublereal deltaS_calgmol = c1term + c2term + a3term + a4term + wterm + wrterm + otterm + otrterm;
1016 
1017  // Convert to Joules / kmol
1018  doublereal deltaS = deltaS_calgmol * 1.0E3 * 4.184;
1019  return deltaS;
1020 }
1021 
1022 
1023 // Internal formula for the calculation of a_g()
1024 /*
1025  * The output of this is in units of Angstroms
1026  */
1027 doublereal PDSS_HKFT::ag(const doublereal temp, const int ifunc) const
1028 {
1029  static doublereal ag_coeff[3] = { -2.037662, 5.747000E-3, -6.557892E-6};
1030  if (ifunc == 0) {
1031  doublereal t2 = temp * temp;
1032  doublereal val = ag_coeff[0] + ag_coeff[1] * temp + ag_coeff[2] * t2;
1033  return val;
1034  } else if (ifunc == 1) {
1035  return ag_coeff[1] + ag_coeff[2] * 2.0 * temp;
1036  }
1037  if (ifunc != 2) {
1038  return 0.0;
1039  }
1040  return ag_coeff[2] * 2.0;
1041 }
1042 
1043 
1044 // Internal formula for the calculation of b_g()
1045 /*
1046  * the output of this is unitless
1047  */
1048 doublereal PDSS_HKFT::bg(const doublereal temp, const int ifunc) const
1049 {
1050  static doublereal bg_coeff[3] = { 6.107361, -1.074377E-2, 1.268348E-5};
1051  if (ifunc == 0) {
1052  doublereal t2 = temp * temp;
1053  doublereal val = bg_coeff[0] + bg_coeff[1] * temp + bg_coeff[2] * t2;
1054  return val;
1055  } else if (ifunc == 1) {
1056  return bg_coeff[1] + bg_coeff[2] * 2.0 * temp;
1057  }
1058  if (ifunc != 2) {
1059  return 0.0;
1060  }
1061  return bg_coeff[2] * 2.0;
1062 }
1063 
1064 
1065 doublereal PDSS_HKFT::f(const doublereal temp, const doublereal pres, const int ifunc) const
1066 {
1067 
1068  static doublereal af_coeff[3] = { 3.666666E1, -0.1504956E-9, 0.5107997E-13};
1069  doublereal TC = temp - 273.15;
1070  doublereal presBar = pres / 1.0E5;
1071 
1072  if (TC < 155.0) {
1073  return 0.0;
1074  }
1075  if (TC > 355.0) {
1076  TC = 355.0;
1077  }
1078  if (presBar > 1000.) {
1079  return 0.0;
1080  }
1081 
1082 
1083  doublereal T1 = (TC-155.0)/300.;
1084  doublereal fac1;
1085 
1086  doublereal p2 = (1000. - presBar) * (1000. - presBar);
1087  doublereal p3 = (1000. - presBar) * p2;
1088  doublereal p4 = p2 * p2;
1089  doublereal fac2 = af_coeff[1] * p3 + af_coeff[2] * p4;
1090  if (ifunc == 0) {
1091  fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
1092  return fac1 * fac2;
1093  } else if (ifunc == 1) {
1094  fac1 = (4.8 * pow(T1,3.8) + 16.0 * af_coeff[0] * pow(T1, 15.0)) / 300.;
1095  return fac1 * fac2;
1096  } else if (ifunc == 2) {
1097  fac1 = (4.8 * 3.8 * pow(T1,2.8) + 16.0 * 15.0 * af_coeff[0] * pow(T1, 14.0)) / (300. * 300.);
1098  return fac1 * fac2;
1099  } else if (ifunc == 3) {
1100  fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
1101  fac2 = - (3.0 * af_coeff[1] * p2 + 4.0 * af_coeff[2] * p3)/ 1.0E5;
1102  return fac1 * fac2;
1103  } else {
1104  throw CanteraError("HKFT_PDSS::gg", "unimplemented");
1105  }
1106  return 0.0;
1107 }
1108 
1109 
1110 doublereal PDSS_HKFT::g(const doublereal temp, const doublereal pres, const int ifunc) const
1111 {
1112  doublereal afunc = ag(temp, 0);
1113  doublereal bfunc = bg(temp, 0);
1114  m_waterSS->setState_TP(temp, pres);
1115  m_densWaterSS = m_waterSS->density();
1116  // density in gm cm-3
1117  doublereal dens = m_densWaterSS * 1.0E-3;
1118  doublereal gval = afunc * pow((1.0-dens), bfunc);
1119  if (dens >= 1.0) {
1120  return 0.0;
1121  }
1122  if (ifunc == 0) {
1123  return gval;
1124 
1125  } else if (ifunc == 1 || ifunc == 2) {
1126  doublereal afuncdT = ag(temp, 1);
1127  doublereal bfuncdT = bg(temp, 1);
1128  doublereal alpha = m_waterSS->thermalExpansionCoeff();
1129 
1130  doublereal fac1 = afuncdT * gval / afunc;
1131  doublereal fac2 = bfuncdT * gval * log(1.0 - dens);
1132  doublereal fac3 = gval * alpha * bfunc * dens / (1.0 - dens);
1133 
1134  doublereal dgdt = fac1 + fac2 + fac3;
1135  if (ifunc == 1) {
1136  return dgdt;
1137  }
1138 
1139  doublereal afuncdT2 = ag(temp, 2);
1140  doublereal bfuncdT2 = bg(temp, 2);
1141 
1142  doublereal dfac1dT = dgdt * afuncdT / afunc + afuncdT2 * gval / afunc
1143  - afuncdT * afuncdT * gval / (afunc * afunc);
1144 
1145  doublereal ddensdT = - alpha * dens;
1146  doublereal dfac2dT = bfuncdT2 * gval * log(1.0 - dens)
1147  + bfuncdT * dgdt * log(1.0 - dens)
1148  - bfuncdT * gval /(1.0 - dens) * ddensdT;
1149 
1150  doublereal dalphadT = m_waterSS->dthermalExpansionCoeffdT();
1151 
1152  doublereal dfac3dT = dgdt * alpha * bfunc * dens / (1.0 - dens)
1153  + gval * dalphadT * bfunc * dens / (1.0 - dens)
1154  + gval * alpha * bfuncdT * dens / (1.0 - dens)
1155  + gval * alpha * bfunc * ddensdT / (1.0 - dens)
1156  + gval * alpha * bfunc * dens / ((1.0 - dens) * (1.0 - dens)) * ddensdT;
1157 
1158  return dfac1dT + dfac2dT + dfac3dT;
1159 
1160  } else if (ifunc == 3) {
1161  doublereal beta = m_waterSS->isothermalCompressibility();
1162 
1163  doublereal dgdp = - bfunc * gval * dens * beta / (1.0 - dens);
1164 
1165  return dgdp;
1166  } else {
1167  throw CanteraError("HKFT_PDSS::g", "unimplemented");
1168  }
1169  return 0.0;
1170 }
1171 
1172 
1173 doublereal PDSS_HKFT::gstar(const doublereal temp, const doublereal pres, const int ifunc) const
1174 {
1175  doublereal gval = g(temp, pres, ifunc);
1176  doublereal fval = f(temp, pres, ifunc);
1177  double res = gval - fval;
1178 #ifdef DEBUG_MODE_NOT
1179  if (ifunc == 2) {
1180  double gval1 = g(temp, pres, 1);
1181  double fval1 = f(temp, pres, 1);
1182  double gval2 = g(temp + 0.001, pres, 1);
1183  double fval2 = f(temp + 0.001, pres, 1);
1184  double gvalT = (gval2 - gval1) / 0.001;
1185  double fvalT = (fval2 - fval1) / 0.001;
1186  if (fabs(gvalT - gval) > 1.0E-9) {
1187  printf("we are here\n");
1188  }
1189  if (fabs(fvalT - fval) > 1.0E-9) {
1190  printf("we are here\n");
1191  }
1192  // return gvalT - fvalT;
1193  }
1194 #endif
1195  return res;
1196 }
1197 
1198 //! Static function to look up Element Free Energies
1199 /*!
1200  *
1201  * This static function looks up the argument string in the
1202  * database above and returns the associated Gibbs Free energies.
1203  *
1204  * @param elemName String. Only the first 3 characters are significant
1205  *
1206  * @return
1207  * Return value contains the Gibbs free energy for that element
1208  *
1209  * @exception CanteraError
1210  * If a match is not found, a CanteraError is thrown as well
1211  */
1212 doublereal PDSS_HKFT::LookupGe(const std::string& elemName)
1213 {
1214  size_t iE = m_tp->elementIndex(elemName);
1215  if (iE == npos) {
1216  throw CanteraError("PDSS_HKFT::LookupGe", "element " + elemName + " not found");
1217  }
1218  doublereal geValue = m_tp->entropyElement298(iE);
1219  if (geValue == ENTROPY298_UNKNOWN) {
1220  throw CanteraError("PDSS_HKFT::LookupGe",
1221  "element " + elemName + " doesn not have a supplied entropy298");
1222  }
1223  geValue *= (-298.15);
1224  return geValue;
1225 }
1226 
1227 void PDSS_HKFT::convertDGFormation()
1228 {
1229  /*
1230  * Ok let's get the element compositions and conversion factors.
1231  */
1232  size_t ne = m_tp->nElements();
1233  doublereal na;
1234  doublereal ge;
1235  string ename;
1236 
1237  doublereal totalSum = 0.0;
1238  for (size_t m = 0; m < ne; m++) {
1239  na = m_tp->nAtoms(m_spindex, m);
1240  if (na > 0.0) {
1241  ename = m_tp->elementName(m);
1242  ge = LookupGe(ename);
1243  totalSum += na * ge;
1244  }
1245  }
1246  // Add in the charge
1247  if (m_charge_j != 0.0) {
1248  ename = "H";
1249  ge = LookupGe(ename);
1250  totalSum -= m_charge_j * ge;
1251  }
1252  // Ok, now do the calculation. Convert to joules kmol-1
1253  doublereal dg = m_deltaG_formation_tr_pr * 4.184 * 1.0E3;
1254  //! Store the result into an internal variable.
1255  m_Mu0_tr_pr = dg + totalSum;
1256 }
1257 
1258 // This utility function reports back the type of
1259 // parameterization and all of the parameters for the
1260 // species, index.
1261 /*
1262  *
1263  * @param index Species index
1264  * @param type Integer type of the standard type
1265  * @param c Vector of coefficients used to set the
1266  * parameters for the standard state.
1267  * @param minTemp output - Minimum temperature
1268  * @param maxTemp output - Maximum temperature
1269  * @param refPressure output - reference pressure (Pa).
1270  *
1271  */
1272 void PDSS_HKFT::reportParams(size_t& kindex, int& type,
1273  doublereal* const c,
1274  doublereal& minTemp,
1275  doublereal& maxTemp,
1276  doublereal& refPressure) const
1277 {
1278 
1279  // Fill in the first part
1280  PDSS::reportParams(kindex, type, c, minTemp, maxTemp,
1281  refPressure);
1282 
1283 
1284  c[0] = m_deltaG_formation_tr_pr;
1285  c[1] = m_deltaH_formation_tr_pr;
1286  c[2] = m_Mu0_tr_pr;
1287  c[3] = m_Entrop_tr_pr;
1288  c[4] = m_a1;
1289  c[5] = m_a2;
1290  c[6] = m_a3;
1291  c[7] = m_a4;
1292  c[8] = m_c1;
1293  c[9] = m_c2;
1294  c[10] = m_omega_pr_tr;
1295 
1296 }
1297 
1298 
1299 
1300 }
Generated by   doxygen 1.8.2