Cantera  2.1.2
Hydrogen.cpp
Go to the documentation of this file.
1 //! @file Hydrogen.cpp
2 #include "Hydrogen.h"
4 
5 using namespace Cantera;
6 
7 namespace tpx
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 double hydrogen::C(int i, double rt, double rt2)
55 {
56  switch (i) {
57  case 0 :
58  return Ahydro[0] * T + Ahydro[1] * sqrt(T) + Ahydro[2] + (Ahydro[3] + Ahydro[4] * rt) * rt;
59  case 1 :
60  return Ahydro[5] * T + Ahydro[6] + rt * (Ahydro[7] + Ahydro[8] * rt);
61  case 2 :
62  return Ahydro[9] * T + Ahydro[10] + Ahydro[11] * rt;
63  case 3 :
64  return Ahydro[12];
65  case 4 :
66  return rt*(Ahydro[13] + Ahydro[14]*rt);
67  case 5 :
68  return Ahydro[15]*rt;
69  case 6 :
70  return rt*(Ahydro[16] + Ahydro[17]*rt);
71  case 7 :
72  return Ahydro[18]*rt2;
73  case 8 :
74  return rt2*(Ahydro[19] + Ahydro[20]*rt);
75  case 9 :
76  return rt2*(Ahydro[21] + Ahydro[22]*rt2);
77  case 10 :
78  return rt2*(Ahydro[23] + Ahydro[24]*rt);
79  case 11 :
80  return rt2*(Ahydro[25] + Ahydro[26]*rt2);
81  case 12 :
82  return rt2*(Ahydro[27] + Ahydro[28]*rt);
83  case 13 :
84  return rt2*(Ahydro[29] + Ahydro[30]*rt + Ahydro[31]*rt2);
85  default :
86  return 0.0;
87  }
88 }
89 
90 double hydrogen::Cprime(int i, double rt, double rt2, double rt3)
91 {
92  switch (i) {
93  case 0 :
94  return Ahydro[0] + 0.5*Ahydro[1]/sqrt(T) - (Ahydro[3] + 2.0*Ahydro[4]*rt)*rt2;
95  case 1 :
96  return Ahydro[5] - rt2*(Ahydro[7] + 2.0*Ahydro[8]*rt);
97  case 2 :
98  return Ahydro[9] - Ahydro[11]*rt2;
99  case 3 :
100  return 0.0;
101  case 4 :
102  return -rt2*(Ahydro[13] + 2.0*Ahydro[14]*rt);
103  case 5 :
104  return -Ahydro[15]*rt2;
105  case 6 :
106  return -rt2*(Ahydro[16] + 2.0*Ahydro[17]*rt);
107  case 7 :
108  return -2.0*Ahydro[18]*rt3;
109  case 8 :
110  return -rt3*(2.0*Ahydro[19] + 3.0*Ahydro[20]*rt);
111  case 9 :
112  return -rt3*(2.0*Ahydro[21] + 4.0*Ahydro[22]*rt2);
113  case 10 :
114  return -rt3*(2.0*Ahydro[23] + 3.0*Ahydro[24]*rt);
115  case 11 :
116  return -rt3*(2.0*Ahydro[25] + 4.0*Ahydro[26]*rt2);
117  case 12 :
118  return -rt3*(2.0*Ahydro[27] + 3.0*Ahydro[28]*rt);
119  case 13 :
120  return -rt3*(2.0*Ahydro[29] + 3.0*Ahydro[30]*rt + 4.0*Ahydro[31]*rt2);
121  default :
122  return 0.0;
123  }
124 }
125 
126 double hydrogen::W(int n, double egrho)
127 {
128  return (n == 0 ? (1.0 - egrho)/(2.0*Gamma) :
129  (n*W(n-1, egrho) - 0.5*pow(Rho,2*n)*egrho)/Gamma);
130 }
131 
132 double hydrogen::H(int i, double egrho)
133 {
134  return (i < 8 ? pow(Rho,i+2) : pow(Rho,2*i-13)*egrho);
135 }
136 
137 double hydrogen::I(int i, double egrho)
138 {
139  return (i < 8 ? pow(Rho,i+1)/double(i+1) : W(i-8, egrho));
140 }
141 
142 double hydrogen::icv(int i, double x, double xlg)
143 {
144  return (i == 0 ? x - 1 : x*pow(xlg,i) - i*icv(i-1,x,xlg));
145 }
146 
147 double hydrogen::up()
148 {
149  double rt = 1.0/T;
150  double rt2 = rt*rt;
151  double rt3 = rt*rt2;
152  double egrho = exp(-Gamma*Rho*Rho);
153  double sum = u0;
154  for (int i=0; i<14; i++) {
155  sum += (C(i, rt, rt2) - T*Cprime(i, rt, rt2, rt3))*I(i, egrho);
156  }
157 
158  // add \int c_{v,0} term
159  sum += Ghydro[0] * (std::min(T, T1) - To);
160  if (T > T1) {
161  double x = std::min(T, T2) / T1;
162  for (int i = 0; i < 12; i++) {
163  sum += Ghydro[i] * T1 * icv(i, x, log(x));
164  }
165  }
166  if (T > T2) {
167  double x = T/T2;
168  for (int i = 0; i < 5; i++) {
169  sum += Ghydro[i+12] * T2 * icv(i, x, log(x));
170  }
171  }
172 
173  return sum + m_energy_offset;
174 }
175 
176 double hydrogen::sp()
177 {
178  double rt = 1.0/T;
179  double rt2 = rt*rt;
180  double rt3 = rt*rt2;
181  double egrho = exp(-Gamma*Rho*Rho);
182  double sum = s0 - R*log(Rho);
183  for (int i=0; i<14; i++) {
184  sum -= Cprime(i, rt, rt2, rt3)*I(i, egrho);
185  }
186 
187  // add \int c_{v,0}/T term
188  sum += Ghydro[0] * log(std::min(T, T1)/ To);
189  if (T > T1) {
190  double xlg = log(std::min(T, T2)/T1);
191  for (int i = 0; i < 12; i++) {
192  sum += Ghydro[i] / (i + 1) * pow(xlg, i+1);
193  }
194  }
195  if (T > T2) {
196  double xlg = log(T/T2);
197  for (int i = 0; i < 5; i++) {
198  sum += Ghydro[i+12] / (i + 1) * pow(xlg, i+1);
199  }
200  }
201 
202  return sum + m_entropy_offset;
203 }
204 
205 double hydrogen::Pp()
206 {
207  double rt = 1.0/T;
208  double rt2 = rt*rt;
209  // double rt3 = rt*rt2;
210  double egrho = exp(-Gamma*Rho*Rho);
211 
212  double P = Rho*R*T;
213  for (int i=0; i<14; i++) {
214  P += C(i, rt, rt2)*H(i, egrho);
215  }
216  return P;
217 }
218 
219 double hydrogen::ldens()
220 {
221  if ((T < Tmn) || (T > Tc)) {
222  throw TPX_Error("hydrogen::ldens",
223  "Temperature out of range. T = " + fp2str(T));
224  }
225  double x=1-T/Tc;
226  double sum;
227  int i;
228  for (i=1, sum=0; i<=6; i++) {
229  sum+=Dhydro[i]*pow(x, 1+double(i-1)/3.0);
230  }
231  return sum+Roc+Dhydro[0]*pow(x,alpha1);
232 }
233 
234 double hydrogen::Psat()
235 {
236  double x = (1.0 - Tt/T)/(1.0 - Tt/Tc);
237  double result;
238  if ((T < Tmn) || (T > Tc)) {
239  throw TPX_Error("hydrogen::Psat",
240  "Temperature out of range. T = " + fp2str(T));
241  }
242  result = Fhydro[0]*x + Fhydro[1]*x*x + Fhydro[2]*x*x*x +
243  Fhydro[3]*x*pow(1-x, alpha);
244  return exp(result)*Pt;
245 }
246 
247 double hydrogen::Tcrit()
248 {
249  return Tc;
250 }
251 double hydrogen::Pcrit()
252 {
253  return Pc;
254 }
255 double hydrogen::Vcrit()
256 {
257  return 1.0/Roc;
258 }
259 double hydrogen::Tmin()
260 {
261  return Tmn;
262 }
263 double hydrogen::Tmax()
264 {
265  return Tmx;
266 }
267 char* hydrogen::name()
268 {
269  return (char*) m_name.c_str();
270 }
271 char* hydrogen::formula()
272 {
273  return (char*) m_formula.c_str();
274 }
275 double hydrogen::MolWt()
276 {
277  return M;
278 }
279 
280 }
std::string fp2str(const double x, const std::string &fmt)
Convert a double into a c++ string.
Definition: stringUtils.cpp:29
Contains declarations for string manipulation functions within Cantera.