Cantera  2.0
CarbonDioxide.cpp
1 /* FILE: CarbonDioxide.cpp
2  * DESCRIPTION:
3  * representation of substance Carbon Dioxide
4  * values and functions are from
5  * "Thermodynamic Properties in SI" bu W.C. Reynolds
6  * AUTHOR: me@rebeccahhunt.com: GCEP, Stanford University
7  *
8  */
9
10 #include "CarbonDioxide.h"
11 #include <math.h>
12 #include <string.h>
13
14 namespace tpx
15 {
16
17 /*
18  * Carbon Dioxide constants
19  */
20 static const double Tmn = 216.54; // [K] minimum temperature for which calculations are valid
21 static const double Tmx = 1500.0; // [K] maximum temperature for which calculations are valid
22 static const double Tc=304.21; // [K] critical temperature
23 static const double Roc=464.00; // [kg/m^3] critical density
24 static const double To=216.54; // [K] reference Temperature
25 static const double R=188.918; // [] gas constant for CO2 J/kg/K
26 static const double Gamma=5.0E-6; // [??]
27 static const double u0=3.2174105E5; // [] internal energy at To
28 static const double s0=2.1396056E3; // [] entropy at To
29 static const double Tp=250; // [K] ??
30 static const double Pc=7.38350E6; // [Pa] critical pressure
31 static const double M=44.01; // [kg/kmol] molar density
32
33 /*
34  * array Acarbdi is used by the function named Pp
35  */
36 static const double Acarbdi[]= {
37  2.2488558E-1,
38  -1.3717965E2,
39  -1.4430214E4,
40  -2.9630491E6,
41  -2.0606039E8,
42  4.5554393E-5,
43  7.7042840E-2,
44  4.0602371E1,
45  4.0029509E-7,
46  -3.9436077E-4,
47  1.2115286E-10,
48  1.0783386E-7,
49  4.3962336E-11,
50  -3.6505545E4,
51  1.9490511E7,
52  -2.9186718E9,
53  2.4358627E-2,
54  -3.7546530E1,
55  1.1898141E4
56 };
57
58
59 /*
60  * array F is used by the function named Psat
61  */
62 static const double F[]= {
63  -6.5412610,
64  -2.7914636E-1,
65  -3.4716202,
66  -3.4989637,
67  -1.9770948E1,
68  1.3922839E2,
69  -2.7670389E2,
70  -7.0510251E3
71 };
72
73
74 /*
75  * array D is used by the function ldens
76  */
77 static const double D[]= {
78  4.6400009E2,
79  6.7938129E2,
80  1.4776836E3,
81  -3.1267676E3,
82  3.6397656E3,
83  -1.3437098E3
84 };
85
86
87 /*
88  * array G is used by the function sp
89  */
90 static const double G[]= {
91  8.726361E3,
92  1.840040E2,
93  1.914025,
94  -1.667825E-3,
95  7.305950E-7,
96  -1.255290E-10,
97 };
98
99 /*
100  * C returns a multiplier in each term of the sum
101  * in P-3, used in conjunction with C in the function Pp
102  * j is used to represent which of the values in the summation to calculate
103  * j=0 is the second additive in the formula in reynolds
104  * j=1 is the third...
105  * (this part does not include the multiplier rho^n)
106  */
107 double CarbonDioxide::C(int j,double Tinverse, double T2inverse, double T3inverse, double T4inverse)
108 {
109  switch (j) {
110  case 0 :
111  return Acarbdi[0]*T +
112  Acarbdi[1] +
113  Acarbdi[2] * Tinverse +
114  Acarbdi[3] * T2inverse +
115  Acarbdi[4] * T3inverse ;
116  case 1 :
117  return Acarbdi[5] *T +
118  Acarbdi[6] +
119  Acarbdi[7] * Tinverse ;
120  case 2 :
121  return Acarbdi[8]*T + Acarbdi[9];
122  case 3 :
123  return Acarbdi[10]*T + Acarbdi[11];
124  case 4 :
125  return Acarbdi[12];
126  case 5 :
127  return Acarbdi[13] *T2inverse +
128  Acarbdi[14] *T3inverse +
129  Acarbdi[15] *T4inverse;
130  case 6 :
131  return Acarbdi[16] *T2inverse +
132  Acarbdi[17] *T3inverse +
133  Acarbdi[18] *T4inverse;
134  default :
135  return 0.0;
136  }
137 }
138
139 /* cprime
140  * derivative of C(i)
141  */
142 inline double CarbonDioxide::Cprime(int j, double T2inverse, double T3inverse, double T4inverse)
143 {
144
145  switch (j) {
146  case 0 :
147  return Acarbdi[0] +
148  - Acarbdi[2] * T2inverse +
149  -2 * Acarbdi[3] * T3inverse +
150  -3 * Acarbdi[4] * T4inverse ;
151  case 1 :
152  return Acarbdi[5] -
153  Acarbdi[7] * T2inverse;
154  case 2 :
155  return Acarbdi[8] ;
156  case 3 :
157  return Acarbdi[10] ;
158  case 4 :
159  return 0;
160  case 5 :
161  return
162  -2 *Acarbdi[13] *T3inverse +
163  -3 *Acarbdi[14] *T4inverse +
164  -4 *Acarbdi[15]* pow(T,-5);
165  case 6 :
166  return
167  -2 *Acarbdi[16] *T3inverse +
168  -3 *Acarbdi[17] *T4inverse +
169  -4 *Acarbdi[18] *pow(T,-5);
170  default :
171  return 0.0;
172  }
173 }
174
175 /*
176  * I = integral from o-rho { 1/(rho^2) * H(i, rho) d rho }
177  * ( see section 2 of Reynolds TPSI )
178  */
179 inline double CarbonDioxide::I(int j, double ergho, double Gamma)
180 {
181  switch (j) {
182
183  case 0:
184  return Rho;
185  case 1:
186  return pow(Rho, 2)/2;
187  case 2:
188  return pow(Rho, 3)/ 3;
189  case 3:
190  return pow(Rho, 4)/ 4;
191  case 4:
192  return pow(Rho, 5)/ 5;
193  case 5:
194  return (1 - ergho) / double(2 * Gamma);
195  case 6:
196  return (1 - ergho * double(Gamma * pow(Rho,2) + double(1)))/ double(2 * Gamma * Gamma);
197  default:
198  return 0.0;
199  }
200 }
201
202
203 /* H returns a multiplier in each term of the sum
204  * in P-3
205  * this is used in conjunction with C in the function Pp
206  * this represents the product rho^n
207  * i=0 is the second additive in the formula in reynolds
208  * i=1 is the third ...
209  */
210 double CarbonDioxide::H(int i, double egrho)
211 {
212  if (i < 5) {
213  return pow(Rho,i+2);
214  } else if (i == 5) {
215  return pow(Rho,3)*egrho;
216  } else if (i == 6) {
217  return pow(Rho,5)*egrho;
218  } else {
219  return 0;
220  }
221 }
222
223 /*
224  * internal energy
225  * see Reynolds eqn (15) section 2
226  * u = (the integral from T to To of co(T)dT) +
227  * sum from i to N ([C(i) - T*Cprime(i)] + uo
228  */
229 double CarbonDioxide::up()
230 {
231
232  double Tinverse = 1.0/T;
233  double T2inverse = pow(T, -2);
234  double T3inverse = pow(T, -3);
235  double T4inverse = pow(T, -4);
236  double egrho = exp(-Gamma*Rho*Rho);
237
238  double sum = 0.0;
239
240  // Equation C-6 integrated
241  sum += G[0]*log(T/To);
242  int i;
243  for (i=1; i<=5; i++) {
244  sum += G[i]*(pow(T,i) - pow(To,i))/double(i);
245  }
246
247
248  for (i=0; i<=6; i++) {
249  sum += I(i,egrho, Gamma) *
250  (C(i, Tinverse, T2inverse, T3inverse, T4inverse) - T*Cprime(i,T2inverse, T3inverse, T4inverse));
251  }
252
253  sum += u0;
254  return sum + m_energy_offset;
255
256 }
257
258 /*
259 * entropy
260  * see Reynolds eqn (16) section 2
261 */
262
263 double CarbonDioxide::sp()
264 {
265  //double Tinverse = 1.0/T;
266  double T2inverse = pow(T, -2);
267  double T3inverse = pow(T, -3);
268  double T4inverse = pow(T, -4);
269  double egrho = exp(-Gamma*Rho*Rho);
270
271  double sum = 0.0;
272
273  for (int i=2; i<=5; i++) {
274  sum += G[i]*(pow(T,i-1) - pow(To,i-1))/double(i-1);
275  }
276
277  sum += G[1]*log(T/To);
278  sum -= G[0]*(1.0/T - 1.0/To);
279
280
281  for (int i=0; i<=6; i++) {
282  sum -= Cprime(i,T2inverse, T3inverse, T4inverse)*I(i,egrho,Gamma);
283  }
284
285  sum += s0 - R*log(Rho);
286
287  return sum + m_entropy_offset;
288 }
289
290
291 /*
292  * Equation P-3 in Reynolds
293  * P - rho - T
294  * returns P (pressure)
295  */
296 double CarbonDioxide::Pp()
297 {
298  double Tinverse = pow(T,-1);
299  double T2inverse = pow(T, -2);
300  double T3inverse = pow(T, -3);
301  double T4inverse = pow(T, -4);
302  double egrho = exp(-Gamma*Rho*Rho);
303
304  double P = Rho*R*T;
305
306  // when i=0 we are on second sum of equation (where rho^2)
307  for (int i=0; i<=6; i++) {
308  P += C(i,Tinverse, T2inverse, T3inverse, T4inverse)*H(i,egrho);
309  }
310  return P;
311 }
312
313
314 /*
315  * Equation S-2 in Reynolds
316  * Pressure at Saturation
317  */
318 double CarbonDioxide::Psat()
319 {
320
321  double log, sum=0,P;
322  if ((T < Tmn) || (T > Tc)) {
323  std::cout << " error in Psat " << TempError << std::endl;
324  set_Err(TempError); // Error("CarbonDioxide::Psat",TempError,T);
325  }
326  for (int i=1; i<=8; i++) {
327  sum += F[i-1] * pow((T/Tp -1),double(i-1));
328  }
329
330  log = ((Tc/T)-1)*sum;
331  P=exp(log)*Pc;
332
333  //cout << "Psat is returning " << P << " at T " << T << " and Pc " << Pc << " and Tp " << Tp << endl;
334  return P;
335
336 }
337
338 /*
339  * Equation D2 in Reynolds
340  * liquid density, of rho_f
341  */
342 double CarbonDioxide::ldens()
343 {
344  double xx=1-(T/Tc), sum=0;
345  if ((T < Tmn) || (T > Tc)) {
346  std::cout << " error in ldens " << TempError << std::endl;
347  set_Err(TempError);
348  }
349  for (int i=1; i<=6; i++) {
350  sum+=D[i-1]*pow(xx,double(i-1)/3.0);
351  }
352
353  return sum;
354 }
355
356 /*
357  * the following functions allow users
358  * to get the properties of CarbonDioxide
359  * that are not dependent on the state
360  */
361 double CarbonDioxide::Tcrit()
362 {
363  return Tc;
364 }
365 double CarbonDioxide::Pcrit()
366 {
367  return Pc;
368 }
369 double CarbonDioxide::Vcrit()
370 {
371  return 1.0/Roc;
372 }
373 double CarbonDioxide::Tmin()
374 {
375  return Tmn;
376 }
377 double CarbonDioxide::Tmax()
378 {
379  return Tmx;
380 }
381 char* CarbonDioxide::name()
382 {
383  return (char*) m_name.c_str() ;
384 }
385 char* CarbonDioxide::formula()
386 {
387  return (char*) m_formula.c_str();
388 }
389 double CarbonDioxide::MolWt()
390 {
391  return M;
392 }
393
394 }
395
396
397