Cantera  2.0
Hydrogen.cpp
1 // Hydrogen
2 
3 #include "Hydrogen.h"
4 #include <math.h>
5 
6 namespace tpx
7 {
8 
9 static const double
10 M = 2.0159,
11 Tmn = 13.8,
12 Tmx = 5000.0,
13 Tc = 32.938,
14 Pc = 1.2838e6,
15 Roc= 31.36,
16 To = 13.8,
17 Tt = 13.8,
18 Pt = 7042.09,
19 R = 4124.299539,
20 Gamma = 1.008854772e-3,
21 u0 = 3.9275114e5,
22 s0 = 2.3900333e4,
23 T1 = 35,
24 T2 = 400,
25 alpha = 1.5814454428, //to be used with psat
26 alpha1 = .3479; //to be used with ldens
27 
28 static const double Ahydro[] = {
29  1.150470519352900e1, 1.055427998826072e3, -1.270685949968568e4,
30  7.287844527295619e4, -7.448780703363973e5, 2.328994151810363e-1,
31  -1.635308393739296e1, 3.730678064960389e3, 6.299667723184813e5,
32  1.210920358305697e-3, 1.753651095884817, -1.367022988058101e2,
33  -6.869936641299885e-3, 3.644494201750974e-2, -2.559784772600182,
34  -4.038855202905836e-4, 1.485396303520942e-6, 4.243613981060742e-4,
35  -2.307910113586888e-6, -6.082192173879582e5, -1.961080967486886e6,
36  -5.786932854076408e2, 2.799129504191752e4, -2.381566558300913e-1,
37  8.918796032452872e-1, -6.985739539036644e-5, -7.339554179182899e-3,
38  -5.597033440289980e-9, 8.842130160884514e-8, -2.655507264539047e-12,
39  -4.544474518140164e-12, 9.818775257001922e-11
40 };
41 static const double Dhydro[]= {
42  4.8645813003e1, -3.4779278180e1, 4.0776538192e2,
43  -1.1719787304e3, 1.6213924400e3, -1.1531096683e3, 3.3825492039e2
44 };
45 static const double Fhydro[]=
46 { 3.05300134164, 2.80810925813, -6.55461216567e-1, 1.59514439374 };
47 static const double Ghydro[]= {
48  6.1934792e3, 2.9490437e2, -1.5401979e3, -4.9176101e3,
49  6.8957165e4, -2.2282185e5, 3.7990059e5, -3.7094216e5, 2.1326792e5,
50  -7.1519411e4, 1.2971743e4, -9.8533014e2, 1.0434776e4,
51  -3.9144179e2, 5.8277696e2, 6.5409163e2, -1.8728847e2
52 };
53 
54 
55 double hydrogen::C(int i, double rt, double rt2)
56 {
57  switch (i) {
58  case 0 :
59  return Ahydro[0] * T + Ahydro[1] * sqrt(T) + Ahydro[2] + (Ahydro[3] + Ahydro[4] * rt) * rt;
60  case 1 :
61  return Ahydro[5] * T + Ahydro[6] + rt * (Ahydro[7] + Ahydro[8] * rt);
62  case 2 :
63  return Ahydro[9] * T + Ahydro[10] + Ahydro[11] * rt;
64  case 3 :
65  return Ahydro[12];
66  case 4 :
67  return rt*(Ahydro[13] + Ahydro[14]*rt);
68  case 5 :
69  return Ahydro[15]*rt;
70  case 6 :
71  return rt*(Ahydro[16] + Ahydro[17]*rt);
72  case 7 :
73  return Ahydro[18]*rt2;
74  case 8 :
75  return rt2*(Ahydro[19] + Ahydro[20]*rt);
76  case 9 :
77  return rt2*(Ahydro[21] + Ahydro[22]*rt2);
78  case 10 :
79  return rt2*(Ahydro[23] + Ahydro[24]*rt);
80  case 11 :
81  return rt2*(Ahydro[25] + Ahydro[26]*rt2);
82  case 12 :
83  return rt2*(Ahydro[27] + Ahydro[28]*rt);
84  case 13 :
85  return rt2*(Ahydro[29] + Ahydro[30]*rt + Ahydro[31]*rt2);
86  default :
87  return 0.0;
88  }
89 }
90 
91 double hydrogen::Cprime(int i, double rt, double rt2, double rt3)
92 {
93  switch (i) {
94  case 0 :
95  return Ahydro[0] + 0.5*Ahydro[1]/sqrt(T) - (Ahydro[3] + 2.0*Ahydro[4]*rt)*rt2;
96  case 1 :
97  return Ahydro[5] - rt2*(Ahydro[7] + 2.0*Ahydro[8]*rt);
98  case 2 :
99  return Ahydro[9] - Ahydro[11]*rt2;
100  case 3 :
101  return 0.0;
102  case 4 :
103  return -rt2*(Ahydro[13] + 2.0*Ahydro[14]*rt);
104  case 5 :
105  return -Ahydro[15]*rt2;
106  case 6 :
107  return -rt2*(Ahydro[16] + 2.0*Ahydro[17]*rt);
108  case 7 :
109  return -2.0*Ahydro[18]*rt3;
110  case 8 :
111  return -rt3*(2.0*Ahydro[19] + 3.0*Ahydro[20]*rt);
112  case 9 :
113  return -rt3*(2.0*Ahydro[21] + 4.0*Ahydro[22]*rt2);
114  case 10 :
115  return -rt3*(2.0*Ahydro[23] + 3.0*Ahydro[24]*rt);
116  case 11 :
117  return -rt3*(2.0*Ahydro[25] + 4.0*Ahydro[26]*rt2);
118  case 12 :
119  return -rt3*(2.0*Ahydro[27] + 3.0*Ahydro[28]*rt);
120  case 13 :
121  return -rt3*(2.0*Ahydro[29] + 3.0*Ahydro[30]*rt + 4.0*Ahydro[31]*rt2);
122  default :
123  return 0.0;
124  }
125 }
126 
127 double hydrogen::W(int n, double egrho)
128 {
129  return (n == 0 ? (1.0 - egrho)/(2.0*Gamma) :
130  (n*W(n-1, egrho) - 0.5*pow(Rho,2*n)*egrho)/Gamma);
131 }
132 
133 double hydrogen::H(int i, double egrho)
134 {
135  return (i < 8 ? pow(Rho,i+2) : pow(Rho,2*i-13)*egrho);
136 }
137 
138 double hydrogen::I(int i, double egrho)
139 {
140  return (i < 8 ? pow(Rho,i+1)/double(i+1) : W(i-8, egrho));
141 }
142 
143 double hydrogen::icv(int i, double x, double xlg)
144 {
145  return (i == 0 ? x - 1 : x*pow(xlg,i) - i*icv(i-1,x,xlg));
146 }
147 
148 double hydrogen::up()
149 {
150  double rt = 1.0/T;
151  double rt2 = rt*rt;
152  double rt3 = rt*rt2;
153  double egrho = exp(-Gamma*Rho*Rho);
154  double sum = u0;
155  for (int i=0; i<14; i++) {
156  sum += (C(i, rt, rt2) - T*Cprime(i, rt, rt2, rt3))*I(i, egrho);
157  }
158 
159  // add \int c_{v,0} term
160  sum += Ghydro[0] * (std::min(T, T1) - To);
161  if (T > T1) {
162  double x = std::min(T, T2) / T1;
163  for (int i = 0; i < 12; i++) {
164  sum += Ghydro[i] * T1 * icv(i, x, log(x));
165  }
166  }
167  if (T > T2) {
168  double x = T/T2;
169  for (int i = 0; i < 5; i++) {
170  sum += Ghydro[i+12] * T2 * icv(i, x, log(x));
171  }
172  }
173 
174  return sum + m_energy_offset;
175 }
176 
177 double hydrogen::sp()
178 {
179  double rt = 1.0/T;
180  double rt2 = rt*rt;
181  double rt3 = rt*rt2;
182  double egrho = exp(-Gamma*Rho*Rho);
183  double sum = s0 - R*log(Rho);
184  for (int i=0; i<14; i++) {
185  sum -= Cprime(i, rt, rt2, rt3)*I(i, egrho);
186  }
187 
188  // add \int c_{v,0}/T term
189  sum += Ghydro[0] * log(std::min(T, T1)/ To);
190  if (T > T1) {
191  double xlg = log(std::min(T, T2)/T1);
192  for (int i = 0; i < 12; i++) {
193  sum += Ghydro[i] / (i + 1) * pow(xlg, i+1);
194  }
195  }
196  if (T > T2) {
197  double xlg = log(T/T2);
198  for (int i = 0; i < 5; i++) {
199  sum += Ghydro[i+12] / (i + 1) * pow(xlg, i+1);
200  }
201  }
202 
203  return sum + m_entropy_offset;
204 }
205 
206 double hydrogen::Pp()
207 {
208  double rt = 1.0/T;
209  double rt2 = rt*rt;
210  // double rt3 = rt*rt2;
211  double egrho = exp(-Gamma*Rho*Rho);
212 
213  double P = Rho*R*T;
214  for (int i=0; i<14; i++) {
215  P += C(i, rt, rt2)*H(i, egrho);
216  }
217  return P;
218 }
219 
220 //equation D4
221 double hydrogen::ldens()
222 {
223  if ((T < Tmn) || (T > Tc)) {
224  set_Err(TempError);
225  }
226  double x=1-T/Tc;
227  double sum, term;
228  int i;
229  for (i=1, sum=0; i<=6; i++) {
230  sum+=Dhydro[i]*pow(x, 1+double(i-1)/3.0);
231  }
232  term = sum+Roc+Dhydro[0]*pow(x,alpha1);
233  return term;
234 }
235 
236 
237 //equation s3
238 double hydrogen::Psat()
239 {
240  double x = (1.0 - Tt/T)/(1.0 - Tt/Tc);
241  double result;
242  if ((T < Tmn) || (T > Tc)) {
243  set_Err(TempError);
244  }
245  result = Fhydro[0]*x + Fhydro[1]*x*x + Fhydro[2]*x*x*x +
246  Fhydro[3]*x*pow(1-x, alpha);
247  return exp(result)*Pt;
248 }
249 
250 double hydrogen::Tcrit()
251 {
252  return Tc;
253 }
254 double hydrogen::Pcrit()
255 {
256  return Pc;
257 }
258 double hydrogen::Vcrit()
259 {
260  return 1.0/Roc;
261 }
262 double hydrogen::Tmin()
263 {
264  return Tmn;
265 }
266 double hydrogen::Tmax()
267 {
268  return Tmx;
269 }
270 char* hydrogen::name()
271 {
272  return (char*) m_name.c_str();
273 }
274 char* hydrogen::formula()
275 {
276  return (char*) m_formula.c_str();
277 }
278 double hydrogen::MolWt()
279 {
280  return M;
281 }
282 
283 }
Generated by   doxygen 1.8.2