1 //! @file Oxygen.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 "Oxygen.h"
7 #include "cantera/base/stringUtils.h"
8
9 using namespace Cantera;
10
11 namespace tpx
12 {
13 static const double
14 M = 31.9994,
15 Tmn = 54.34,
16 Tmx = 2000.0,
17 Tc = 154.581,
18 Pc = 5.0429e6,
19 Roc = 436.15,
20 R = 2.59820853437877e2,
21 Gamma = 5.46895508389297e-6,
22 alpha = 1.91576,
23 beta = 2239.18105,
24 u0 = 198884.2435,
25 s0 = 668.542976;
26
27 static const double Aoxy[] = {
28 -4.26396872798684e-1, 3.48334938784107e1, -5.77516910418738e2,
29 2.40961751553325e4, -1.23332307855543e6, 3.73585286319658e-4,
30 -1.70178244046465e-1 ,-3.33226903068473e-4, 8.61334799901291e3,
31 -6.80394661057309e-7, 7.09583347162704e-4, -5.73905688255053e-2,
32 -1.92123080811409e-7, 3.11764722329504e-8, -8.09463854745591e-6,
33 -2.22562296356501e-11, 9.18401045361994e-15, 5.75758417511114e-12,
34 -2.10752269644774e-15, 3.62884761272184e3, -1.23317754317110e6,
35 -5.03800414800672e-2, 3.30686173177055e2, -5.26259633964252e-8 ,
36 5.53075442383100e-6, -2.71042853363688e-13, -1.65732450675251e-9 ,
37 -5.82711196409204e-20, 4.42953322148281e-17 ,-2.95529679136244e-25,
38 -1.92361786708846e-23, 9.43758410350413e-23
39 };
40
41 static const double Foxy[] = {
42 -5.581932039e2, -1.0966262185e2, -8.3456211630e-2,
43 2.6603644330e-3, 1.6875023830e-5, -2.1262477120e-7,
44 9.5741096780e-10, -1.6617640450e-12, 2.7545605710e1
45 };
46
47 static const double Doxy[] =
48 { 4.3615175e2, 7.5897189e2, -4.2576866e2, 2.3487106e3, -3.0474660e3, 1.4850169e3 };
49
50 static const double Goxy[] = {
51 -1.29442711174062e6, 5.98231747005341e4, -8.97850772730944e2,
52 6.55236176900400e2, -1.13131252131570e-2,
53 3.4981070244228e-6, 4.21065222886885e-9, 2.67997030050139e2
54 };
55
C(int i,double rt,double rt2)56 double oxygen::C(int i, double rt, double rt2)
57 {
58 switch (i) {
59 case 0:
60 return Aoxy[0] * T + Aoxy[1] * sqrt(T) + Aoxy[2] + (Aoxy[3] + Aoxy[4] * rt) * rt;
61 case 1:
62 return Aoxy[5] * T + Aoxy[6] + rt * (Aoxy[7] + Aoxy[8] * rt);
63 case 2:
64 return Aoxy[9] * T + Aoxy[10] + Aoxy[11] * rt;
65 case 3:
66 return Aoxy[12];
67 case 4:
68 return rt*(Aoxy[13] + Aoxy[14]*rt);
69 case 5:
70 return Aoxy[15]*rt;
71 case 6:
72 return rt*(Aoxy[16] + Aoxy[17]*rt);
73 case 7:
74 return Aoxy[18]*rt2;
75 case 8:
76 return rt2*(Aoxy[19] + Aoxy[20]*rt);
77 case 9:
78 return rt2*(Aoxy[21] + Aoxy[22]*rt2);
79 case 10:
80 return rt2*(Aoxy[23] + Aoxy[24]*rt);
81 case 11:
82 return rt2*(Aoxy[25] + Aoxy[26]*rt2);
83 case 12:
84 return rt2*(Aoxy[27] + Aoxy[28]*rt);
85 case 13:
86 return rt2*(Aoxy[29] + Aoxy[30]*rt + Aoxy[31]*rt2);
87 default:
88 return 0.0;
89 }
90 }
91
Cprime(int i,double rt,double rt2,double rt3)92 double oxygen::Cprime(int i, double rt, double rt2, double rt3)
93 {
94 switch (i) {
95 case 0:
96 return Aoxy[0] + 0.5*Aoxy[1]/sqrt(T) - (Aoxy[3] + 2.0*Aoxy[4]*rt)*rt2;
97 case 1:
98 return Aoxy[5] - rt2*(Aoxy[7] + 2.0*Aoxy[8]*rt);
99 case 2:
100 return Aoxy[9] - Aoxy[11]*rt2;
101 case 3:
102 return 0.0;
103 case 4:
104 return -rt2*(Aoxy[13] + 2.0*Aoxy[14]*rt);
105 case 5:
106 return -Aoxy[15]*rt2;
107 case 6:
108 return -rt2*(Aoxy[16] + 2.0*Aoxy[17]*rt);
109 case 7:
110 return -2.0*Aoxy[18]*rt3;
111 case 8:
112 return -rt3*(2.0*Aoxy[19] + 3.0*Aoxy[20]*rt);
113 case 9:
114 return -rt3*(2.0*Aoxy[21] + 4.0*Aoxy[22]*rt2);
115 case 10:
116 return -rt3*(2.0*Aoxy[23] + 3.0*Aoxy[24]*rt);
117 case 11:
118 return -rt3*(2.0*Aoxy[25] + 4.0*Aoxy[26]*rt2);
119 case 12:
120 return -rt3*(2.0*Aoxy[27] + 3.0*Aoxy[28]*rt);
121 case 13:
122 return -rt3*(2.0*Aoxy[29] + 3.0*Aoxy[30]*rt + 4.0*Aoxy[31]*rt2);
123 default:
124 return 0.0;
125 }
126 }
127
W(int n,double egrho)128 double oxygen::W(int n, double egrho)
129 {
130 return (n == 0 ? (1.0 - egrho)/(2.0*Gamma) :
131 (n*W(n-1, egrho) - 0.5*pow(Rho,2*n)*egrho)/Gamma);
132 }
133
H(int i,double egrho)134 double oxygen::H(int i, double egrho)
135 {
136 return (i < 8 ? pow(Rho,i+2) : pow(Rho,2*i-13)*egrho);
137 }
138
I(int i,double egrho)139 double oxygen::I(int i, double egrho)
140 {
141 return (i < 8 ? pow(Rho,i+1)/double(i+1) : W(i-8, egrho));
142 }
143
up()144 double oxygen::up()
145 {
146 double rt = 1.0/T;
147 double rt2 = rt*rt;
148 double rt3 = rt*rt2;
149 double egrho = exp(-Gamma*Rho*Rho);
150
151 double sum = 0.0;
152 for (int i=0; i<14; i++) {
153 sum += (C(i,rt,rt2) - T*Cprime(i,rt,rt2,rt3))*I(i,egrho);
154 }
155 sum += (((0.25*Goxy[6]*T + Goxy[5]/3.0)*T + 0.5*Goxy[4])*T + Goxy[3])*T + Goxy[2]*log(T)
156 - (Goxy[1] + 0.5*Goxy[0]*rt)*rt + Goxy[7]*beta/(exp(beta*rt) - 1.0) + u0;
157 return sum + m_energy_offset;
158 }
159
sp()160 double oxygen::sp()
161 {
162 double rt = 1.0/T;
163 double rt2 = rt*rt;
164 double rt3 = rt*rt2;
165 double egrho = exp(-Gamma*Rho*Rho);
166
167 double sum = 0.0;
168 sum = s0 - R*log(Rho);
169 for (int i=0; i<14; i++) {
170 sum -= Cprime(i,rt,rt2,rt3)*I(i,egrho);
171 }
172 sum += (((Goxy[6]/3.0)*T + 0.5*Goxy[5])*T + Goxy[4])*T + Goxy[3]*log(T)
173 -((Goxy[0]*rt/3.0 + 0.5*Goxy[1])*rt + Goxy[2])*rt
174 + Goxy[7]*(beta*rt + beta*rt/(exp(beta*rt) - 1.0)
175 - log(exp(beta*rt) - 1.0));
176 return sum + m_entropy_offset;
177 }
178
Pp()179 double oxygen::Pp()
180 {
181 double rt = 1.0/T;
182 double rt2 = rt*rt;
183 double egrho = exp(-Gamma*Rho*Rho);
184
185 double P = Rho*R*T;
186 for (int i=0; i<14; i++) {
187 P += C(i,rt,rt2)*H(i,egrho);
188 }
189 return P;
190 }
191
Psat()192 double oxygen::Psat()
193 {
194 double lnp;
195 int i;
196 if ((T < Tmn) || (T > Tc)) {
197 throw CanteraError("oxygen::Psat",
198 "Temperature out of range. T = {}", T);
199 }
200 for (i=0, lnp=0; i<=7; i++) {
201 if (i==3) {
202 lnp+=Foxy[i]*pow(Tc-T, alpha);
203 } else {
204 lnp+=Foxy[i]*pow(T,i-1);
205 }
206 }
207 lnp+=Foxy[8]*log(T);
208 return exp(lnp);
209 }
210
ldens()211 double oxygen::ldens()
212 {
213 double xx=1-T/Tc, sum=0;
214 if ((T < Tmn) || (T > Tc)) {
215 throw CanteraError("oxygen::ldens",
216 "Temperature out of range. T = {}", T);
217 }
218 for (int i=0; i<=5; i++) {
219 sum+=Doxy[i]*pow(xx,double(i)/3.0);
220 }
221 return sum;
222 }
223
Tcrit()224 double oxygen::Tcrit()
225 {
226 return Tc;
227 }
Pcrit()228 double oxygen::Pcrit()
229 {
230 return Pc;
231 }
Vcrit()232 double oxygen::Vcrit()
233 {
234 return 1.0/Roc;
235 }
Tmin()236 double oxygen::Tmin()
237 {
238 return Tmn;
239 }
Tmax()240 double oxygen::Tmax()
241 {
242 return Tmx;
243 }
MolWt()244 double oxygen::MolWt()
245 {
246 return M;
247 }
248
249 }
250