Cantera  2.0
thermoFunctions.cpp
Go to the documentation of this file.
1 /**
2  * @file thermoFunctions.cpp
3  * File containing thermo evaluation functions for NASA polynomials,
4  * which are used in testing the interpolations.
5  */
6
7 // Copyright 2001 California Institute of Technology
8
9 #include <math.h>
10 #include "thermoFunctions.h"
11 #include <iostream>
12 using namespace std;
13
14 namespace ckr
15 {
16
17 /**
18  * non-dimensional heat capacity (\f$C_p/R \f$) at constant P for
19  * one species @param t temperature @param s species object
20  */
21 double cp(double t, const Species& s)
22 {
23  if (s.thermoFormatType == 1) {
24  const vector_fp* cpc;
25  int ireg = -1;
26  for (int i = 0; i < s.nTempRegions; i++) {
27  if (t <= s.maxTemps[i]) {
28  ireg = i;
29  break;
30  }
31  }
32  cpc = s.region_coeffs[ireg];
33  const vector_fp& c = *cpc;
34  double cp0r = c[0]/(t*t) + c[1]/t + c[2] + c[3]*t + c[4]*t*t
35  + c[5]*t*t*t + c[6]*t*t*t*t;
36  return cp0r;
37  }
38  const vector_fp* cpc;
39  if (t > s.tmid) {
40  cpc = &s.highCoeffs;
41  } else {
42  cpc = &s.lowCoeffs;
43  }
44  const vector_fp& c = *cpc;
45  double cp0r = c[0] + c[1]*t + c[2]*t*t + c[3]*t*t*t + c[4]*t*t*t*t;
46  return cp0r;
47 }
48
49
50
51 /**
52  * enthalpy in Kelvin (\f$H/R \f$) for
53  * one species. @param t temperature @param s species object
54  */
55 double enthalpy(double t, const Species& s)
56 {
57  if (s.thermoFormatType == 1) {
58  const vector_fp* cpc;
59  int ireg = -1;
60  for (int i = 0; i < s.nTempRegions; i++) {
61  if (t <= s.maxTemps[i]) {
62  ireg = i;
63  break;
64  }
65  }
66  cpc = s.region_coeffs[ireg];
67  const vector_fp& c = *cpc;
68  double h0rt = -c[0]/(t*t) + c[1]*log(t)/t
69  + c[2] + 0.5*c[3]*t + c[4]*t*t/3.0 + 0.25*c[5]*t*t*t
70  + 0.2*c[6]*t*t*t*t + c[7]/t;
71  return t*h0rt;
72  }
73  const vector_fp* cp;
74  if (t > s.tmid) {
75  cp = &s.highCoeffs;
76  } else {
77  cp = &s.lowCoeffs;
78  }
79  const vector_fp& c = *cp;
80  double h0rt = c[0] + 0.5*c[1]*t + c[2]*t*t/3.0 + 0.25*c[3]*t*t*t
81  + 0.2*c[4]*t*t*t*t + c[5]/t;
82  return t*h0rt;
83 }
84
85
86 /**
87  * non-dimensional entropy (\f$S/R \f$) for
88  * one species @param t temperature @param s species object
89  */
90 double entropy(double t, const Species& s)
91 {
92  if (s.thermoFormatType == 1) {
93  const vector_fp* cpc;
94  int ireg = -1;
95  for (int i = 0; i < s.nTempRegions; i++) {
96  if (t <= s.maxTemps[i]) {
97  ireg = i;
98  break;
99  }
100  }
101  cpc = s.region_coeffs[ireg];
102  const vector_fp& c = *cpc;
103  double s0r = -0.5*c[0]/(t*t) - c[1]/t
104  + c[2]*log(t) + c[3]*t + 0.5*c[4]*t*t + c[5]*t*t*t/3.0
105  + 0.25*c[6]*t*t*t*t + c[8];
106  return t*s0r;
107  }
108  const vector_fp* cp;
109  if (t > s.tmid) {
110  cp = &s.highCoeffs;
111  } else {
112  cp = &s.lowCoeffs;
113  }
114  const vector_fp& c = *cp;
115  double s0r = c[0]*log(t) + c[1]*t + 0.5*c[2]*t*t + c[3]*t*t*t/3.0
116  + 0.25*c[4]*t*t*t*t + c[6];
117  return t*s0r;
118 }
119
120 /**
121  * Gibbs function in Kelvin (\f$G/R \f$) for
122  * one species. @param t temperature @param s species object
123  */
124 double gibbs(double t, const Species& s)
125 {
126  if (s.thermoFormatType == 1) {
127  double s0r = entropy(t, s);
128  double h0r = enthalpy(t, s);
129  return (h0r - s0r * t);
130  }
131  const vector_fp* cp;
132  if (t > s.tmid) {
133  cp = &s.highCoeffs;
134  } else {
135  cp = &s.lowCoeffs;
136  }
137  const vector_fp& c = *cp;
138  double h0rt = c[0] + 0.5*c[1]*t + c[2]*t*t/3.0 + 0.25*c[3]*t*t*t
139  + 0.2*c[4]*t*t*t*t + c[5]/t;
140  double s0r = c[0]*log(t) + c[1]*t + 0.5*c[2]*t*t + c[3]*t*t*t/3.0
141  + 0.25*c[4]*t*t*t*t + c[6];
142  return t*(h0rt - s0r);
143 }
144
145 }
146