Cantera  2.0
Nitrogen.cpp
1 // Nitrogen
2 
3 #include "Nitrogen.h"
4 #include <math.h>
5 #include <iostream>
6 using namespace std;
7 
8 namespace tpx
9 {
10 
11 
12 static const double M = 28.01348,
13  Tmn = 63.15,
14  Tmx = 2000.0,
15  Tc = 126.200,
16  Pc = 3.4e6,
17  Roc = 314.03,
18  To = 63.15,
19  R = 2.96790515164171e2,
20  Gamma = 7.13602531283233e-6,
21  alpha = 1.95,
22  beta = 3353.40610,
23  u0 = 150877.551,
24  s0 = 214.9352518;
25 
26 
27 static const double Ann[] = {
28  1.75889959256970e-1, 1.38197604384933e1, -3.14918412133921e2,
29  4.40300150239380e3, -5.45358971644916e5, 4.84413320182919e-4,
30  -5.18964416491365e-2, 6.57265859197103e-4 ,8.51299771713314e4 ,
31  1.33459405162578e-8, 3.83381319826746e-4, -8.35421151028455e-2,
32  2.84874912286101e-7,
33  -2.38296116270360e-7, -1.48321912935764e-4, 5.62605853190540e-10,
34  -2.98201050924595e-13, 9.85319087685241e-11, -1.92002176056468e-14,
35  -7.82250103373122e4, -5.51801778744598e5, -5.72781957607352e-1,
36  3.25760529488327e2, -1.34659309828737e-6, -1.92036423064911e-5,
37  -3.94564337674524e-12,-2.44388245328965e-9, -1.50970602460077e-18,
38  1.25854885346038e-16,-8.34271144923969e-24, -1.17299202018417e-22,
39  9.06544823455730e-22
40 };
41 
42 static const double Fnn[]= {
43  8.3944094440e3, -1.8785191705e3, -7.2822291650,
44  1.0228509660e-2, 5.5560638250e-4,
45  -5.9445446620e-6, 2.7154339320e-8,
46  -4.8795359040e-11, 5.0953608240e2
47 };
48 
49 static const double Dnn[] = {
50  3.1402991e2, 4.4111015e2, 9.4622994e2 ,
51  -2.9067111e3, 4.4785979e3, -2.2746914e3
52 };
53 
54 static const double Gnn[] = {
55  -2.18203473713518e5, 1.01573580096247e4, -1.65504721657240e2,
56  7.43175999190430e2, -5.14605623546025e-3,
57  5.18347156760489e-6, -1.05922170493616e-9, 2.98389393363817e2
58 };
59 
60 
61 //equation P4
62 double nitrogen::C(int i, double rt, double rt2)
63 {
64  switch (i) {
65  case 0 :
66  return Ann[0] * T + Ann[1] * sqrt(T)
67  + Ann[2] + (Ann[3] + Ann[4] * rt) * rt;
68  case 1 :
69  return Ann[5] * T + Ann[6] + rt * (Ann[7] + Ann[8] * rt);
70  case 2 :
71  return Ann[9] * T + Ann[10] + Ann[11] * rt;
72  case 3 :
73  return Ann[12];
74  case 4 :
75  return rt*(Ann[13] + Ann[14]*rt);
76  case 5 :
77  return Ann[15]*rt;
78  case 6 :
79  return rt*(Ann[16] + Ann[17]*rt);
80  case 7 :
81  return Ann[18]*rt2;
82  case 8 :
83  return rt2*(Ann[19] + Ann[20]*rt);
84  case 9 :
85  return rt2*(Ann[21] + Ann[22]*rt2);
86  case 10 :
87  return rt2*(Ann[23] + Ann[24]*rt);
88  case 11 :
89  return rt2*(Ann[25] + Ann[26]*rt2);
90  case 12 :
91  return rt2*(Ann[27] + Ann[28]*rt);
92  case 13 :
93  return rt2*(Ann[29] + Ann[30]*rt + Ann[31]*rt2);
94  default :
95  return 0.0;
96  }
97 }
98 
99 double nitrogen::Cprime(int i, double rt, double rt2, double rt3)
100 {
101  switch (i) {
102  case 0 :
103  return Ann[0] + 0.5*Ann[1]/sqrt(T) - (Ann[3] + 2.0*Ann[4]*rt)*rt2;
104  case 1 :
105  return Ann[5] - rt2*(Ann[7] + 2.0*Ann[8]*rt);
106  case 2 :
107  return Ann[9] - Ann[11]*rt2;
108  case 3 :
109  return 0.0;
110  case 4 :
111  return -rt2*(Ann[13] + 2.0*Ann[14]*rt);
112  case 5 :
113  return -Ann[15]*rt2;
114  case 6 :
115  return -rt2*(Ann[16] + 2.0*Ann[17]*rt);
116  case 7 :
117  return -2.0*Ann[18]*rt3;
118  case 8 :
119  return -rt3*(2.0*Ann[19] + 3.0*Ann[20]*rt);
120  case 9 :
121  return -rt3*(2.0*Ann[21] + 4.0*Ann[22]*rt2);
122  case 10 :
123  return -rt3*(2.0*Ann[23] + 3.0*Ann[24]*rt);
124  case 11 :
125  return -rt3*(2.0*Ann[25] + 4.0*Ann[26]*rt2);
126  case 12 :
127  return -rt3*(2.0*Ann[27] + 3.0*Ann[28]*rt);
128  case 13 :
129  return -rt3*(2.0*Ann[29] + 3.0*Ann[30]*rt + 4.0*Ann[31]*rt2);
130  default :
131  return 0.0;
132  }
133 }
134 
135 double nitrogen::W(int n, double egrho)
136 {
137  return (n == 0 ? (1.0 - egrho)/(2.0*Gamma) :
138  (n*W(n-1, egrho) - 0.5*pow(Rho,2*n)*egrho)/Gamma);
139 }
140 
141 double nitrogen::H(int i, double egrho)
142 {
143  return (i < 8 ? pow(Rho,i+2) : pow(Rho,2*i-13)*egrho);
144 }
145 
146 double nitrogen::I(int i, double egrho)
147 {
148  return (i < 8 ? pow(Rho,i+1)/double(i+1) : W(i-8, egrho));
149 }
150 
151 double nitrogen::up()
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 
158  double sum = 0.0;
159  for (int i=0; i<14; i++) {
160  sum += (C(i,rt,rt2) - T*Cprime(i,rt,rt2,rt3))*I(i,egrho);
161  }
162 
163  sum += (((0.25*Gnn[6]*T + Gnn[5]/3.0)*T
164  + 0.5*Gnn[4])*T + Gnn[3])*T + Gnn[2]*log(T)
165  - (Gnn[1] + 0.5*Gnn[0]*rt)*rt
166  + Gnn[7]*beta/(exp(beta*rt) - 1.0) + u0
167  + m_energy_offset;
168 
169  return sum;
170 }
171 
172 
173 double nitrogen::sp()
174 {
175  double rt = 1.0/T;
176  double rt2 = rt*rt;
177  double rt3 = rt*rt2;
178  double egrho = exp(-Gamma*Rho*Rho);
179 
180  double sum = 0.0;
181  sum = s0 + m_entropy_offset - R*log(Rho);
182  for (int i=0; i<14; i++) {
183  sum -= Cprime(i,rt,rt2,rt3)*I(i,egrho);
184  }
185 
186  sum += (((Gnn[6]/3.0)*T + 0.5*Gnn[5])*T + Gnn[4])*T + Gnn[3]*log(T)
187  -((Gnn[0]*rt/3.0 + 0.5*Gnn[1])*rt + Gnn[2])*rt
188  + Gnn[7]*(beta*rt + beta*rt/(exp(beta*rt) - 1.0)
189  - log(exp(beta*rt) - 1.0));
190  return sum;
191 }
192 
193 double nitrogen::Pp()
194 {
195  double rt = 1.0/T;
196  double rt2 = rt*rt;
197  double egrho = exp(-Gamma*Rho*Rho);
198 
199  double P = Rho*R*T;
200  for (int i=0; i<14; i++) {
201  P += C(i,rt,rt2)*H(i,egrho);
202  }
203  return P;
204 }
205 
206 //equation s4
207 double nitrogen::Psat()
208 {
209  double lnp;
210  int i;
211  if ((T < Tmn) || (T > Tc)) {
212  set_Err(TempError);
213  }
214  for (i=0, lnp=0; i<=7; i++) {
215  if (i==3) {
216  lnp+=Fnn[i]*pow(Tc-T, alpha);
217  } else {
218  lnp+=Fnn[i]*pow(T,i-1);
219  }
220  }
221  lnp+=Fnn[8]*log(T);
222  return exp(lnp);
223 }
224 
225 //equation D2
226 double nitrogen::ldens()
227 {
228  double xx=1-T/Tc, sum=0;
229  if ((T < Tmn) || (T > Tc)) {
230  set_Err(TempError);
231  }
232  for (int i=0; i<=5; i++) {
233  sum+=Dnn[i]*pow(xx,double(i)/3.0);
234  }
235  return sum;
236 }
237 
238 double nitrogen::Tcrit()
239 {
240  return Tc;
241 }
242 double nitrogen::Pcrit()
243 {
244  return Pc;
245 }
246 double nitrogen::Vcrit()
247 {
248  return 1.0/Roc;
249 }
250 double nitrogen::Tmin()
251 {
252  return Tmn;
253 }
254 double nitrogen::Tmax()
255 {
256  return Tmx;
257 }
258 char* nitrogen::name()
259 {
260  return (char*) m_name.c_str();
261 }
262 char* nitrogen::formula()
263 {
264  return (char*) m_formula.c_str();
265 }
266 double nitrogen::MolWt()
267 {
268  return M;
269 }
270 }