Cantera  2.0
Mu0Poly.cpp
Go to the documentation of this file.
1 /**
2  * @file Mu0Poly.cpp
3  * Definitions for a single-species standard state object derived
4  * from \link Cantera::SpeciesThermoInterpType SpeciesThermoInterpType\endlink based
5  * on a piecewise constant mu0 interpolation
6  * (see \ref spthermo and class \link Cantera::Mu0Poly Mu0Poly\endlink).
7  */
8 #include "cantera/thermo/Mu0Poly.h"
9 #include "cantera/base/ctexceptions.h"
10 #include "cantera/thermo/speciesThermoTypes.h"
11 #include "cantera/thermo/SpeciesThermo.h"
12 #include "cantera/base/xml.h"
13 #include "cantera/base/ctml.h"
14 
15 using namespace std;
16 using namespace ctml;
17 
18 namespace Cantera
19 {
20 
21 
22 Mu0Poly::Mu0Poly() : m_numIntervals(0),
23  m_H298(0.0),
24  m_lowT(0.0),
25  m_highT(0.0),
26  m_Pref(0.0),
27  m_index(0)
28 {
29 }
30 
31 /*
32  * Mu0Poly():
33  *
34  * In the constructor, we calculate and store the
35  * piecewise linear approximation to the thermodynamic
36  * functions.
37  *
38  * coeffs[0] = number of points (integer)
39  * 1 = H298(J/kmol)
40  * 2 = T1 (Kelvin)
41  * 3 = mu1 (J/kmol)
42  * 4 = T2 (Kelvin)
43  * 5 = mu2 (J/kmol)
44  * 6 = T3 (Kelvin)
45  * 7 = mu3 (J/kmol)
46  * ........
47  */
48 Mu0Poly::Mu0Poly(size_t n, doublereal tlow, doublereal thigh,
49  doublereal pref,
50  const doublereal* coeffs) :
51  m_numIntervals(0),
52  m_H298(0.0),
53  m_lowT(tlow),
54  m_highT(thigh),
55  m_Pref(pref),
56  m_index(n)
57 {
58 
59  processCoeffs(coeffs);
60 }
61 
62 
63 Mu0Poly::Mu0Poly(const Mu0Poly& b)
64  : m_numIntervals(b.m_numIntervals),
65  m_H298(b.m_H298),
66  m_t0_int(b.m_t0_int),
67  m_mu0_R_int(b.m_mu0_R_int),
68  m_h0_R_int(b.m_h0_R_int),
69  m_s0_R_int(b.m_s0_R_int),
70  m_cp0_R_int(b.m_cp0_R_int),
71  m_lowT(b.m_lowT),
72  m_highT(b.m_highT),
73  m_Pref(b.m_Pref),
74  m_index(b.m_index)
75 {
76 }
77 
78 Mu0Poly& Mu0Poly::operator=(const Mu0Poly& b)
79 {
80  if (&b != this) {
81  m_numIntervals = b.m_numIntervals;
82  m_H298 = b.m_H298;
83  m_t0_int = b.m_t0_int;
84  m_mu0_R_int = b.m_mu0_R_int;
85  m_h0_R_int = b.m_h0_R_int;
86  m_s0_R_int = b.m_s0_R_int;
87  m_cp0_R_int = b.m_cp0_R_int;
88  m_lowT = b.m_lowT;
89  m_highT = b.m_highT;
90  m_Pref = b.m_Pref;
91  m_index = b.m_index;
92  }
93  return *this;
94 }
95 
96 /*
97  * Destructor:
98  */
99 Mu0Poly::~Mu0Poly()
100 {
101 }
102 
103 SpeciesThermoInterpType*
104 Mu0Poly::duplMyselfAsSpeciesThermoInterpType() const
105 {
106  Mu0Poly* mp = new Mu0Poly(*this);
107  return (SpeciesThermoInterpType*) mp;
108 }
109 
110 doublereal Mu0Poly::minTemp() const
111 {
112  return m_lowT;
113 }
114 doublereal Mu0Poly::maxTemp() const
115 {
116  return m_highT;
117 }
118 doublereal Mu0Poly::refPressure() const
119 {
120  return m_Pref;
121 }
122 
123 /*
124  * updateProperties is the main workhorse program.
125  * Given a temperature (*tt), it calculates the thermodynamic
126  * functions H/RT, S_R, and cp_R, and returns the answer.
127  *
128  * Note, it returns an answer by inserting the values into the
129  * index position, m_index in vectors of H/RT, S_R, and cp_R.
130  *
131  *
132  * Input
133  * -------
134  * *tt = Temperature (Kelvin)
135  *
136  */
137 void Mu0Poly::
138 updateProperties(const doublereal* tt, doublereal* cp_R,
139  doublereal* h_RT, doublereal* s_R) const
140 {
141  size_t j = m_numIntervals;
142  double T = *tt;
143  for (size_t i = 0; i < m_numIntervals; i++) {
144  double T2 = m_t0_int[i+1];
145  if (T <=T2) {
146  j = i;
147  break;
148  }
149  }
150  double T1 = m_t0_int[j];
151  double cp_Rj = m_cp0_R_int[j];
152 
153  doublereal rt = 1.0/T;
154  cp_R[m_index] = cp_Rj;
155  h_RT[m_index] = rt*(m_h0_R_int[j] + (T - T1) * cp_Rj);
156  s_R[m_index] = m_s0_R_int[j] + cp_Rj * (log(T/T1));
157 }
158 
159 void Mu0Poly::
160 updatePropertiesTemp(const doublereal T,
161  doublereal* cp_R,
162  doublereal* h_RT,
163  doublereal* s_R) const
164 {
165  updateProperties(&T, cp_R, h_RT, s_R);
166 }
167 
168 /*
169  * report all of the parameters that make up this
170  * interpolation.
171  *
172  *
173  */
174 void Mu0Poly::reportParameters(size_t& n, int& type,
175  doublereal& tlow, doublereal& thigh,
176  doublereal& pref,
177  doublereal* const coeffs) const
178 {
179  n = m_index;
180  type = MU0_INTERP;
181  tlow = m_lowT;
182  thigh = m_highT;
183  pref = m_Pref;
184  coeffs[0] = int(m_numIntervals)+1;
185  coeffs[1] = m_H298 * GasConstant;
186  int j = 2;
187  for (size_t i = 0; i < m_numIntervals+1; i++) {
188  coeffs[j] = m_t0_int[i];
189  coeffs[j+1] = m_mu0_R_int[i] * GasConstant;
190  j += 2;
191  }
192 }
193 
194 void Mu0Poly::modifyParameters(doublereal* coeffs)
195 {
196  processCoeffs(coeffs);
197 }
198 
199 /*
200  * Install a Mu0 polynomial thermodynamic reference state property
201  * parameterization for species k into a SpeciesThermo instance,
202  * getting the information from an XML database.
203  */
204 void installMu0ThermoFromXML(std::string speciesName,
205  SpeciesThermo& sp, size_t k,
206  const XML_Node* Mu0Node_ptr)
207 {
208 
209  doublereal tmin, tmax;
210  bool dimensionlessMu0Values = false;
211  const XML_Node& Mu0Node = *Mu0Node_ptr;
212 
213  tmin = fpValue(Mu0Node["Tmin"]);
214  tmax = fpValue(Mu0Node["Tmax"]);
215  doublereal pref = fpValue(Mu0Node["Pref"]);
216 
217  doublereal h298 = 0.0;
218  if (Mu0Node.hasChild("H298")) {
219  h298 = getFloat(Mu0Node, "H298", "actEnergy");
220  }
221 
222  size_t numPoints = 1;
223  if (Mu0Node.hasChild("numPoints")) {
224  numPoints = getInteger(Mu0Node, "numPoints");
225  }
226 
227  vector_fp cValues(numPoints);
228  const XML_Node* valNode_ptr =
229  getByTitle(const_cast<XML_Node&>(Mu0Node), "Mu0Values");
230  if (!valNode_ptr) {
231  throw CanteraError("installMu0ThermoFromXML",
232  "missing required while processing "
233  + speciesName);
234  }
235  getFloatArray(*valNode_ptr, cValues, true, "actEnergy");
236  /*
237  * Check to see whether the Mu0's were input in a dimensionless
238  * form. If they were, then the assumed temperature needs to be
239  * adjusted from the assumed T = 273.15
240  */
241  string uuu = (*valNode_ptr)["units"];
242  if (uuu == "Dimensionless") {
243  dimensionlessMu0Values = true;
244  }
245  size_t ns = cValues.size();
246  if (ns != numPoints) {
247  throw CanteraError("installMu0ThermoFromXML",
248  "numPoints inconsistent while processing "
249  + speciesName);
250  }
251 
252  vector_fp cTemperatures(numPoints);
253  const XML_Node* tempNode_ptr =
254  getByTitle(const_cast<XML_Node&>(Mu0Node), "Mu0Temperatures");
255  if (!tempNode_ptr) {
256  throw CanteraError("installMu0ThermoFromXML",
257  "missing required while processing + "
258  + speciesName);
259  }
260  getFloatArray(*tempNode_ptr, cTemperatures, false);
261  ns = cTemperatures.size();
262  if (ns != numPoints) {
263  throw CanteraError("installMu0ThermoFromXML",
264  "numPoints inconsistent while processing "
265  + speciesName);
266  }
267 
268  /*
269  * Fix up dimensionless Mu0 values if input
270  */
271  if (dimensionlessMu0Values) {
272  for (size_t i = 0; i < numPoints; i++) {
273  cValues[i] *= cTemperatures[i] / 273.15;
274  }
275  }
276 
277 
278  vector_fp c(2 + 2 * numPoints);
279 
280  c[0] = static_cast<double>(numPoints);
281  c[1] = h298;
282  for (size_t i = 0; i < numPoints; i++) {
283  c[2+i*2] = cTemperatures[i];
284  c[2+i*2+1] = cValues[i];
285  }
286 
287  sp.install(speciesName, k, MU0_INTERP, &c[0], tmin, tmax, pref);
288 }
289 
290 /*
291  * Mu0Poly():
292  *
293  * In the constructor, we calculate and store the
294  * piecewise linear approximation to the thermodynamic
295  * functions.
296  *
297  * coeffs[0] = number of points (integer)
298  * 1 = H298(J/kmol)
299  * 2 = T1 (Kelvin)
300  * 3 = mu1 (J/kmol)
301  * 4 = T2 (Kelvin)
302  * 5 = mu2 (J/kmol)
303  * 6 = T3 (Kelvin)
304  * 7 = mu3 (J/kmol)
305  * ........
306  */
307 void Mu0Poly::processCoeffs(const doublereal* coeffs)
308 {
309 
310  size_t i, iindex;
311  double T1, T2;
312  size_t nPoints = (size_t) coeffs[0];
313  if (nPoints < 2) {
314  throw CanteraError("Mu0Poly",
315  "nPoints must be >= 2");
316  }
317  m_numIntervals = nPoints - 1;
318  m_H298 = coeffs[1] / GasConstant;
319  size_t iT298 = 0;
320  /*
321  * Resize according to the number of points
322  */
323  m_t0_int.resize(nPoints);
324  m_h0_R_int.resize(nPoints);
325  m_s0_R_int.resize(nPoints);
326  m_cp0_R_int.resize(nPoints);
327  m_mu0_R_int.resize(nPoints);
328  /*
329  * Calculate the T298 interval and make sure that
330  * the temperatures are strictly monotonic.
331  * Also distribute the data into the internal arrays.
332  */
333  bool ifound = false;
334  for (i = 0, iindex = 2; i < nPoints; i++) {
335  T1 = coeffs[iindex];
336  m_t0_int[i] = T1;
337  m_mu0_R_int[i] = coeffs[iindex+1] / GasConstant;
338  if (T1 == 298.15) {
339  iT298 = i;
340  ifound = true;
341  }
342  if (i < nPoints - 1) {
343  T2 = coeffs[iindex+2];
344  if (T2 <= T1) {
345  throw CanteraError("Mu0Poly",
346  "Temperatures are not monotonic increasing");
347  }
348  }
349  iindex += 2;
350  }
351  if (!ifound) {
352  throw CanteraError("Mu0Poly",
353  "One temperature has to be 298.15");
354  }
355 
356  /*
357  * Starting from the interval with T298, we go up
358  */
359  doublereal mu2, s1, s2, h1, h2, cpi, deltaMu, deltaT;
360  T1 = m_t0_int[iT298];
361  doublereal mu1 = m_mu0_R_int[iT298];
362  m_h0_R_int[iT298] = m_H298;
363  m_s0_R_int[iT298] = - (mu1 - m_h0_R_int[iT298]) / T1;
364  for (i = iT298; i < m_numIntervals; i++) {
365  T1 = m_t0_int[i];
366  s1 = m_s0_R_int[i];
367  h1 = m_h0_R_int[i];
368  mu1 = m_mu0_R_int[i];
369  T2 = m_t0_int[i+1];
370  mu2 = m_mu0_R_int[i+1];
371  deltaMu = mu2 - mu1;
372  deltaT = T2 - T1;
373  cpi = (deltaMu - T1 * s1 + T2 * s1) / (deltaT - T2 * log(T2/T1));
374  h2 = h1 + cpi * deltaT;
375  s2 = s1 + cpi * log(T2/T1);
376  m_cp0_R_int[i] = cpi;
377  m_h0_R_int[i+1] = h2;
378  m_s0_R_int[i+1] = s2;
379  m_cp0_R_int[i+1] = cpi;
380  }
381 
382  /*
383  * Starting from the interval with T298, we go down
384  */
385  if (iT298 != 0) {
386  T2 = m_t0_int[iT298];
387  mu2 = m_mu0_R_int[iT298];
388  m_h0_R_int[iT298] = m_H298;
389  m_s0_R_int[iT298] = - (mu2 - m_h0_R_int[iT298]) / T2;
390  for (i = iT298 - 1; i != npos; i--) {
391  T1 = m_t0_int[i];
392  mu1 = m_mu0_R_int[i];
393  T2 = m_t0_int[i+1];
394  mu2 = m_mu0_R_int[i+1];
395  s2 = m_s0_R_int[i+1];
396  h2 = m_h0_R_int[i+1];
397  deltaMu = mu2 - mu1;
398  deltaT = T2 - T1;
399  cpi = (deltaMu - T1 * s2 + T2 * s2) / (deltaT - T1 * log(T2/T1));
400  h1 = h2 - cpi * deltaT;
401  s1 = s2 - cpi * log(T2/T1);
402  m_cp0_R_int[i] = cpi;
403  m_h0_R_int[i] = h1;
404  m_s0_R_int[i] = s1;
405  if (i == (m_numIntervals-1)) {
406  m_cp0_R_int[i+1] = cpi;
407  }
408  }
409  }
410 #ifdef DEBUG_HKM_NOT
411  printf(" Temp mu0(J/kmol) cp0(J/kmol/K) "
412  " h0(J/kmol) s0(J/kmol/K) \n");
413  for (i = 0; i < nPoints; i++) {
414  printf("%12.3g %12.5g %12.5g %12.5g %12.5g\n",
415  m_t0_int[i], m_mu0_R_int[i] * GasConstant,
416  m_cp0_R_int[i]* GasConstant,
417  m_h0_R_int[i]* GasConstant,
418  m_s0_R_int[i]* GasConstant);
419  fflush(stdout);
420  }
421 #endif
422 }
423 
424 }
425 
426 
427 
428 
429 
Generated by   doxygen 1.8.2