Cantera  2.0
ShomatePoly.h
Go to the documentation of this file.
1 /**
2  * @file ShomatePoly.h
3  * Header for a single-species standard state object derived
5  * on the Shomate temperature polynomial form applied to one temperature region
8  * Shomate polynomial expressions.
9  */
10 // Copyright 2001 California Institute of Technology
11
12
13 #ifndef CT_SHOMATEPOLY1_H
14 #define CT_SHOMATEPOLY1_H
15
17
18 namespace Cantera
19 {
20
21
22 //! The Shomate polynomial parameterization for one temperature range
23 //! for one species
24 /*!
25  *
26  * Seven coefficients \f$(A,\dots,G)\f$ are used to represent
27  * \f$c_p^0(T)\f$, \f$h^0(T)\f$, and \f$s^0(T) \f$ as
28  * polynomials in the temperature, \f$T \f$ :
29  *
30  * \f[
31  * \tilde{c}_p^0(T) = A + B t + C t^2 + D t^3 + \frac{E}{t^2}
32  * \f]
33  * \f[
34  * \tilde{h}^0(T) = A t + \frac{B t^2}{2} + \frac{C t^3}{3}
35  + \frac{D t^4}{4} - \frac{E}{t} + F.
36  * \f]
37  * \f[
38  * \tilde{s}^0(T) = A\ln t + B t + \frac{C t^2}{2}
39  + \frac{D t^3}{3} - \frac{E}{2t^2} + G.
40  * \f]
41  *
42  * In the above expressions, the thermodynamic polynomials are expressed
43  * in dimensional units, but the temperature,\f$t \f$, is divided by 1000. The
44  * following dimensions are assumed in the above expressions:
45  *
46  * - \f$\tilde{c}_p^0(T)\f$ = Heat Capacity (J/gmol*K)
47  * - \f$\tilde{h}^0(T) \f$ = standard Enthalpy (kJ/gmol)
48  * - \f$\tilde{s}^0(T) \f$= standard Entropy (J/gmol*K)
49  * - \f$t \f$= temperature (K) / 1000.
50  *
52  * http://webbook.nist.gov/
53  *
54  * Before being used within Cantera, the dimensions must be adjusted to those
55  * used by Cantera (i.e., Joules and kmol).
56  *
57
58  * @ingroup spthermo
59  */
61 {
62
63 public:
64
65  //! Empty constructor
67  : m_lowT(0.0), m_highT(0.0),
68  m_Pref(0.0), m_index(0) {}
69
70  //! Constructor used in templated instantiations
71  /*!
72  * @param n Species index
73  * @param tlow Minimum temperature
74  * @param thigh Maximum temperature
75  * @param pref reference pressure (Pa).
76  * @param coeffs Vector of coefficients used to set the
77  * parameters for the standard state for species n.
78  * There are 7 coefficients for the Shomate polynomial:
79  * - c[0] = \f$A \f$
80  * - c[1] = \f$B \f$
81  * - c[2] = \f$C \f$
82  * - c[3] = \f$D \f$
83  * - c[4] = \f$E \f$
84  * - c[5] = \f$F \f$
85  * - c[6] = \f$G \f$
86  *
87  * See the class description for the polynomial representation of the
88  * thermo functions in terms of \f$A, \dots, G \f$.
89  */
90  ShomatePoly(size_t n, doublereal tlow, doublereal thigh, doublereal pref,
91  const doublereal* coeffs) :
92  m_lowT(tlow),
93  m_highT(thigh),
94  m_Pref(pref),
95  m_index(n) {
96  m_coeff.resize(7);
97  std::copy(coeffs, coeffs + 7, m_coeff.begin());
98  }
99
100  //! copy constructor
101  /*!
102  * @param b object to be copied
103  */
105  m_lowT(b.m_lowT),
106  m_highT(b.m_highT),
107  m_Pref(b.m_Pref),
108  m_coeff(vector_fp(7)),
109  m_index(b.m_index) {
110  std::copy(b.m_coeff.begin(),
111  b.m_coeff.begin() + 7,
112  m_coeff.begin());
113  }
114
115  //! Assignment operator
116  /*!
117  * @param b
118  */
120  if (&b != this) {
121  m_lowT = b.m_lowT;
122  m_highT = b.m_highT;
123  m_Pref = b.m_Pref;
124  m_index = b.m_index;
125  m_coeff.resize(7);
126  std::copy(b.m_coeff.begin(),
127  b.m_coeff.begin() + 7,
128  m_coeff.begin());
129  }
130  return *this;
131  }
132
133  //! Destructor
134  virtual ~ShomatePoly() {}
135
136  //! Duplicator from the base class
137  virtual SpeciesThermoInterpType*
139  ShomatePoly* sp = new ShomatePoly(*this);
140  return (SpeciesThermoInterpType*) sp;
141  }
142
143  //! Returns the minimum temperature that the thermo
144  //! parameterization is valid
145  virtual doublereal minTemp() const {
146  return m_lowT;
147  }
148
149  //! Returns the maximum temperature that the thermo
150  //! parameterization is valid
151  virtual doublereal maxTemp() const {
152  return m_highT;
153  }
154
155  //! Returns the reference pressure (Pa)
156  virtual doublereal refPressure() const {
157  return m_Pref;
158  }
159
160  //! Returns an integer representing the type of parameterization
161  virtual int reportType() const {
162  return SHOMATE;
163  }
164
165  //! Returns an integer representing the species index
166  virtual size_t speciesIndex() const {
167  return m_index;
168  }
169
170
171  //! Update the properties for this species, given a temperature polynomial
172  /*!
173  * This method is called with a pointer to an array containing the functions of
174  * temperature needed by this parameterization, and three pointers to arrays where the
175  * computed property values should be written. This method updates only one value in
176  * each array.
177  *
178  * tt is T/1000.
179  * m_t[0] = tt;
180  * m_t[1] = tt*tt;
181  * m_t[2] = m_t[1]*tt;
182  * m_t[3] = 1.0/m_t[1];
183  * m_t[4] = log(tt);
184  * m_t[5] = 1.0/GasConstant;
185  * m_t[6] = 1.0/(GasConstant * T);
186  *
187  * @param tt Vector of temperature polynomials
188  * @param cp_R Vector of Dimensionless heat capacities.
189  * (length m_kk).
190  * @param h_RT Vector of Dimensionless enthalpies.
191  * (length m_kk).
192  * @param s_R Vector of Dimensionless entropies.
193  * (length m_kk).
194  */
195  virtual void updateProperties(const doublereal* tt,
196  doublereal* cp_R, doublereal* h_RT,
197  doublereal* s_R) const {
198
199  doublereal A = m_coeff[0];
200  doublereal Bt = m_coeff[1]*tt[0];
201  doublereal Ct2 = m_coeff[2]*tt[1];
202  doublereal Dt3 = m_coeff[3]*tt[2];
203  doublereal Etm2 = m_coeff[4]*tt[3];
204  doublereal F = m_coeff[5];
205  doublereal G = m_coeff[6];
206
207  doublereal cp, h, s;
208  cp = A + Bt + Ct2 + Dt3 + Etm2;
209  h = tt[0]*(A + 0.5*Bt + OneThird*Ct2 + 0.25*Dt3 - Etm2) + F;
210  s = A*tt[4] + Bt + 0.5*Ct2 + OneThird*Dt3 - 0.5*Etm2 + G;
211
212  /*
213  * Shomate polynomials parameterizes assuming units of
214  * J/(gmol*K) for cp_r and s_R and kJ/(gmol) for h.
215  * However, Cantera assumes default MKS units of
216  * J/(kmol*K). This requires us to multiply cp and s
217  * by 1.e3 and h by 1.e6, before we then nondimensionlize
218  * the results by dividing by (GasConstant * T),
219  * where GasConstant has units of J/(kmol * K).
220  */
221  cp_R[m_index] = 1.e3 * cp * tt[5];
222  h_RT[m_index] = 1.e6 * h * tt[6];
223  s_R[m_index] = 1.e3 * s * tt[5];
224  }
225
226  //! Compute the reference-state property of one species
227  /*!
228  * Given temperature T in K, this method updates the values of
229  * the non-dimensional heat capacity at constant pressure,
230  * enthalpy, and entropy, at the reference pressure, Pref
231  * of one of the species. The species index is used
232  * to reference into the cp_R, h_RT, and s_R arrays.
233  *
234  * @param temp 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 updatePropertiesTemp(const doublereal temp,
243  doublereal* cp_R, doublereal* h_RT,
244  doublereal* s_R) const {
245  double tPoly[7];
246  doublereal tt = 1.e-3*temp;
247  tPoly[0] = tt;
248  tPoly[1] = tt * tt;
249  tPoly[2] = tPoly[1] * tt;
250  tPoly[3] = 1.0/tPoly[1];
251  tPoly[4] = std::log(tt);
252  tPoly[5] = 1.0/GasConstant;
253  tPoly[6] = 1.0/(GasConstant * temp);
254  updateProperties(tPoly, cp_R, h_RT, s_R);
255  }
256
257  //!This utility function reports back the type of
258  //! parameterization and all of the parameters for the
259  //! species, index.
260  /*!
261  * All parameters are output variables
262  *
263  * @param n Species index
264  * @param type Integer type of the standard type
265  * @param tlow output - Minimum temperature
266  * @param thigh output - Maximum temperature
267  * @param pref output - reference pressure (Pa).
268  * @param coeffs Vector of coefficients used to set the
269  * parameters for the standard state.
270  */
271  virtual void reportParameters(size_t& n, int& type,
272  doublereal& tlow, doublereal& thigh,
273  doublereal& pref,
274  doublereal* const coeffs) const {
275  n = m_index;
276  type = SHOMATE;
277  tlow = m_lowT;
278  thigh = m_highT;
279  pref = m_Pref;
280  for (int i = 0; i < 7; i++) {
281  coeffs[i] = m_coeff[i];
282  }
283  }
284
285  //! Modify parameters for the standard state
286  /*!
287  * @param coeffs Vector of coefficients used to set the
288  * parameters for the standard state.
289  */
290  virtual void modifyParameters(doublereal* coeffs) {
291  if (m_coeff.size() != 7) {
292  throw CanteraError("modifyParameters",
293  "modifying something that hasn't been initialized");
294  }
295  std::copy(coeffs, coeffs + 7, m_coeff.begin());
296  }
297
298 #ifdef H298MODIFY_CAPABILITY
299
300  //! Report the 298 K Heat of Formation of the standard state of one species (J kmol-1)
301  /*!
302  * The 298K Heat of Formation is defined as the enthalpy change to create the standard state
303  * of the species from its constituent elements in their standard states at 298 K and 1 bar.
304  *
305  * @param h298 If this is nonnull, the current value of the Heat of Formation at 298K and 1 bar for
306  * species m_index is returned in h298[m_index].
307  * @return Returns the current value of the Heat of Formation at 298K and 1 bar for
308  * species m_index.
309  */
310  virtual doublereal reportHf298(doublereal* const h298 = 0) const {
311
312  double tPoly[4];
313  doublereal tt = 1.e-3*298.15;
314  tPoly[0] = tt;
315  tPoly[1] = tt * tt;
316  tPoly[2] = tPoly[1] * tt;
317  tPoly[3] = 1.0/tPoly[1];
318
319  doublereal A = m_coeff[0];
320  doublereal Bt = m_coeff[1]*tPoly[0];
321  doublereal Ct2 = m_coeff[2]*tPoly[1];
322  doublereal Dt3 = m_coeff[3]*tPoly[2];
323  doublereal Etm2 = m_coeff[4]*tPoly[3];
324  doublereal F = m_coeff[5];
325
326  doublereal h = tPoly[0]*(A + 0.5*Bt + OneThird*Ct2 + 0.25*Dt3 - Etm2) + F;
327
328  double hh = 1.e6 * h;
329  if (h298) {
330  h298[m_index] = 1.e6 * h;
331  }
332  return hh;
333  }
334
335  //! Modify the value of the 298 K Heat of Formation of one species in the phase (J kmol-1)
336  /*!
337  * The 298K heat of formation is defined as the enthalpy change to create the standard state
338  * of the species from its constituent elements in their standard states at 298 K and 1 bar.
339  *
340  * @param k Species k
341  * @param Hf298New Specify the new value of the Heat of Formation at 298K and 1 bar
342  */
343  virtual void modifyOneHf298(const int k, const doublereal Hf298New) {
344  doublereal hnow = reportHf298();
345  doublereal delH = Hf298New - hnow;
346  m_coeff[5] += delH / 1.0E6;
347  }
348
349 #endif
350
351 protected:
352  //! Minimum temperature for which the parameterization is valid (Kelvin)
353  doublereal m_lowT;
354  //! Maximum temperature for which the parameterization is valid (Kelvin)
355  doublereal m_highT;
356  //! Reference pressure (Pa)
357  doublereal m_Pref;
358  //! Array of coeffcients
360  //! Species Index
361  size_t m_index;
362
363 private:
364
365 };
366
367 //! The Shomate polynomial parameterization for two temperature ranges
368 //! for one species
369 /*!
370  *
371  * Seven coefficients \f$(A,\dots,G)\f$ are used to represent
372  * \f$c_p^0(T)\f$, \f$h^0(T)\f$, and \f$s^0(T) \f$ as
373  * polynomials in the temperature, \f$T \f$, in one temperature region:
374  *
375  * \f[
376  * \tilde{c}_p^0(T) = A + B t + C t^2 + D t^3 + \frac{E}{t^2}
377  * \f]
378  * \f[
379  * \tilde{h}^0(T) = A t + \frac{B t^2}{2} + \frac{C t^3}{3}
380  + \frac{D t^4}{4} - \frac{E}{t} + F.
381  * \f]
382  * \f[
383  * \tilde{s}^0(T) = A\ln t + B t + \frac{C t^2}{2}
384  + \frac{D t^3}{3} - \frac{E}{2t^2} + G.
385  * \f]
386  *
387  * In the above expressions, the thermodynamic polynomials are expressed
388  * in dimensional units, but the temperature,\f$t \f$, is divided by 1000. The
389  * following dimensions are assumed in the above expressions:
390  *
391  * - \f$\tilde{c}_p^0(T)\f$ = Heat Capacity (J/gmol*K)
392  * - \f$\tilde{h}^0(T) \f$ = standard Enthalpy (kJ/gmol)
393  * - \f$\tilde{s}^0(T) \f$= standard Entropy (J/gmol*K)
394  * - \f$t \f$= temperature (K) / 1000.
395  *
397  * http://webbook.nist.gov/
398  *
399  * Before being used within Cantera, the dimensions must be adjusted to those
400  * used by Cantera (i.e., Joules and kmol).
401  *
402  * This function uses two temperature regions, each with a Shomate polynomial
403  * representation to represent the thermo functions. There are 15 coefficients,
404  * therefore, in this representation. The first coefficient is the midrange
405  * temperature.
406  *
407  *
408  * @ingroup spthermo
409  */
411 {
412 public:
413
414  //! Empty constructor
416  : m_lowT(0.0),
417  m_midT(0.0),
418  m_highT(0.0),
419  m_Pref(0.0),
420  msp_low(0),
421  msp_high(0),
422  m_index(0) {
423  m_coeff.resize(15);
424  }
425
426  //! Constructor used in templated instantiations
427  /*!
428  * @param n Species index
429  * @param tlow Minimum temperature
430  * @param thigh Maximum temperature
431  * @param pref reference pressure (Pa).
432  * @param coeffs Vector of coefficients used to set the
433  * parameters for the standard state.
434  * There are 15 coefficients for the 2-zone Shomate polynomial.
435  * The first coefficient is the value of Tmid. The next 7
436  * coefficients are the low temperature range Shomate coefficients.
437  * The last 7 are the high temperature range Shomate coefficients.
438  */
439  ShomatePoly2(size_t n, doublereal tlow, doublereal thigh, doublereal pref,
440  const doublereal* coeffs) :
441  m_lowT(tlow),
442  m_midT(0.0),
443  m_highT(thigh),
444  m_Pref(pref),
445  msp_low(0),
446  msp_high(0),
447  m_index(n) {
448  m_coeff.resize(15);
449  std::copy(coeffs, coeffs + 15, m_coeff.begin());
450  m_midT = coeffs[0];
451  msp_low = new ShomatePoly(n, tlow, m_midT, pref, coeffs+1);
452  msp_high = new ShomatePoly(n, m_midT, thigh, pref, coeffs+8);
453  }
454
455  //! Copy constructor
456  /*!
457  * @param b object to be copied.
458  */
460  m_lowT(b.m_lowT),
461  m_midT(b.m_midT),
462  m_highT(b.m_highT),
463  m_Pref(b.m_Pref),
464  msp_low(0),
465  msp_high(0),
466  m_coeff(vector_fp(15)),
467  m_index(b.m_index) {
468  std::copy(b.m_coeff.begin(),
469  b.m_coeff.begin() + 15,
470  m_coeff.begin());
472  m_Pref, &m_coeff[1]);
474  m_Pref, &m_coeff[8]);
475  }
476
477  //! Assignment operator
478  /*!
479  * @param b object to be copied.
480  */
482  if (&b != this) {
483  m_lowT = b.m_lowT;
484  m_midT = b.m_midT;
485  m_highT = b.m_highT;
486  m_Pref = b.m_Pref;
487  m_index = b.m_index;
488  std::copy(b.m_coeff.begin(),
489  b.m_coeff.begin() + 15,
490  m_coeff.begin());
491  if (msp_low) {
492  delete msp_low;
493  }
494  if (msp_high) {
495  delete msp_high;
496  }
498  m_Pref, &m_coeff[1]);
500  m_Pref, &m_coeff[8]);
501  }
502  return *this;
503  }
504
505  //! Destructor
506  virtual ~ShomatePoly2() {
507  delete msp_low;
508  delete msp_high;
509  }
510
511
512  //! duplicator
513  virtual SpeciesThermoInterpType*
515  ShomatePoly2* sp = new ShomatePoly2(*this);
516  return (SpeciesThermoInterpType*) sp;
517  }
518
519  //! Returns the minimum temperature that the thermo
520  //! parameterization is valid
521  virtual doublereal minTemp() const {
522  return m_lowT;
523  }
524
525  //! Returns the maximum temperature that the thermo
526  //! parameterization is valid
527  virtual doublereal maxTemp() const {
528  return m_highT;
529  }
530
531  //! Returns the reference pressure (Pa)
532  virtual doublereal refPressure() const {
533  return m_Pref;
534  }
535
536  //! Returns an integer representing the type of parameterization
537  virtual int reportType() const {
538  return SHOMATE2;
539  }
540
541  //! Returns an integer representing the species index
542  virtual size_t speciesIndex() const {
543  return m_index;
544  }
545
546  //! Update the properties for this species, given a temperature polynomial
547  /*!
548  * This method is called with a pointer to an array containing the functions of
549  * temperature needed by this parameterization, and three pointers to arrays where the
550  * computed property values should be written. This method updates only one value in
551  * each array.
552  *
553  * Temperature Polynomial:
554  * tt[0] = t;
555  * tt[1] = t*t;
556  * tt[2] = m_t[1]*t;
557  * tt[3] = m_t[2]*t;
558  * tt[4] = 1.0/t;
559  * tt[5] = std::log(t);
560  *
561  * @param tt vector of temperature polynomials
562  * @param cp_R Vector of Dimensionless heat capacities.
563  * (length m_kk).
564  * @param h_RT Vector of Dimensionless enthalpies.
565  * (length m_kk).
566  * @param s_R Vector of Dimensionless entropies.
567  * (length m_kk).
568  */
569  virtual void updateProperties(const doublereal* tt,
570  doublereal* cp_R, doublereal* h_RT,
571  doublereal* s_R) const {
572  double T = 1000 * tt[0];
573  if (T <= m_midT) {
574  msp_low->updateProperties(tt, cp_R, h_RT, s_R);
575  } else {
576  msp_high->updateProperties(tt, cp_R, h_RT, s_R);
577  }
578
579  }
580
581  //! Compute the reference-state property of one species
582  /*!
583  * Given temperature T in K, this method updates the values of
584  * the non-dimensional heat capacity at constant pressure,
585  * enthalpy, and entropy, at the reference pressure, Pref
586  * of one of the species. The species index is used
587  * to reference into the cp_R, h_RT, and s_R arrays.
588  *
589  * @param temp Temperature (Kelvin)
590  * @param cp_R Vector of Dimensionless heat capacities.
591  * (length m_kk).
592  * @param h_RT Vector of Dimensionless enthalpies.
593  * (length m_kk).
594  * @param s_R Vector of Dimensionless entropies.
595  * (length m_kk).
596  */
597  virtual void updatePropertiesTemp(const doublereal temp,
598  doublereal* cp_R,
599  doublereal* h_RT,
600  doublereal* s_R) const {
601  if (temp <= m_midT) {
602  msp_low->updatePropertiesTemp(temp, cp_R, h_RT, s_R);
603  } else {
604  msp_high->updatePropertiesTemp(temp, cp_R, h_RT, s_R);
605  }
606  }
607
608  //!This utility function reports back the type of
609  //! parameterization and all of the parameters for the
610  //! species, index.
611  /*!
612  * All parameters are output variables
613  *
614  * @param n Species index
615  * @param type Integer type of the standard type
616  * @param tlow output - Minimum temperature
617  * @param thigh output - Maximum temperature
618  * @param pref output - reference pressure (Pa).
619  * @param coeffs Vector of coefficients used to set the
620  * parameters for the standard state.
621  */
622  virtual void reportParameters(size_t& n, int& type,
623  doublereal& tlow, doublereal& thigh,
624  doublereal& pref,
625  doublereal* const coeffs) const {
626  n = m_index;
627  type = SHOMATE2;
628  tlow = m_lowT;
629  thigh = m_highT;
630  pref = m_Pref;
631  for (int i = 0; i < 15; i++) {
632  coeffs[i] = m_coeff[i];
633  }
634  }
635
636  //! Modify parameters for the standard state
637  /*!
638  * Here, we take the tact that we will just regenerate the
639  * object.
640  *
641  * @param coeffs Vector of coefficients used to set the
642  * parameters for the standard state.
643  */
644  virtual void modifyParameters(doublereal* coeffs) {
645  delete msp_low;
646  delete msp_high;
647  std::copy(coeffs, coeffs + 15, m_coeff.begin());
648  m_midT = coeffs[0];
649  msp_low = new ShomatePoly(m_index, m_lowT, m_midT, m_Pref, coeffs+1);
650  msp_high = new ShomatePoly(m_index, m_midT, m_highT, m_Pref, coeffs+8);
651  }
652
653 #ifdef H298MODIFY_CAPABILITY
654
655  virtual doublereal reportHf298(doublereal* const h298 = 0) const {
656  doublereal h;
657  if (298.15 <= m_midT) {
658  h = msp_low->reportHf298(h298);
659  } else {
660  h = msp_high->reportHf298(h298);
661  }
662  if (h298) {
663  h298[m_index] = h;
664  }
665  return h;
666  }
667
668  virtual void modifyOneHf298(const int k, const doublereal Hf298New) {
669  if (k != m_index) {
670  return;
671  }
672
673  doublereal h298now = reportHf298(0);
674  doublereal delH = Hf298New - h298now;
675  double h = msp_low->reportHf298(0);
676  double hnew = h + delH;
677  msp_low->modifyOneHf298(k, hnew);
678  h = msp_high->reportHf298(0);
679  hnew = h + delH;
680  msp_high->modifyOneHf298(k, hnew);
681  }
682
683 #endif
684
685 protected:
686  //! Minimum temperature the representation is valid(kelvin)
687  doublereal m_lowT;
688  //! Midrange temperature (kelvin)
689  doublereal m_midT;
690  //! Maximum temperature the representation is valid (kelvin)
691  doublereal m_highT;
692  //! Reference pressure (Pascal)
693  doublereal m_Pref;
694  //! Pointer to the Shomate polynomial for the low temperature region.
696  //! Pointer to the Shomate polynomial for the high temperature region.
698  //! Array of the original coefficients.
700  //! Species index
701  size_t m_index;
702 };
703 }
704
705 #endif