Cantera  3.1.0a1
Loading...
Searching...
No Matches
Hydrogen.cpp
Go to the documentation of this file.
1//! @file Hydrogen.cpp
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
6#include "Hydrogen.h"
8
9using namespace Cantera;
10
11namespace tpx
12{
13static const double
14M = 2.0159,
15Tmn = 13.8,
16Tmx = 5000.0,
17Tc = 32.938,
18Pc = 1.2838e6,
19Roc= 31.36,
20To = 13.8,
21Tt = 13.8,
22Pt = 7042.09,
23R = 4124.299539,
24Gamma = 1.008854772e-3,
25u0 = 3.9275114e5,
26s0 = 2.3900333e4,
27T1 = 35,
28T2 = 400,
29alpha = 1.5814454428, //to be used with psat
30alpha1 = .3479; //to be used with ldens
31
32static const double Ahydro[] = {
33 1.150470519352900e1, 1.055427998826072e3, -1.270685949968568e4,
34 7.287844527295619e4, -7.448780703363973e5, 2.328994151810363e-1,
35 -1.635308393739296e1, 3.730678064960389e3, 6.299667723184813e5,
36 1.210920358305697e-3, 1.753651095884817, -1.367022988058101e2,
37 -6.869936641299885e-3, 3.644494201750974e-2, -2.559784772600182,
38 -4.038855202905836e-4, 1.485396303520942e-6, 4.243613981060742e-4,
39 -2.307910113586888e-6, -6.082192173879582e5, -1.961080967486886e6,
40 -5.786932854076408e2, 2.799129504191752e4, -2.381566558300913e-1,
41 8.918796032452872e-1, -6.985739539036644e-5, -7.339554179182899e-3,
42 -5.597033440289980e-9, 8.842130160884514e-8, -2.655507264539047e-12,
43 -4.544474518140164e-12, 9.818775257001922e-11
44};
45static const double Dhydro[]= {
46 4.8645813003e1, -3.4779278180e1, 4.0776538192e2,
47 -1.1719787304e3, 1.6213924400e3, -1.1531096683e3, 3.3825492039e2
48};
49static const double Fhydro[]=
50{ 3.05300134164, 2.80810925813, -6.55461216567e-1, 1.59514439374 };
51static const double Ghydro[]= {
52 6.1934792e3, 2.9490437e2, -1.5401979e3, -4.9176101e3,
53 6.8957165e4, -2.2282185e5, 3.7990059e5, -3.7094216e5, 2.1326792e5,
54 -7.1519411e4, 1.2971743e4, -9.8533014e2, 1.0434776e4,
55 -3.9144179e2, 5.8277696e2, 6.5409163e2, -1.8728847e2
56};
57
58double hydrogen::C(int i, double rt, double rt2)
59{
60 switch (i) {
61 case 0:
62 return Ahydro[0] * T + Ahydro[1] * sqrt(T) + Ahydro[2] + (Ahydro[3] + Ahydro[4] * rt) * rt;
63 case 1:
64 return Ahydro[5] * T + Ahydro[6] + rt * (Ahydro[7] + Ahydro[8] * rt);
65 case 2:
66 return Ahydro[9] * T + Ahydro[10] + Ahydro[11] * rt;
67 case 3:
68 return Ahydro[12];
69 case 4:
70 return rt*(Ahydro[13] + Ahydro[14]*rt);
71 case 5:
72 return Ahydro[15]*rt;
73 case 6:
74 return rt*(Ahydro[16] + Ahydro[17]*rt);
75 case 7:
76 return Ahydro[18]*rt2;
77 case 8:
78 return rt2*(Ahydro[19] + Ahydro[20]*rt);
79 case 9:
80 return rt2*(Ahydro[21] + Ahydro[22]*rt2);
81 case 10:
82 return rt2*(Ahydro[23] + Ahydro[24]*rt);
83 case 11:
84 return rt2*(Ahydro[25] + Ahydro[26]*rt2);
85 case 12:
86 return rt2*(Ahydro[27] + Ahydro[28]*rt);
87 case 13:
88 return rt2*(Ahydro[29] + Ahydro[30]*rt + Ahydro[31]*rt2);
89 default:
90 return 0.0;
91 }
92}
93
94double hydrogen::Cprime(int i, double rt, double rt2, double rt3)
95{
96 switch (i) {
97 case 0:
98 return Ahydro[0] + 0.5*Ahydro[1]/sqrt(T) - (Ahydro[3] + 2.0*Ahydro[4]*rt)*rt2;
99 case 1:
100 return Ahydro[5] - rt2*(Ahydro[7] + 2.0*Ahydro[8]*rt);
101 case 2:
102 return Ahydro[9] - Ahydro[11]*rt2;
103 case 3:
104 return 0.0;
105 case 4:
106 return -rt2*(Ahydro[13] + 2.0*Ahydro[14]*rt);
107 case 5:
108 return -Ahydro[15]*rt2;
109 case 6:
110 return -rt2*(Ahydro[16] + 2.0*Ahydro[17]*rt);
111 case 7:
112 return -2.0*Ahydro[18]*rt3;
113 case 8:
114 return -rt3*(2.0*Ahydro[19] + 3.0*Ahydro[20]*rt);
115 case 9:
116 return -rt3*(2.0*Ahydro[21] + 4.0*Ahydro[22]*rt2);
117 case 10:
118 return -rt3*(2.0*Ahydro[23] + 3.0*Ahydro[24]*rt);
119 case 11:
120 return -rt3*(2.0*Ahydro[25] + 4.0*Ahydro[26]*rt2);
121 case 12:
122 return -rt3*(2.0*Ahydro[27] + 3.0*Ahydro[28]*rt);
123 case 13:
124 return -rt3*(2.0*Ahydro[29] + 3.0*Ahydro[30]*rt + 4.0*Ahydro[31]*rt2);
125 default:
126 return 0.0;
127 }
128}
129
130double hydrogen::W(int n, double egrho)
131{
132 return (n == 0 ? (1.0 - egrho)/(2.0*Gamma) :
133 (n*W(n-1, egrho) - 0.5*pow(Rho,2*n)*egrho)/Gamma);
134}
135
136double hydrogen::H(int i, double egrho)
137{
138 return (i < 8 ? pow(Rho,i+2) : pow(Rho,2*i-13)*egrho);
139}
140
141double hydrogen::I(int i, double egrho)
142{
143 return (i < 8 ? pow(Rho,i+1)/double(i+1) : W(i-8, egrho));
144}
145
146double hydrogen::icv(int i, double x, double xlg)
147{
148 return (i == 0 ? x - 1 : x*pow(xlg,i) - i*icv(i-1,x,xlg));
149}
150
152{
153 double rt = 1.0/T;
154 double rt2 = rt*rt;
155 double rt3 = rt*rt2;
156 double egrho = exp(-Gamma*Rho*Rho);
157 double sum = u0;
158 for (int i=0; i<14; i++) {
159 sum += (C(i, rt, rt2) - T*Cprime(i, rt, rt2, rt3))*I(i, egrho);
160 }
161
162 // add \int c_{v,0} term
163 sum += Ghydro[0] * (std::min(T, T1) - To);
164 if (T > T1) {
165 double x = std::min(T, T2) / T1;
166 for (int i = 0; i < 12; i++) {
167 sum += Ghydro[i] * T1 * icv(i, x, log(x));
168 }
169 }
170 if (T > T2) {
171 double x = T/T2;
172 for (int i = 0; i < 5; i++) {
173 sum += Ghydro[i+12] * T2 * icv(i, x, log(x));
174 }
175 }
176 return sum + m_energy_offset;
177}
178
180{
181 double rt = 1.0/T;
182 double rt2 = rt*rt;
183 double rt3 = rt*rt2;
184 double egrho = exp(-Gamma*Rho*Rho);
185 double sum = s0 - R*log(Rho);
186 for (int i=0; i<14; i++) {
187 sum -= Cprime(i, rt, rt2, rt3)*I(i, egrho);
188 }
189
190 // add \int c_{v,0}/T term
191 sum += Ghydro[0] * log(std::min(T, T1)/ To);
192 if (T > T1) {
193 double xlg = log(std::min(T, T2)/T1);
194 for (int i = 0; i < 12; i++) {
195 sum += Ghydro[i] / (i + 1) * pow(xlg, i+1);
196 }
197 }
198 if (T > T2) {
199 double xlg = log(T/T2);
200 for (int i = 0; i < 5; i++) {
201 sum += Ghydro[i+12] / (i + 1) * pow(xlg, i+1);
202 }
203 }
204 return sum + m_entropy_offset;
205}
206
207double hydrogen::Pp()
208{
209 double rt = 1.0/T;
210 double rt2 = rt*rt;
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
221{
222 if ((T < Tmn) || (T > Tc)) {
223 throw CanteraError("hydrogen::ldens",
224 "Temperature out of range. T = {}", T);
225 }
226 double x=1-T/Tc;
227 double sum;
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 return sum+Roc+Dhydro[0]*pow(x,alpha1);
233}
234
236{
237 double x = (1.0 - Tt/T)/(1.0 - Tt/Tc);
238 double result;
239 if ((T < Tmn) || (T > Tc)) {
240 throw CanteraError("hydrogen::Psat",
241 "Temperature out of range. T = {}", T);
242 }
243 result = Fhydro[0]*x + Fhydro[1]*x*x + Fhydro[2]*x*x*x +
244 Fhydro[3]*x*pow(1-x, alpha);
245 return exp(result)*Pt;
246}
247
249{
250 return Tc;
251}
253{
254 return Pc;
255}
257{
258 return 1.0/Roc;
259}
261{
262 return Tmn;
263}
265{
266 return Tmx;
267}
269{
270 return M;
271}
272
273}
Base class for exceptions thrown by Cantera classes.
double x()
Vapor mass fraction.
Definition Sub.cpp:262
double P()
Pressure [Pa].
Definition Sub.cpp:35
double Tmax() override
Maximum temperature for which the equation of state is valid.
Definition Hydrogen.cpp:264
double up() override
Internal energy of a single-phase state.
Definition Hydrogen.cpp:151
double ldens() override
Liquid density. Equation D4 in Reynolds TPSI.
Definition Hydrogen.cpp:220
double Tmin() override
Minimum temperature for which the equation of state is valid.
Definition Hydrogen.cpp:260
double Tcrit() override
Critical temperature [K].
Definition Hydrogen.cpp:248
double sp() override
Entropy of a single-phase state.
Definition Hydrogen.cpp:179
double MolWt() override
Molecular weight [kg/kmol].
Definition Hydrogen.cpp:268
double Vcrit() override
Critical specific volume [m^3/kg].
Definition Hydrogen.cpp:256
double Pcrit() override
Critical pressure [Pa].
Definition Hydrogen.cpp:252
double Psat() override
Saturation pressure. Equation s3 in Reynolds TPSI.
Definition Hydrogen.cpp:235
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
Contains declarations for string manipulation functions within Cantera.