Cantera  2.0
ShomateThermo.h
Go to the documentation of this file.
1 /**
2  * @file ShomateThermo.h
3  * Header for the 2 regions Shomate polynomial
4  * for multiple species in a phase, derived from the
5  * \link Cantera::SpeciesThermo SpeciesThermo\endlink base class (see \ref mgrsrefcalc and
7  */
8 // Copyright 2001 California Institute of Technology
9
10
11 #ifndef CT_SHOMATETHERMO_H
12 #define CT_SHOMATETHERMO_H
13
15 #include "ShomatePoly.h"
17
18 namespace Cantera
19 {
20
21 //! A species thermodynamic property manager for the Shomate polynomial parameterization.
22 /*!
23  * This is the parameterization used
24  * in the NIST Chemistry WebBook (http://webbook.nist.gov/chemistry)
25  * The parameterization assumes there are two temperature regions
26  * each with its own Shomate polynomial representation, for each
27  * species in the phase.
28  *
29  * \f[
30  * \tilde{c}_p^0(T) = A + B t + C t^2 + D t^3 + \frac{E}{t^2}
31  * \f]
32  * \f[
33  * \tilde{h}^0(T) = A t + \frac{B t^2}{2} + \frac{C t^3}{3}
34  + \frac{D t^4}{4} - \frac{E}{t} + F.
35  * \f]
36  * \f[
37  * \tilde{s}^0(T) = A\ln t + B t + \frac{C t^2}{2}
38  + \frac{D t^3}{3} - \frac{E}{2t^2} + G.
39  * \f]
40  *
41  * In the above expressions, the thermodynamic polynomials are expressed
42  * in dimensional units, but the temperature,\f$t \f$, is divided by 1000. The
43  * following dimensions are assumed in the above expressions:
44  *
45  * - \f$\tilde{c}_p^0(T)\f$ = Heat Capacity (J/gmol*K)
46  * - \f$\tilde{h}^0(T) \f$ = standard Enthalpy (kJ/gmol)
47  * - \f$\tilde{s}^0(T) \f$= standard Entropy (J/gmol*K)
48  * - \f$t \f$= temperature (K) / 1000.
49  *
50  * Note, the polynomial data (i.e., A, ... , G) is entered in dimensional
51  * form.
52  *
53  * This is in contrast to the NASA database polynomials which are entered in
54  * nondimensional form (i.e., NASA parameterizes C_p/R, while Shomate
55  * parameterizes C_p assuming units of J/gmol*K - and kJ/gmol*K for H).
56  * Note, also that the H - H_298.15 equation has units of kJ/gmol, because of
57  * the implicit integration of (t = T 1000), which provides a
58  * multiplier of 1000 to the Enthalpy equation.
59  *
60  * @ingroup mgrsrefcalc
61  */
63 {
64
65 public:
66
67  //! Initialized to the type of parameterization
68  /*!
69  * Note, this value is used in some template functions
70  */
71  const int ID;
72
73  //! constructor
75  ID(SHOMATE),
76  m_tlow_max(0.0),
77  m_thigh_min(1.e30),
78  m_p0(-1.0),
79  m_ngroups(0) {
80  m_t.resize(7);
81  }
82
83  //! destructor
84  virtual ~ShomateThermo() {}
85
86  //! Copy Constructor
87  /*!
88  * @param right Object to be copied
89  */
90  ShomateThermo(const ShomateThermo& right) :
91  ID(SHOMATE),
92  m_tlow_max(0.0),
93  m_thigh_min(1.e30),
94  m_p0(-1.0),
95  m_ngroups(0) {
96  *this = operator=(right);
97  }
98
99  //! Assignment Operator
100  /*!
101  * @param right Object to be copied
102  */
104  if (&right == this) {
105  return *this;
106  }
107
108  m_high = right.m_high;
109  m_low = right.m_low;
110  m_index = right.m_index;
111  m_tmid = right.m_tmid;
112  m_tlow_max = right.m_tlow_max;
113  m_thigh_min = right.m_thigh_min;
114  m_tlow = right.m_tlow;
115  m_thigh = right.m_thigh;
116  m_p0 = right.m_p0;
117  m_ngroups = right.m_ngroups;
118  m_t = right.m_t;
119  m_group_map = right.m_group_map;
121
122  return *this;
123  }
124
125
126  //! Duplication routine for objects which inherit from
127  //! %SpeciesThermo
128  /*!
129  * This virtual routine can be used to duplicate %SpeciesThermo objects
130  * inherited from %SpeciesThermo even if the application only has
131  * a pointer to %SpeciesThermo to work with.
132  * ->commented out because we first need to add copy constructors
133  * and assignment operators to all of the derived classes.
134  */
136  ShomateThermo* st = new ShomateThermo(*this);
137  return (SpeciesThermo*) st;
138  }
139
140  //! Install a new species thermodynamic property
141  //! parameterization for one species using Shomate polynomials
142  //!
143  /*!
144  * Two temperature regions are assumed.
145  *
146  * @param name Name of the species
147  * @param index Species index
148  * @param type int flag specifying the type of parameterization to be
149  * installed.
150  * @param c Vector of coefficients for the parameterization.
151  * There are 15 coefficients for the 2-zone Shomate polynomial.
152  * The first coefficient is the value of Tmid. The next 7
153  * coefficients are the low temperature range Shomate coefficients.
154  * The last 7 are the high temperature range Shomate coefficients.
155  *
156  * @param minTemp minimum temperature for which this parameterization
157  * is valid.
158  * @param maxTemp maximum temperature for which this parameterization
159  * is valid.
160  * @param refPressure standard-state pressure for this
161  * parameterization.
162  *
163  * @see ShomatePoly
164  * @see ShomatePoly2
165  */
166  virtual void install(std::string name, size_t index, int type,
167  const doublereal* c,
168  doublereal minTemp, doublereal maxTemp,
169  doublereal refPressure) {
170  int imid = int(c[0]); // midpoint temp converted to integer
171  int igrp = m_index[imid]; // has this value been seen before?
172  if (igrp == 0) { // if not, prepare new group
173  std::vector<ShomatePoly> v;
174  m_high.push_back(v);
175  m_low.push_back(v);
176  m_tmid.push_back(c[0]);
177  m_index[imid] = igrp = static_cast<int>(m_high.size());
178  m_ngroups++;
179  }
180  m_group_map[index] = igrp;
181  m_posInGroup_map[index] = (int) m_low[igrp-1].size();
182  doublereal tlow = minTemp;
183  doublereal tmid = c[0];
184  doublereal thigh = maxTemp;
185
186  const doublereal* clow = c + 1;
187  const doublereal* chigh = c + 8;
188  m_high[igrp-1].push_back(ShomatePoly(index, tmid, thigh,
189  refPressure, chigh));
190  m_low[igrp-1].push_back(ShomatePoly(index, tlow, tmid,
191  refPressure, clow));
192  if (tlow > m_tlow_max) {
193  m_tlow_max = tlow;
194  }
195  if (thigh < m_thigh_min) {
196  m_thigh_min = thigh;
197  }
198
199  if (m_tlow.size() < index + 1) {
200  m_tlow.resize(index + 1, tlow);
201  m_thigh.resize(index + 1, thigh);
202  }
203  m_tlow[index] = tlow;
204  m_thigh[index] = thigh;
205
206  if (m_p0 < 0.0) {
207  m_p0 = refPressure;
208  } else if (fabs(m_p0 - refPressure) > 0.1) {
209  std::string logmsg = " ERROR ShomateThermo: New Species, " + name
210  + ", has a different reference pressure, "
211  + fp2str(refPressure) + ", than existing reference pressure, " + fp2str(m_p0) + "\n";
212  writelog(logmsg);
213  logmsg = " This is now a fatal error\n";
214  writelog(logmsg);
215  throw CanteraError("install()", "Species have different reference pressures");
216  }
217  m_p0 = refPressure;
218
219  }
220
221  //! Install a new species thermodynamic property
222  //! parameterization for one species.
223  /*!
224  * @param stit_ptr Pointer to the SpeciesThermoInterpType object
225  * This will set up the thermo for one species
226  */
227  virtual void install_STIT(SpeciesThermoInterpType* stit_ptr) {
228  throw CanteraError("install_STIT", "not implemented");
229  }
230
231  //! Like update(), but only updates the single species k.
232  /*!
233  * @param k species index
234  * @param t Temperature (Kelvin)
235  * @param cp_R Vector of Dimensionless heat capacities.
236  * (length m_kk).
237  * @param h_RT Vector of Dimensionless enthalpies.
238  * (length m_kk).
239  * @param s_R Vector of Dimensionless entropies.
240  * (length m_kk).
241  */
242  virtual void update_one(size_t k, doublereal t, doublereal* cp_R,
243  doublereal* h_RT, doublereal* s_R) const {
244
245  doublereal tt = 1.e-3*t;
246  m_t[0] = tt;
247  m_t[1] = tt*tt;
248  m_t[2] = m_t[1]*tt;
249  m_t[3] = 1.0/m_t[1];
250  m_t[4] = log(tt);
251  m_t[5] = 1.0/GasConstant;
252  m_t[6] = 1.0/(GasConstant * t);
253
254  size_t grp = m_group_map[k];
255  size_t pos = m_posInGroup_map[k];
256  const std::vector<ShomatePoly> &mlg = m_low[grp-1];
257  const ShomatePoly* nlow = &(mlg[pos]);
258
259  doublereal tmid = nlow->maxTemp();
260  if (t < tmid) {
261  nlow->updateProperties(&m_t[0], cp_R, h_RT, s_R);
262  } else {
263  const std::vector<ShomatePoly> &mhg = m_high[grp-1];
264  const ShomatePoly* nhigh = &(mhg[pos]);
265  nhigh->updateProperties(&m_t[0], cp_R, h_RT, s_R);
266  }
267  }
268
269  //! Compute the reference-state properties for all species.
270  /*!
271  * Given temperature T in K, this method updates the values of
272  * the non-dimensional heat capacity at constant pressure,
273  * enthalpy, and entropy, at the reference pressure, Pref
274  * of each of the standard states.
275  *
276  * @param t Temperature (Kelvin)
277  * @param cp_R Vector of Dimensionless heat capacities.
278  * (length m_kk).
279  * @param h_RT Vector of Dimensionless enthalpies.
280  * (length m_kk).
281  * @param s_R Vector of Dimensionless entropies.
282  * (length m_kk).
283  */
284  virtual void update(doublereal t, doublereal* cp_R,
285  doublereal* h_RT, doublereal* s_R) const {
286  int i;
287
288  doublereal tt = 1.e-3*t;
289  m_t[0] = tt;
290  m_t[1] = tt*tt;
291  m_t[2] = m_t[1]*tt;
292  m_t[3] = 1.0/m_t[1];
293  m_t[4] = log(tt);
294  m_t[5] = 1.0/GasConstant;
295  m_t[6] = 1.0/(GasConstant * t);
296
297  std::vector<ShomatePoly>::const_iterator _begin, _end;
298  for (i = 0; i != m_ngroups; i++) {
299  if (t > m_tmid[i]) {
300  _begin = m_high[i].begin();
301  _end = m_high[i].end();
302  } else {
303  _begin = m_low[i].begin();
304  _end = m_low[i].end();
305  }
306  for (; _begin != _end; ++_begin) {
307  _begin->updateProperties(&m_t[0], cp_R, h_RT, s_R);
308  }
309  }
310  }
311
312  //! Minimum temperature.
313  /*!
314  * If no argument is supplied, this
315  * method returns the minimum temperature for which \e all
316  * parameterizations are valid. If an integer index k is
317  * supplied, then the value returned is the minimum
318  * temperature for species k in the phase.
319  *
320  * @param k Species index
321  */
322  virtual doublereal minTemp(size_t k=npos) const {
323  if (k == npos) {
324  return m_tlow_max;
325  } else {
326  return m_tlow[k];
327  }
328  }
329
330  //! Maximum temperature.
331  /*!
332  * If no argument is supplied, this
333  * method returns the maximum temperature for which \e all
334  * parameterizations are valid. If an integer index k is
335  * supplied, then the value returned is the maximum
336  * temperature for parameterization k.
337  *
338  * @param k species index
339  */
340  virtual doublereal maxTemp(size_t k=npos) const {
341  if (k == npos) {
342  return m_thigh_min;
343  } else {
344  return m_thigh[k];
345  }
346  }
347
348  //! The reference-state pressure for species k.
349  /*!
350  *
351  * returns the reference state pressure in Pascals for
352  * species k. If k is left out of the argument list,
353  * it returns the reference state pressure for the first
354  * species.
355  * Note that some SpeciesThermo implementations, such
356  * as those for ideal gases, require that all species
357  * in the same phase have the same reference state pressures.
358  *
359  * @param k species index
360  */
361  virtual doublereal refPressure(size_t k=npos) const {
362  return m_p0;
363  }
364
365  //! This utility function reports the type of parameterization
366  //! used for the species with index number index.
367  /*!
368  *
369  * @param index Species index
370  */
371  virtual int reportType(size_t index) const {
372  return SHOMATE;
373  }
374
375  /*!
376  * This utility function reports back the type of
377  * parameterization and all of the parameters for the
378  * species, index.
379  *
380  * @param index Species index
381  * @param type Integer type of the standard type
382  * @param c Vector of coefficients used to set the
383  * parameters for the standard state.
384  *
385  * @param minTemp output - Minimum temperature
386  * @param maxTemp output - Maximum temperature
387  * @param refPressure output - reference pressure (Pa).
388  */
389  virtual void reportParams(size_t index, int& type,
390  doublereal* const c,
391  doublereal& minTemp,
392  doublereal& maxTemp,
393  doublereal& refPressure) const {
394  type = reportType(index);
395  if (type == SHOMATE) {
396  size_t grp = m_group_map[index];
397  size_t pos = m_posInGroup_map[index];
398  int itype = SHOMATE;
399  const std::vector<ShomatePoly> &mlg = m_low[grp-1];
400  const std::vector<ShomatePoly> &mhg = m_high[grp-1];
401  const ShomatePoly* lowPoly = &(mlg[pos]);
402  const ShomatePoly* highPoly = &(mhg[pos]);
403  doublereal tmid = lowPoly->maxTemp();
404  c[0] = tmid;
405  size_t n;
406  double ttemp;
407  lowPoly->reportParameters(n, itype, minTemp, ttemp, refPressure,
408  c + 1);
409  if (n != index) {
410  throw CanteraError(" ", "confused");
411  }
412  if (itype != SHOMATE && itype != SHOMATE1) {
413  throw CanteraError(" ", "confused");
414  }
415  highPoly->reportParameters(n, itype, ttemp, maxTemp,
416  refPressure, c + 8);
417  if (n != index) {
418  throw CanteraError(" ", "confused");
419  }
420  if (itype != SHOMATE && itype != SHOMATE1) {
421  throw CanteraError(" ", "confused");
422  }
423  } else {
424  throw CanteraError(" ", "confused");
425  }
426  }
427
428  //! Modify parameters for the standard state
429  /*!
430  * @param index Species index
431  * @param c Vector of coefficients used to set the
432  * parameters for the standard state.
433  */
434  virtual void modifyParams(size_t index, doublereal* c) {
435  int type = reportType(index);
436  if (type == SHOMATE) {
437  size_t grp = m_group_map[index];
438  size_t pos = m_posInGroup_map[index];
439  std::vector<ShomatePoly> &mlg = m_low[grp-1];
440  std::vector<ShomatePoly> &mhg = m_high[grp-1];
441  ShomatePoly* lowPoly = &(mlg[pos]);
442  ShomatePoly* highPoly = &(mhg[pos]);
443  doublereal tmid = lowPoly->maxTemp();
444  if (fabs(c[0] - tmid) > 0.001) {
445  throw CanteraError("modifyParams", "can't change mid temp");
446  }
447
448  lowPoly->modifyParameters(c + 1);
449
450  highPoly->modifyParameters(c + 8);
451
452  } else {
453  throw CanteraError(" ", "confused");
454  }
455  }
456
457 #ifdef H298MODIFY_CAPABILITY
458
459  virtual doublereal reportOneHf298(int k) const {
460  doublereal h;
461  doublereal t = 298.15;
462
463  int grp = m_group_map[k];
464  int pos = m_posInGroup_map[k];
465  const std::vector<ShomatePoly> &mlg = m_low[grp-1];
466  const ShomatePoly* nlow = &(mlg[pos]);
467
468  doublereal tmid = nlow->maxTemp();
469  if (t <= tmid) {
470  h = nlow->reportHf298();
471  } else {
472  const std::vector<ShomatePoly> &mhg = m_high[grp-1];
473  const ShomatePoly* nhigh = &(mhg[pos]);
474  h = nhigh->reportHf298();
475  }
476  return h;
477  }
478
479  virtual void modifyOneHf298(const int k, const doublereal Hf298New) {
480
481  int grp = m_group_map[k];
482  int pos = m_posInGroup_map[k];
483  std::vector<ShomatePoly> &mlg = m_low[grp-1];
484  ShomatePoly* nlow = &(mlg[pos]);
485  std::vector<ShomatePoly> &mhg = m_high[grp-1];
486  ShomatePoly* nhigh = &(mhg[pos]);
487  doublereal tmid = nlow->maxTemp();
488
489  double hnow = reportOneHf298(k);
490  double delH = Hf298New - hnow;
491  if (298.15 <= tmid) {
492  nlow->modifyOneHf298(k, Hf298New);
493  double h = nhigh->reportHf298(0);
494  double hnew = h + delH;
495  nhigh->modifyOneHf298(k, hnew);
496  } else {
497  nhigh->modifyOneHf298(k, Hf298New);
498  double h = nlow->reportHf298(0);
499  double hnew = h + delH;
500  nlow->modifyOneHf298(k, hnew);
501  }
502
503  }
504
505
506 #endif
507 protected:
508
509  //! Vector of vector of NasaPoly1's for the high temp region.
510  /*!
511  * This is the high temp region representation.
512  * The first Length is equal to the number of groups.
513  * The second vector is equal to the number of species
514  * in that particular group.
515  */
516  std::vector<std::vector<ShomatePoly> > m_high;
517
518  //! Vector of vector of NasaPoly1's for the low temp region.
519  /*!
520  * This is the low temp region representation.
521  * The first Length is equal to the number of groups.
522  * The second vector is equal to the number of species
523  * in that particular group.
524  */
525  std::vector<std::vector<ShomatePoly> > m_low;
526
527  //! Map between the midpoint temperature, as an int, to the group number
528  /*!
529  * Length is equal to the number of groups. Only used in the setup.
530  */
531  std::map<int, int> m_index;
532
533  //! Vector of log temperature limits
534  /*!
535  * Length is equal to the number of groups.
536  */
538
539  //! Maximum value of the low temperature limit
540  doublereal m_tlow_max;
541
542  //! Minimum value of the high temperature limit
543  doublereal m_thigh_min;
544
545  //! Vector of low temperature limits (species index)
546  /*!
547  * Length is equal to number of species
548  */
550
551  //! Vector of low temperature limits (species index)
552  /*!
553  * Length is equal to number of species
554  */
556
557  //! Reference pressure (Pa)
558  /*!
559  * all species must have the same reference pressure.
560  */
561  doublereal m_p0;
562
563  //! number of groups
565
566  //! Vector of temperature polynomials
567  mutable vector_fp m_t;
568
569  /*!
570  * This map takes as its index, the species index in the phase.
571  * It returns the group index, where the temperature polynomials
572  * for that species are stored. group indices start at 1,
573  * so a decrement is always performed to access vectors.
574  */
575  mutable std::map<size_t, size_t> m_group_map;
576
577  /*!
578  * This map takes as its index, the species index in the phase.
579  * It returns the position index within the group, where the
580  * temperature polynomials for that species are stored.
581  */
582  mutable std::map<size_t, size_t> m_posInGroup_map;
583 };
584
585 }
586
587 #endif