Cantera  2.0
lk.cpp
1 // Lee-Kesler equation of state
2 
3 #include "lk.h"
4 #include <math.h>
5 
6 namespace tpx
7 {
8 
9 static double b[2][4] = {{0.1181193, 0.265728, 0.154790, 0.030323},
10  {0.2026579, 0.331511, 0.027655, 0.203488}
11 };
12 static double c[2][4] = {{0.0236744, 0.0186984, 0.0, 0.042724},
13  {0.0313385, 0.0503618, 0.016901, 0.041577}
14 };
15 static double d[2][2] = {{1.55488e-5, 6.23689e-5},{4.8736e-5, 0.740336e-5}};
16 static double beta[2] = {0.65392, 1.226};
17 static double gamma[2] = {0.060167, 0.03754};
18 
19 //--------------------------- member functions ------------------
20 
21 double leekesler::W(int n, double egrho, double Gamma)
22 {
23  return (n == 0 ? (1.0 - egrho)/(2.0*Gamma) :
24  (n*W(n-1, egrho, Gamma) - 0.5*pow(Rho,2*n)*egrho)/Gamma);
25 }
26 
27 double leekesler::up()
28 {
29  return -(8314.3/Mw)*T*(1.0 + T*I()/Tcr); // + h_0(T)
30 }
31 
32 double leekesler::hdep()
33 {
34  double tr = T/Tcr;
35  return tr*tr*I() + (1.0 - z())*tr;
36 }
37 
38 double leekesler::sdep()
39 {
40  double tr = T/Tcr;
41  return tr*I() + J() - log(z());
42 }
43 
44 double leekesler::sp()
45 {
46  const double Pref = 101325.0;
47  double rgas = 8314.3/Mw;
48  return rgas*(log(Pref/(Rho*rgas*T)) - (T/Tcr)*I() - J());
49 }
50 
51 double leekesler::I() // \int_0^\rho_r (1/\rho_r)(dZ/dT_r) d\rho_r
52 {
53  double Bp, Cp, Dp;
54  double rtr = Tcr/T;
55  double rtr2 = rtr*rtr;
56  double rvr = 8314.3*Tcr*Rho/(Pcr*Mw); // 1/v_r^\prime
57  double rvr2 = rvr*rvr;
58  double egrho;
59 
60  egrho = exp(-gamma[Isr]*rvr2);
61  Bp = rtr2*b[Isr][1] + 2.0*rtr*rtr2*b[Isr][2] + 3.0*rtr2*rtr2*b[Isr][3];
62  Cp = rtr2*c[Isr][1] - 3.0*c[Isr][2]*rtr2*rtr2;
63  Dp = -d[Isr][1]*rtr2;
64  double r = Bp*rvr + 0.5*rvr2*Cp + 0.2*pow(rvr,5)*Dp
65  - 3.0*c[Isr][3]*rtr2*rtr2*(beta[Isr]*W(0,egrho,gamma[Isr])
66  + gamma[Isr]*W(1,egrho,gamma[Isr]));
67  return r;
68 }
69 
70 double leekesler::J() // \int_0^\rho_r (1/\rho_r)(Z - 1) d\rho_r
71 {
72  double BB, CC, DD;
73  double rtr = Tcr/T;
74  double rtr2 = rtr*rtr;
75  double rvr = 8314.3*Tcr*Rho/(Pcr*Mw); // 1/v_r^\prime
76  double rvr2 = rvr*rvr;
77  double egrho;
78 
79  egrho = exp(-gamma[Isr]*rvr2);
80  BB = b[Isr][0] - rtr*(b[Isr][1]
81  + rtr*(b[Isr][2] + rtr*b[Isr][3]));
82  CC = c[Isr][0] - rtr*(c[Isr][1] - c[Isr][2]*rtr*rtr);
83  DD = d[Isr][0] + d[Isr][1]*rtr;
84  double r = BB*rvr + 0.5*rvr2*CC + 0.2*pow(rvr,5)*DD
85  + c[Isr][3]*rtr2*rtr*(beta[Isr]*W(0,egrho,gamma[Isr])
86  + gamma[Isr]*W(1,egrho,gamma[Isr]));
87  return r;
88 }
89 
90 double leekesler::z()
91 {
92  double zz, rvr2, BB, CC, DD, EE;
93  double rtr = Tcr/T; // 1/T_r
94  double rvr = Rho*8314.3*Tcr/(Pcr*Mw);
95  rvr2 = rvr*rvr;
96  BB = b[Isr][0] - rtr*(b[Isr][1]
97  + rtr*(b[Isr][2] + rtr*b[Isr][3]));
98  CC = c[Isr][0] - rtr*(c[Isr][1] - c[Isr][2]*rtr*rtr);
99  DD = d[Isr][0] + d[Isr][1]*rtr;
100  EE = exp(-gamma[Isr]*rvr2);
101 
102  zz = 1.0 + BB*rvr + CC*rvr2 + DD*pow(rvr,5)
103  + c[Isr][3]*pow(rtr,3)*rvr2*
104  (beta[Isr] + gamma[Isr]*rvr2)*EE;
105  return zz;
106 }
107 
108 
109 double leekesler::Pp()
110 {
111  return 8314.3*z()*Rho*T/Mw;
112 }
113 
114 double leekesler::Psat()
115 {
116  double tr = 1.0 - Tcr/T;
117  double lpr;
118 
119  if (Isr == 0) {
120  lpr = 5.395743797*tr + 0.05524287*tr*tr + 0.06853005*tr*tr*tr;
121  } else {
122  lpr = 7.259961465*tr - 0.549206092*tr*tr + 0.177581752*tr*tr*tr;
123  }
124  return Pcr*exp(lpr);
125 }
126 
127 double leekesler::ldens()
128 {
129  double x = 1.0 - T/Tcr;
130 
131  // for simple fluid
132  double rho_r;
133  if (Isr == 0) {
134  rho_r = 5.2307 + 15.16*x - 21.9778*x*x + 18.767*x*x*x;
135  } else {
136  rho_r = 6.166930606 + 17.42866964*x - 18.62589833*x*x
137  + 11.73957224*x*x*x;
138  rho_r *= 1.0;
139  }
140  return Pcr*rho_r*Mw/(8314.3*Tcr);
141 }
142 
143 double leekesler::Tcrit()
144 {
145  return Tcr;
146 }
147 double leekesler::Pcrit()
148 {
149  return Pcr;
150 }
151 double leekesler::Vcrit()
152 {
153  return 0.2901*8314.3*Tcr/(Pcr*Mw);
154 }
155 double leekesler::Tmin()
156 {
157  return -100.0;
158 }
159 double leekesler::Tmax()
160 {
161  return 10000.0;
162 }
163 char* leekesler::name()
164 {
165  return (char*) m_name.c_str();
166 }
167 char* leekesler::formula()
168 {
169  return (char*) m_formula.c_str();
170 }
171 double leekesler::MolWt()
172 {
173  return Mw;
174 }
175 }