Cantera  2.0
Heptane.cpp
1 /* FILE: Heptane.cpp
2  * DESCRIPTION:
3  * representation of substance Heptane
4  * values and functions are from
5  * "Thermodynamic Properties in SI" bu W.C. Reynolds
6  * AUTHOR: jrh@stanford.edu: GCEP, Stanford University
7  *
8  */
9
10 #include "Heptane.h"
11 #include <math.h>
12 #include <string.h>
13
14 namespace tpx
15 {
16
17 /*
18  * Heptane constants
19  */
20 static const double Tmn = 182.56; // [K] minimum temperature for which calculations are valid
21 static const double Tmx = 1000.0; // [K] maximum temperature for which calculations are valid
22 static const double Tc=537.68; // [K] critical temperature
23 static const double Roc=197.60; // [kg/m^3] critical density
24 static const double To=300; // [K] reference Temperature
25 static const double R=82.99504; // [J/(kg*K)] gas constant (for this substance)
26 static const double Gamma=9.611604E-6; // [??]
27 static const double u0=3.4058439E5; // [] internal energy at To
28 static const double s0=1.1080254E3; // [] entropy at To
29 static const double Tp=400; // [K] ??
30 static const double Pc=2.6199E6; // [Pa] critical pressure
31 static const double M=100.20; // [kg/kmol] molar density
32
33 /*
34  * array Ahept is used by the function Pp
35  */
36 static const double Ahept[]= {
37  2.246032E-3,
38  2.082990E2,
39  5.085746E7,
40  3.566396E9,
41  1.622168E9,
42  1.065237E-5,
43  5.987922E-1,
44  7.736602,
45  1.929386E5,
46  5.291379E-9
47 };
48
49
50 /*
51  * array F is used by Psat
52  */
53 static const double F[]= {
54  -7.2298764,
55  3.8607475E-1,
56  -3.4216472,
57  4.6274432E-1,
58  -9.7926124,
59  -4.2058094E1,
60  7.5468678E1,
61  3.1758992E2
62 };
63
64
65 /*
66  * array D is used by the function ldens
67  */
68 static const double D[]= {
69  1.9760405E2,
70  8.9451237E2,
71  -1.1462908E3,
72  1.7996947E3,
73  -1.7250843E3,
74  9.7088329E2
75 };
76
77
78 /*
79  * array G is used by the function sp
80  */
81 static const double G[]= {
82  1.1925213E5,
83  -7.7231363E2,
84  7.4463527,
85  -3.0888167E-3,
86  0.0,
87  0.0
88 };
89
90
91 /*
92  * C returns a multiplier in each term of the sum
93  * in P-2, used in conjunction with C in the function Pp
94  * j is used to represent which of the values in the summation to calculate
95  * j=0 is the second additive in the formula in reynolds
96  * j=1 is the third...
97  */
98 double Heptane::C(int j,double Tinverse, double T2inverse, double T3inverse, double T4inverse)
99 {
100  switch (j) {
101  case 0 :
102  return Ahept[0] * R * T -
103  Ahept[1] -
104  Ahept[2] * T2inverse +
105  Ahept[3] * T3inverse -
106  Ahept[4] * T4inverse;
107  case 1 :
108  return Ahept[5] * R * T -
109  Ahept[6] -
110  Ahept[7] * Tinverse;
111  case 2 :
112  return Ahept[9] * (Ahept[6] + Ahept[7] * Tinverse);
113  case 3 :
114  return Ahept[8] * T2inverse;
115  default :
116  return 0.0;
117  }
118 }
119
120
121 /* cprime
122  * derivative of C(i)
123  */
124 inline double Heptane::Cprime(int j, double T2inverse, double T3inverse, double T4inverse)
125 {
126  switch (j) {
127  case 0 :
128  return Ahept[0] * R -
129  -2 * Ahept[2] * T3inverse +
130  -3 * Ahept[3] * T4inverse -
131  -4 * Ahept[4] * pow(T, -5.0);
132  case 1 :
133  return Ahept[5] * R -
134  -1 * Ahept[7] * T2inverse;
135  case 2 :
136  return Ahept[9] * (-1 * Ahept[7] * T2inverse);
137  case 3 :
138  return -2 * Ahept[8] * T3inverse;
139  default :
140  return 0.0;
141  }
142 }
143
144
145 /*
146  * I = integral from o-rho { 1/(rho^2) * H(i, rho) d rho }
147  * ( see section 2 of Reynolds TPSI )
148  */
149 inline double Heptane::I(int j, double ergho, double Gamma)
150 {
151  switch (j) {
152  case 0:
153  return Rho;
154  case 1:
155  return Rho * Rho / 2;
156  case 2:
157  return pow(Rho, 5.0)/ 5;
158  case 3:
159  return 1 / Gamma - (Gamma * Rho * Rho + 2) * ergho / (2 * Gamma);
160  default:
161  return 0.0;
162  }
163 }
164
165
166 /* H returns a multiplier in each term of the sum
167  * in P-2
168  * this is used in conjunction with C in the function Pp
169  * this represents the product rho^n
170  * i=0 is the second additive in the formula in reynolds
171  * i=1 is the third ...
172  */
173 double Heptane::H(int i, double egrho)
174 {
175  if (i < 2) {
176  return pow(Rho,i+2);
177  } else if (i == 2) {
178  return pow(Rho,6.0);
179  } else if (i == 3) {
180  return pow(Rho,3) * (1 + Gamma * Rho * Rho) * egrho;
181  } else {
182  return 0;
183  }
184 }
185
186
187 /*
188  * internal energy
189  * see Reynolds eqn (15) section 2
190  * u = (the integral from T to To of co(T)dT) +
191  * sum from i to N ([C(i) - T*Cprime(i)] + uo
192  */
193 double Heptane::up()
194 {
195  double Tinverse = 1.0/T;
196  double T2inverse = pow(T, -2);
197  double T3inverse = pow(T, -3);
198  double T4inverse = pow(T, -4);
199  double egrho = exp(-Gamma*Rho*Rho);
200
201  double sum = 0.0;
202  int i;
203  for (i=1; i<=5; i++) {
204  sum += G[i]*(pow(T,i) - pow(To,i))/double(i);
205  }
206
207  sum += G[0]*log(T/To);
208
209  for (i=0; i<=6; i++) {
210  sum += (C(i, Tinverse, T2inverse, T3inverse, T4inverse) - T*Cprime(i,T2inverse, T3inverse, T4inverse))*I(i,egrho, Gamma);
211  }
212
213  sum += u0;
214
215  return sum + m_energy_offset;
216 }
217
218
219 /*
220  * entropy
221  * see Reynolds eqn (16) section 2
222  */
223 double Heptane::sp()
224 {
225  double T2inverse = pow(T, -2);
226  double T3inverse = pow(T, -3);
227  double T4inverse = pow(T, -4);
228  double egrho = exp(-Gamma*Rho*Rho);
229
230  double sum = 0.0;
231
232  for (int i=2; i<=5; i++) {
233  sum += G[i]*(pow(T,i-1) - pow(To,i-1))/double(i-1);
234  }
235
236  sum += G[1]*log(T/To);
237  sum -= G[0]*(1.0/T - 1.0/To);
238
239  for (int i=0; i<=6; i++) {
240  sum -= Cprime(i,T2inverse, T3inverse, T4inverse)*I(i,egrho, Gamma);
241  }
242
243  sum += s0 - R*log(Rho);
244
245  return sum + m_entropy_offset;
246 }
247
248
249 /*
250  * Equation P-2 in Reynolds
251  * P - rho - T
252  * returns P (pressure)
253  */
254 double Heptane::Pp()
255 {
256  double Tinverse = pow(T,-1);
257  double T2inverse = pow(T, -2);
258  double T3inverse = pow(T, -3);
259  double T4inverse = pow(T, -4);
260  double egrho = exp(-Gamma*Rho*Rho);
261
262  double P = Rho*R*T;
263
264  for (int i=0; i<=3; i++) {
265  P += C(i,Tinverse, T2inverse, T3inverse, T4inverse)*H(i,egrho);
266  }
267
268  return P;
269 }
270
271
272 /*
273  * Equation S-2 in Reynolds
274  * Pressure at Saturation
275  */
276 double Heptane::Psat()
277 {
278  double log, sum=0,P;
279  if ((T < Tmn) || (T > Tc)) {
280
281  set_Err(TempError); // Error("Heptane::Psat",TempError,T);
282  }
283  for (int i=1; i<=8; i++) {
284  sum += F[i-1] * pow((T/Tp -1),double(i-1));
285  }
286
287  log = ((Tc/T)-1)*sum;
288  P=exp(log)*Pc;
289
290  return P;
291 }
292
293
294 /*
295  * Equation D2 in Reynolds
296  * liquid density, of rho_f
297  */
298 double Heptane::ldens()
299 {
300  double xx=1-(T/Tc), sum=0;
301  if ((T < Tmn) || (T > Tc)) {
302  set_Err(TempError);
303  }
304  for (int i=1; i<=6; i++) {
305  sum+=D[i-1]*pow(xx,double(i-1)/3.0);
306  }
307
308  return sum;
309 }
310
311
312 /*
313  * the following functions allow users
314  * to get the properties of Heptane
315  * that are not dependent on the state
316  */
317 double Heptane::Tcrit()
318 {
319  return Tc;
320 }
321 double Heptane::Pcrit()
322 {
323  return Pc;
324 }
325 double Heptane::Vcrit()
326 {
327  return 1.0/Roc;
328 }
329 double Heptane::Tmin()
330 {
331  return Tmn;
332 }
333 double Heptane::Tmax()
334 {
335  return Tmx;
336 }
337 char* Heptane::name()
338 {
339  return (char*) m_name.c_str();
340 }
341 char* Heptane::formula()
342 {
343  return (char*) m_formula.c_str();
344 }
345 double Heptane::MolWt()
346 {
347  return M;
348 }
349 }