Cantera  2.0
Nasa9PolyMultiTempRegion.cpp
Go to the documentation of this file.
1 /**
2  * @file Nasa9PolyMultiTempRegion.cpp
3  * Definitions for a single-species standard state object derived
6  * on the NASA 9 coefficient temperature polynomial form
7  * applied to one temperature region
8  * (see \ref spthermo and class
10  *
11  * This parameterization has one NASA temperature region.
12  */
13 // Copyright 2007 Sandia National Laboratories
14
15 #include "cantera/base/global.h"
18
19 using namespace std;
20
21 namespace Cantera
22 {
23
24 // The NASA 9 polynomial parameterization for a single species
25 // encompassing multiple temperature regions.
26 /*
27  * This parameterization expresses the heat capacity via a
28  * 7 coefficient polynomial.
29  * Note that this is the form used in the
30  * 2002 NASA equilibrium program. A reference to the form is
31  * provided below:
32  *
33  * "NASA Glenn Coefficients for Calculating Thermodynamic
34  * Properties of Individual Species,"
35  * B. J. McBride, M. J. Zehe, S. Gordon
36  * NASA/TP-2002-211556, Sept. 2002
37  *
38  * Nine coefficients \f$(a_0,\dots,a_6)\f$ are used to represent
39  * \f$C_p^0(T)\f$, \f$H^0(T)\f$, and \f$S^0(T) \f$ as
40  * polynomials in \f$T \f$ :
41  * \f[
42  * \frac{c_p(T)}{R} = a_0 T^{-2} + a_1 T^{-1} + a_2 + a_3 T
43  * + a_4 T^2 + a_5 T^3 + a_6 T^4
44  * \f]
45  *
46  * \f[
47  * \frac{H^0(T)}{RT} = - a_0 T^{-2} + a_1 \frac{\ln(T)}{T} + a_2
48  * + a_3 T + a_4 T^2 + a_5 T^3 + a_6 T^4 + \frac{a_7}{T}
49  * \f]
50  *
51  * \f[
52  * \frac{s^0(T)}{R} = - \frac{a_0}{2} T^{-2} - a_1 T^{-1} + a_2 \ln(T)
53  + + a_3 T \frac{a_4}{2} T^2 + \frac{a_5}{3} T^3
54  * + \frac{a_6}{4} T^4 + a_8
55  * \f]
56  *
57  * The standard state is assumed to be an ideal gas at the
58  * standard pressure of 1 bar, for gases.
59  * For condensed species, the standard state is the
60  * pure crystalline or liquid substance at the standard
61  * pressure of 1 atm.
62  *
63  * These NASA representations may have multiple temperature regions
64  * through the use of this %Nasa9PolyMultiTempRegion object, which uses
65  * multiple copies of the Nasa9Poly1 object to handle multiple temperature
66  * regions.
67  *
68  * @ingroup spthermo
69  */
70
71
72 //! Empty constructor
73 Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion() :
74  m_lowT(0.0),
75  m_highT(0.0),
76  m_Pref(0.0),
77  m_index(0),
78  m_numTempRegions(0),
79  m_currRegion(0)
80 {
81 }
82
83
84 // Constructor used in templated instantiations
85 /*
86  * @param regionPts Vector of pointers to Nasa9Poly1 objects. These
87  * objects all refer to the temperature regions for the
88  * same species. The vector must be in increasing
89  * temperature region format. Together they
90  * represent the reference temperature parameterization
91  * for a single species.
92  */
94 Nasa9PolyMultiTempRegion(std::vector<Cantera::Nasa9Poly1*> &regionPts) :
95  m_lowT(0.0),
96  m_highT(0.0),
97  m_Pref(0.0),
98  m_index(0),
99  m_numTempRegions(0),
100  m_currRegion(0)
101 {
102  m_numTempRegions = regionPts.size();
103  // Do a shallow copy of the pointers. From now on, we will
104  // own these pointers and be responsible for deleting them.
105  m_regionPts = regionPts;
107  m_lowT = m_regionPts[0]->minTemp();
108  m_highT = m_regionPts[m_numTempRegions-1]->maxTemp();
109  m_Pref = m_regionPts[0]->refPressure();
110  m_index = m_regionPts[0]->speciesIndex();
111  for (size_t i = 0; i < m_numTempRegions; i++) {
112  m_lowerTempBounds[i] = m_regionPts[i]->minTemp();
113  if (m_regionPts[i]->speciesIndex() != m_index) {
114  throw CanteraError("Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion",
115  "m_index inconsistency");
116  }
117  if (fabs(m_regionPts[i]->refPressure() - m_Pref) > 0.0001) {
118  throw CanteraError("Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion",
119  "refPressure inconsistency");
120  }
121  if (i > 0) {
122  if (m_lowerTempBounds[i-1] >= m_lowerTempBounds[i]) {
123  throw CanteraError("Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion",
124  "minTemp bounds inconsistency");
125  }
126  if (fabs(m_regionPts[i-1]->maxTemp() - m_lowerTempBounds[i]) > 0.0001) {
127  throw CanteraError("Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion",
128  "Temp bounds inconsistency");
129  }
130  }
131  }
132 }
133
134 // copy constructor
135 /*
136  * @param b object to be copied
137  */
140  m_lowT(b.m_lowT),
141  m_highT(b.m_highT),
142  m_Pref(b.m_Pref),
143  m_index(b.m_index),
144  m_numTempRegions(b.m_numTempRegions),
145  m_lowerTempBounds(b.m_lowerTempBounds),
146  m_currRegion(b.m_currRegion)
147 {
149  for (size_t i = 0; i < m_numTempRegions; i++) {
150  Nasa9Poly1* dptr = b.m_regionPts[i];
151  m_regionPts[i] = new Nasa9Poly1(*dptr);
152  }
153 }
154
155 // assignment operator
156 /*
157  * @param b object to be copied
158  */
161 {
162  if (&b != this) {
163  for (size_t i = 0; i < m_numTempRegions; i++) {
164  delete m_regionPts[i];
165  m_regionPts[i] = 0;
166  }
167  m_lowT = b.m_lowT;
168  m_highT = b.m_highT;
169  m_Pref = b.m_Pref;
170  m_index = b.m_index;
171  m_numTempRegions = b.m_numTempRegions;
174  m_regionPts.resize(m_numTempRegions);
175  for (size_t i = 0; i < m_numTempRegions; i++) {
176  m_regionPts[i] = new Nasa9Poly1(*(b.m_regionPts[i]));
177  }
178  }
179  return *this;
180 }
181
182 // Destructor
184 {
185  for (size_t i = 0; i < m_numTempRegions; i++) {
186  delete m_regionPts[i];
187  m_regionPts[i] = 0;
188  }
189 }
190
191 // duplicator
194 {
196  return (SpeciesThermoInterpType*) np;
197 }
198
199 // Returns the minimum temperature that the thermo
200 // parameterization is valid
202 {
203  return m_lowT;
204 }
205
206 // Returns the maximum temperature that the thermo
207 // parameterization is valid
209 {
210  return m_highT;
211 }
212
213 // Returns the reference pressure (Pa)
215 {
216  return m_Pref;
217 }
218
219 // Returns an integer representing the type of parameterization
221 {
222  return NASA9MULTITEMP;
223 }
224
225
226 // Returns an integer representing the species index
228 {
229  return m_index;
230 }
231
232
233 // Update the properties for this species, given a temperature polynomial
234 /*
235  * This method is called with a pointer to an array containing the functions of
236  * temperature needed by this parameterization, and three pointers to arrays where the
237  * computed property values should be written. This method updates only one value in
238  * each array.
239  *
240  * Temperature Polynomial:
241  * tt[0] = t;
242  * tt[1] = t*t;
243  * tt[2] = t*t*t;
244  * tt[3] = t*t*t*t;
245  * tt[4] = 1.0/t;
246  * tt[5] = 1.0/(t*t);
247  * tt[6] = std::log(t);
248  *
249  * @param tt vector of temperature polynomials
250  * @param cp_R Vector of Dimensionless heat capacities.
251  * (length m_kk).
252  * @param h_RT Vector of Dimensionless enthalpies.
253  * (length m_kk).
254  * @param s_R Vector of Dimensionless entropies.
255  * (length m_kk).
256  */
258  doublereal* cp_R,
259  doublereal* h_RT,
260  doublereal* s_R) const
261 {
262  // Let's put some additional debugging here.
263  // This is an external routine
264 #ifdef DEBUG_HKM
265  double temp = tt[0];
266  if (temp < m_regionPts[m_currRegion]->minTemp()) {
267  if (m_currRegion != 0) {
268  throw CanteraError("Nasa9PolyMultiTempRegion::updateProperties",
269  "region problem");
270  }
271  }
272  if (temp > m_regionPts[m_currRegion]->maxTemp()) {
273  if (m_currRegion != m_numTempRegions - 1) {
274  throw CanteraError("Nasa9PolyMultiTempRegion::updateProperties",
275  "region problem");
276  }
277  }
278 #endif
279  (m_regionPts[m_currRegion])->updateProperties(tt, cp_R, h_RT, s_R);
280 }
281
282
283 // Compute the reference-state property of one species
284 /*
285  * Given temperature T in K, this method updates the values of
286  * the non-dimensional heat capacity at constant pressure,
287  * enthalpy, and entropy, at the reference pressure, Pref
288  * of one of the species. The species index is used
289  * to reference into the cp_R, h_RT, and s_R arrays.
290  *
291  * Temperature Polynomial:
292  * tt[0] = t;
293  * tt[1] = t*t;
294  * tt[2] = t*t*t;
295  * tt[3] = t*t*t*t;
296  * tt[4] = 1.0/t;
297  * tt[5] = 1.0/(t*t);
298  * tt[6] = std::log(t);
299  *
300  * @param temp Temperature (Kelvin)
301  * @param cp_R Vector of Dimensionless heat capacities.
302  * (length m_kk).
303  * @param h_RT Vector of Dimensionless enthalpies.
304  * (length m_kk).
305  * @param s_R Vector of Dimensionless entropies.
306  * (length m_kk).
307  */
309  doublereal* cp_R, doublereal* h_RT,
310  doublereal* s_R) const
311 {
312  double tPoly[7];
313  tPoly[0] = temp;
314  tPoly[1] = temp * temp;
315  tPoly[2] = tPoly[1] * temp;
316  tPoly[3] = tPoly[2] * temp;
317  tPoly[4] = 1.0 / temp;
318  tPoly[5] = tPoly[4] / temp;
319  tPoly[6] = std::log(temp);
320  // Now find the region
321  m_currRegion = 0;
322  for (size_t i = 1; i < m_numTempRegions; i++) {
323  if (temp < m_lowerTempBounds[i]) {
324  break;
325  }
326  m_currRegion++;
327  }
328
329  updateProperties(tPoly, cp_R, h_RT, s_R);
330 }
331
332 //This utility function reports back the type of
333 // parameterization and all of the parameters for the
334 // species, index.
335 /*
336  * All parameters are output variables
337  *
338  * @param n Species index
339  * @param type Integer type of the standard type
340  * @param tlow output - Minimum temperature
341  * @param thigh output - Maximum temperature
342  * @param pref output - reference pressure (Pa).
343  * @param coeffs Vector of coefficients used to set the
344  * parameters for the standard state.
345  */
347  doublereal& tlow, doublereal& thigh,
348  doublereal& pref,
349  doublereal* const coeffs) const
350 {
351  n = m_index;
352  type = NASA9MULTITEMP;
353  tlow = m_lowT;
354  thigh = m_highT;
355  pref = m_Pref;
356  double ctmp[12];
357  coeffs[0] = double(m_numTempRegions);
358  int index = 1;
359  size_t n_tmp = 0;
360  int type_tmp = 0;
361  double pref_tmp = 0.0;
362  for (size_t iReg = 0; iReg < m_numTempRegions; iReg++) {
363  m_regionPts[iReg]->reportParameters(n_tmp, type_tmp,
364  coeffs[index], coeffs[index+1],
365  pref_tmp, ctmp);
366  for (int i = 0; i < 9; i++) {
367  coeffs[index+2+i] = ctmp[3+i];
368  }
369  index += 11;
370  }
371
372 }
373
374 // Modify parameters for the standard state
375 /*
376  * @param coeffs Vector of coefficients used to set the
377  * parameters for the standard state.
378  */
380 {
381  int index = 3;
382  for (size_t iReg = 0; iReg < m_numTempRegions; iReg++) {
383  m_regionPts[iReg]->modifyParameters(coeffs + index);
384  index += 11;
385  }
386 }
387
388 }