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