1 // Reaktoro is a unified framework for modeling chemically reactive systems.
2 //
3 // Copyright (C) 2014-2015 Allan Leal
4 //
5 // This program is free software: you can redistribute it and/or modify
6 // it under the terms of the GNU General Public License as published by
7 // the Free Software Foundation, either version 3 of the License, or
8 // (at your option) any later version.
9 //
10 // This program is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU General Public License for more details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with this program. If not, see <http://www.gnu.org/licenses/>.
17
18 #include "WaterHelmholtzStateHGK.hpp"
19 #include "WaterHelmholtzState.hpp"
20
21 namespace ThermoFun {
22 namespace {
23
24 // Reference temperature of water in units of K
25 const double referenceTemperature = 647.27;
26
27 // Reference density of water in units of kg/m3
28 const double referenceDensity = 317.763;
29
30 // Reference volume of water in units of m3/kg
31 const double referenceVolume = 1.0/referenceDensity;
32
33 // Reference pressure of water in units of Pa
34 const double referencePressure = 22.115e+06;
35
36 // Reference viscosity of water in units of Pa*s
37 const double referenceViscosity = 55.071e-06;
38
39 // Reference thermal conductivity of water in units of W/(K*m)
40 const double referenceConductivity = 0.49450;
41
42 // Reference surface tension of water in units of N/m
43 const double referenceSurfTension = 235.8e-03;
44
45 // Reference constant for Helmholtz function in units of J/kg
46 const double referenceHelmholtz = 69595.89;
47
48 // Reference constant for entropy specific heat in units of J/(kg*K)
49 const double referenceEntropy = 107.5222;
50
51 // Reference constant for sound speed in units of m/s
52 const double referenceSoundSpeed = 263.810;
53
54 const double A0[] =
55 {
56 -0.130840393653E+2,
57 -0.857020420940E+2,
58 0.765192919131E-2,
59 -0.620600116069E+0,
60 -0.106924329402E+2,
61 -0.280671377296E+1,
62 0.119843634845E+3,
63 -0.823907389256E+2,
64 0.555864146443E+2,
65 -0.310698122980E+2,
66 0.136200239305E+2,
67 -0.457116129409E+1,
68 0.115382128188E+1,
69 -0.214242224683E+0,
70 0.282800597384E-1,
71 -0.250384152737E-2,
72 0.132952679669E-3,
73 -0.319277411208E-5
74 };
75
76 const double A1[] =
77 {
78 0.15383053E+1,
79 -0.81048367E+0,
80 -0.68305748E+1,
81 0.00000000E+0,
82 0.86756271E+0
83 };
84
85 const double A20 = 0.42923415E+1;
86
87 const double yc[] =
88 {
89 0.59402227E-1,
90 -0.28128238E-1,
91 0.56826674E-3,
92 -0.27987451E-3
93 };
94
95 const double z0 = 0.317763E+0;
96
97 const int ki[] =
98 {
99 1, 1, 1, 1, 2, 2, 2, 2,
100 3, 3, 3, 3, 4, 4, 4, 4,
101 5, 5, 5, 5, 6, 6, 6, 6,
102 7, 7, 7, 7, 9, 9, 9, 9,
103 3, 3, 1, 5
104 };
105
106 const int li[] =
107 {
108 1, 2, 4, 6, 1, 2, 4, 6,
109 1, 2, 4, 6, 1, 2, 4, 6,
110 1, 2, 4, 6, 1, 2, 4, 6,
111 1, 2, 4, 6, 1, 2, 4, 6,
112 0, 3, 3, 3
113 };
114
115 const double A3[] =
116 {
117 -0.76221190138079E+1,
118 0.32661493707555E+2,
119 0.11305763156821E+2,
120 -0.10015404767712E+1,
121 0.12830064355028E+3,
122 -0.28371416789846E+3,
123 0.24256279839182E+3,
124 -0.99357645626725E+2,
125 -0.12275453013171E+4,
126 0.23077622506234E+4,
127 -0.16352219929859E+4,
128 0.58436648297764E+3,
129 0.42365441415641E+4,
130 -0.78027526961828E+4,
131 0.38855645739589E+4,
132 -0.91225112529381E+3,
133 -0.90143895703666E+4,
134 0.15196214817734E+5,
135 -0.39616651358508E+4,
136 -0.72027511617558E+3,
137 0.11147126705990E+5,
138 -0.17412065252210E+5,
139 0.99918281207782E+3,
140 0.33504807153854E+4,
141 -0.64752644922631E+4,
142 0.98323730907847E+4,
143 0.83877854108422E+3,
144 -0.27919349903103E+4,
145 0.11112410081192E+4,
146 -0.17287587261807E+4,
147 -0.36233262795423E+3,
148 0.61139429010144E+3,
149 0.32968064728562E+2,
150 0.10411239605066E+3,
151 -0.38225874712590E+2,
152 -0.20307478607599E+3
153 };
154
155 const int mi[] = { 2, 2, 2, 4 };
156
157 const int ni[] = { 0, 2, 0, 0 };
158
159 const double alpha[] = { 34, 40, 30, 1050 };
160
161 const double beta[] = { 20000, 20000, 40000, 25 };
162
163 const double ri[] =
164 {
165 0.10038928E+1,
166 0.10038928E+1,
167 0.10038928E+1,
168 0.48778492E+1
169 };
170
171 const double ti[] =
172 {
173 0.98876821E+0,
174 0.98876821E+0,
175 0.99124013E+0,
176 0.41713659E+0
177 };
178
179 const double A4[] =
180 {
181 -0.32329494E-2,
182 -0.24139355E-1,
183 0.79027651E-3,
184 -0.13362857E+1
185 };
186
calculateWaterHelmholtzStateHGK0(Reaktoro_::ThermoScalar t,Reaktoro_::ThermoScalar)187 auto calculateWaterHelmholtzStateHGK0(Reaktoro_::ThermoScalar t, Reaktoro_::ThermoScalar /*d*/) -> WaterHelmholtzState
188 {
189 WaterHelmholtzState s;
190
191 const auto ln_t = log(t);
192
193 s.helmholtz = (A0[0] + A0[1] * t) * ln_t;
194 s.helmholtzT = A0[0]/t + A0[1]*(ln_t + 1);
195 s.helmholtzTT = -A0[0]/(t*t) + A0[1]/t;
196 s.helmholtzTTT = 2*A0[0]/(t*t*t) - A0[1]/(t*t);
197
198 for(int i = 2; i <= 17; ++i)
199 {
200 const auto aux = A0[i] * pow(t, i - 4);
201
202 s.helmholtz += aux;
203 s.helmholtzT += aux * (i - 4)/t;
204 s.helmholtzTT += aux * (i - 4)*(i - 5)/(t*t);
205 s.helmholtzTTT += aux * (i - 4)*(i - 5)*(i - 6)/(t*t*t);
206 }
207
208 return s;
209 }
210
calculateWaterHelmholtzStateHGK1(Reaktoro_::ThermoScalar t,Reaktoro_::ThermoScalar d)211 auto calculateWaterHelmholtzStateHGK1(Reaktoro_::ThermoScalar t, Reaktoro_::ThermoScalar d) -> WaterHelmholtzState
212 {
213 WaterHelmholtzState s;
214
215 for(int i = 0; i <= 4; ++i)
216 {
217 const auto aux = d * A1[i] * pow(t, 1 - i);
218
219 s.helmholtz += aux;
220 s.helmholtzT -= aux * (i - 1)/t;
221 s.helmholtzTT += aux * (i - 1)*i/(t*t);
222 s.helmholtzTTT -= aux * (i - 1)*i*(i + 1)/(t*t*t);
223 }
224
225 s.helmholtzD = s.helmholtz/d;
226 s.helmholtzTD = s.helmholtzT/d;
227 s.helmholtzTTD = s.helmholtzTT/d;
228
229 return s;
230 }
231
calculateWaterHelmholtzStateHGK2(Reaktoro_::ThermoScalar t,Reaktoro_::ThermoScalar d)232 auto calculateWaterHelmholtzStateHGK2(Reaktoro_::ThermoScalar t, Reaktoro_::ThermoScalar d) -> WaterHelmholtzState
233 {
234 WaterHelmholtzState s;
235
236 const auto t3 = pow(t, -3);
237 const auto t5 = pow(t, -5);
238 const auto ln_t = log(t);
239
240 const auto y = d * (yc[0] + yc[1]*ln_t + yc[2]*t3 + yc[3]*t5);
241 const auto y_r = y/d;
242 const auto y_t = d * (yc[1] - 3.0*yc[2]*t3 - 5.0*yc[3]*t5)/t;
243 const auto y_rr = 0.0;
244 const auto y_tt = d * (-yc[1] + 12.0*yc[2]*t3 + 30.0*yc[3]*t5)/(t*t);
245 const auto y_rt = y_t/d;
246 const auto y_rrr = 0.0;
247 const auto y_rrt = y_rt/d - y_t/(d*d);
248 const auto y_rtt = y_tt/d;
249 const auto y_ttt = d * (2*yc[1] - 60*yc[2]*t3 - 210*yc[3]*t5)/(t*t*t);
250
251 const auto x = 1.0/(1.0 - y);
252 const auto x2 = x * x;
253 const auto x_r = y_r * x2;
254 const auto x_t = y_t * x2;
255 const auto x_rr = y_rr * x2 + 2.0 * y_r * x_r * x;
256 const auto x_tt = y_tt * x2 + 2.0 * y_t * x_t * x;
257 const auto x_rt = y_rt * x2 + 2.0 * y_r * x_t * x;
258 const auto x_rrr = y_rrr*x2 + 4*y_rr*x_r*x + 2*y_r*(x_rr*x + x_r*x_r);
259 const auto x_rrt = y_rrt*x2 + 2*(y_rt*x_r + y_rr*x_t) + 2*y_r*(x_rt*x + x_r*x_t);
260 const auto x_rtt = y_rtt*x2 + 4*y_rt*x_t*x + 2*y_r*(x_tt*x + x_t*x_t);
261 const auto x_ttt = y_ttt*x2 + 4*y_tt*x_t*x + 2*y_t*(x_tt*x + x_t*x_t);
262
263 const auto u = log(d * x);
264 const auto u_r = x_r/x + 1.0/d;
265 const auto u_t = x_t/x;
266 const auto u_rr = x_rr/x - x_r*x_r/(x*x) - 1.0/(d*d);
267 const auto u_rt = x_rt/x - x_r*x_t/(x*x);
268 const auto u_tt = x_tt/x - x_t*x_t/(x*x);
269 const auto u_rrr = x_rrr/x - 3*x_rr*x_r/(x*x) + 2*x_r*x_r*x_r/(x*x*x) + 2/(d*d*d);
270 const auto u_rrt = x_rrt/x - (2*x_rt*x_r + x_rr*x_t)/(x*x) + 2*x_r*x_r*x_t/(x*x*x);
271 const auto u_rtt = x_rtt/x - (2*x_rt*x_t + x_tt*x_r)/(x*x) + 2*x_t*x_t*x_r/(x*x*x);
272 const auto u_ttt = x_ttt/x - 3*x_tt*x_t/(x*x) + 2*x_t*x_t*x_t/(x*x*x);
273
274 const auto c1 = -130.0/3.0;
275 const auto c2 = 169.0/6.0;
276 const auto c3 = -14.0;
277
278 s.helmholtz = A20 * t * (u + c1*x + c2*x*x + c3*y);
279 s.helmholtzD = A20 * t * (u_r + c1*x_r + 2*c2*x*x_r + c3*y_r);
280 s.helmholtzT = A20 * t * (u_t + c1*x_t + 2*c2*x*x_t + c3*y_t) + s.helmholtz/t;
281 s.helmholtzDD = A20 * t * (u_rr + c1*x_rr + 2*c2*(x*x_rr + x_r*x_r) + c3*y_rr);
282 s.helmholtzTD = A20 * t * (u_rt + c1*x_rt + 2*c2*(x*x_rt + x_r*x_t) + c3*y_rt) + s.helmholtzD/t;
283 s.helmholtzTT = A20 * t * (u_tt + c1*x_tt + 2*c2*(x*x_tt + x_t*x_t) + c3*y_tt) + 2*(s.helmholtzT/t - s.helmholtz/(t*t));
284 s.helmholtzDDD = A20 * t * (u_rrr + c1*x_rrr + 2*c2*(3*x_r*x_rr + x*x_rrr) + c3*y_rrr);
285 s.helmholtzTDD = A20 * t * (u_rrt + c1*x_rrt + 2*c2*(x_t*x_rr + 2*x_r*x_rt + x*x_rrt) + c3*y_rrt) + s.helmholtzDD/t;
286 s.helmholtzTTD = A20 * t * (u_rtt + c1*x_rtt + 2*c2*(x_r*x_tt + 2*x_t*x_rt + x*x_rtt) + c3*y_rtt) + 2*(s.helmholtzTD - s.helmholtzD/t)/t;
287 s.helmholtzTTT = A20 * t * (u_ttt + c1*x_ttt + 2*c2*(3*x_t*x_tt + x*x_ttt) + c3*y_ttt) + 3*(s.helmholtzTT - 2*s.helmholtzT/t + s.helmholtz/(t*t))/t;
288
289 return s;
290 }
291
calculateWaterHelmholtzStateHGK3(Reaktoro_::ThermoScalar t,Reaktoro_::ThermoScalar d)292 auto calculateWaterHelmholtzStateHGK3(Reaktoro_::ThermoScalar t, Reaktoro_::ThermoScalar d) -> WaterHelmholtzState
293 {
294 WaterHelmholtzState s;
295
296 const auto z = 1 - exp(-z0 * d);
297 const auto z_r = z0 * (1 - z);
298 const auto z_rr = -z0 * z_r;
299 const auto z_rrr = -z0 * z_rr;
300
301 for(int i = 0; i <= 35; ++i)
302 {
303 const auto lambda = A3[i] * pow(t, -li[i]) * pow(z, ki[i]);
304 const auto lambda_r = ki[i]*z_r*lambda/z;
305 const auto lambda_t = -li[i]*lambda/t;
306 const auto lambda_rr = lambda_r*(z_rr/z_r + lambda_r/lambda - z_r/z);
307 const auto lambda_rt = lambda_r*lambda_t/lambda;
308 const auto lambda_tt = lambda_t*(lambda_t/lambda - 1.0/t);
309 const auto lambda_rrr = lambda_rr*(z_rr/z_r + lambda_r/lambda - z_r/z) + lambda_r*(z_rrr/z_r - pow(z_rr/z_r, 2) + lambda_rr/lambda - pow(lambda_r/lambda, 2) - z_rr/z + pow(z_r/z, 2));
310 auto l_pow = pow(lambda_t/lambda, 2); if (l_pow == 0) {l_pow.ddt = 0; l_pow.ddp = 0;}
311 const auto lambda_rrt = -pow(lambda_r/lambda, 2)*lambda_t + (lambda_rr*lambda_t + lambda_rt*lambda_r)/lambda;
312 const auto lambda_rtt = -l_pow*lambda_r + (lambda_tt*lambda_r + lambda_rt*lambda_t)/lambda;
313 const auto lambda_ttt = lambda_tt * (lambda_t/lambda - 1.0/t) + lambda_t*(lambda_tt/lambda - l_pow + 1.0/(t*t));
314
315 s.helmholtz += lambda;
316 s.helmholtzD += lambda_r;
317 s.helmholtzT += lambda_t;
318 s.helmholtzDD += lambda_rr;
319 s.helmholtzTD += lambda_rt;
320 s.helmholtzTT += lambda_tt;
321 s.helmholtzDDD += lambda_rrr;
322 s.helmholtzTDD += lambda_rrt;
323 s.helmholtzTTD += lambda_rtt;
324 s.helmholtzTTT += lambda_ttt;
325 }
326
327 return s;
328 }
329
calculateWaterHelmholtzStateHGK4(Reaktoro_::ThermoScalar t,Reaktoro_::ThermoScalar d)330 auto calculateWaterHelmholtzStateHGK4(Reaktoro_::ThermoScalar t, Reaktoro_::ThermoScalar d) -> WaterHelmholtzState
331 {
332 WaterHelmholtzState s;
333
334 for(int i = 0; i <= 3; ++i)
335 {
336 const auto delta = (d - ri[i])/ri[i];
337 const auto tau = (t - ti[i])/ti[i];
338 const auto delta_r = 1.0/ri[i];
339 const auto tau_t = 1.0/ti[i];
340
341 const auto delta_m = pow(delta, mi[i]);
342 const auto delta_n = pow(delta, ni[i]);
343
344 const auto psi = (ni[i] - alpha[i]*mi[i]*delta_m)*delta_r/delta;
345 const auto psi_r = -(ni[i] + alpha[i]*mi[i]*(mi[i] - 1)*delta_m)*pow(delta_r/delta, 2);
346 const auto psi_rr = (2*ni[i] - alpha[i]*mi[i]*(mi[i] - 1)*(mi[i] - 2)*delta_m)*pow(delta_r/delta, 3);
347
348 const auto theta = A4[i]*delta_n*exp(-alpha[i]*delta_m - beta[i]*tau*tau);
349 const auto theta_r = psi*theta;
350 const auto theta_t = -2*beta[i]*tau*tau_t*theta;
351 const auto theta_rr = psi_r*theta + psi*theta_r;
352 const auto theta_tt = 2*beta[i]*(2*beta[i]*tau*tau - 1)*tau_t*tau_t*theta;
353 const auto theta_rt = -2*beta[i]*tau*tau_t*theta_r;
354 const auto theta_rrr = psi_rr*theta + 2*psi_r*theta_r + psi*theta_rr;
355 const auto theta_rrt = psi_r*theta_r + psi*theta_rt;
356 const auto theta_rtt = psi*theta_tt;
357 const auto theta_ttt = -2*beta[i]*(2*tau_t*tau_t*theta_t + tau*tau_t*theta_tt);
358
359 s.helmholtz += theta;
360 s.helmholtzD += theta_r;
361 s.helmholtzT += theta_t;
362 s.helmholtzDD += theta_rr;
363 s.helmholtzTD += theta_rt;
364 s.helmholtzTT += theta_tt;
365 s.helmholtzDDD += theta_rrr;
366 s.helmholtzTDD += theta_rrt;
367 s.helmholtzTTD += theta_rtt;
368 s.helmholtzTTT += theta_ttt;
369 }
370
371 return s;
372 }
373
374 } // namespace
375
waterHelmholtzStateHGK(Reaktoro_::Temperature T,Reaktoro_::ThermoScalar D)376 auto waterHelmholtzStateHGK(Reaktoro_::Temperature T, Reaktoro_::ThermoScalar D) -> WaterHelmholtzState
377 {
378 // The dimensionless temperature and density
379 const auto t = T/referenceTemperature;
380 const auto r = D/referenceDensity;
381
382 // Compute the contributions from each auxiliary Helmholtz state
383 WaterHelmholtzState aux0, aux1, aux2, aux3, aux4, res;
384 aux0 = calculateWaterHelmholtzStateHGK0(t, r);
385 aux1 = calculateWaterHelmholtzStateHGK1(t, r);
386 aux2 = calculateWaterHelmholtzStateHGK2(t, r);
387 aux3 = calculateWaterHelmholtzStateHGK3(t, r);
388 aux4 = calculateWaterHelmholtzStateHGK4(t, r);
389
390 // Assemble the contributions from each auxiliary Helmholtz state
391 res.helmholtz = aux0.helmholtz + aux1.helmholtz + aux2.helmholtz + aux3.helmholtz + aux4.helmholtz;
392 res.helmholtzD = aux0.helmholtzD + aux1.helmholtzD + aux2.helmholtzD + aux3.helmholtzD + aux4.helmholtzD;
393 res.helmholtzT = aux0.helmholtzT + aux1.helmholtzT + aux2.helmholtzT + aux3.helmholtzT + aux4.helmholtzT;
394 res.helmholtzDD = aux0.helmholtzDD + aux1.helmholtzDD + aux2.helmholtzDD + aux3.helmholtzDD + aux4.helmholtzDD;
395 res.helmholtzTD = aux0.helmholtzTD + aux1.helmholtzTD + aux2.helmholtzTD + aux3.helmholtzTD + aux4.helmholtzTD;
396 res.helmholtzTT = aux0.helmholtzTT + aux1.helmholtzTT + aux2.helmholtzTT + aux3.helmholtzTT + aux4.helmholtzTT;
397 res.helmholtzDDD = aux0.helmholtzDDD + aux1.helmholtzDDD + aux2.helmholtzDDD + aux3.helmholtzDDD + aux4.helmholtzDDD;
398 res.helmholtzTDD = aux0.helmholtzTDD + aux1.helmholtzTDD + aux2.helmholtzTDD + aux3.helmholtzTDD + aux4.helmholtzTDD;
399 res.helmholtzTTD = aux0.helmholtzTTD + aux1.helmholtzTTD + aux2.helmholtzTTD + aux3.helmholtzTTD + aux4.helmholtzTTD;
400 res.helmholtzTTT = aux0.helmholtzTTT + aux1.helmholtzTTT + aux2.helmholtzTTT + aux3.helmholtzTTT + aux4.helmholtzTTT;
401
402 // Convert the Helmholtz free energy of water and its derivatives to dimensioned form
403 res.helmholtz *= referenceHelmholtz;
404 res.helmholtzD *= referenceHelmholtz/referenceDensity;
405 res.helmholtzT *= referenceHelmholtz/referenceTemperature;
406 res.helmholtzDD *= referenceHelmholtz/referenceDensity/referenceDensity;
407 res.helmholtzTD *= referenceHelmholtz/referenceDensity/referenceTemperature;
408 res.helmholtzTT *= referenceHelmholtz/referenceTemperature/referenceTemperature;
409 res.helmholtzDDD *= referenceHelmholtz/referenceDensity/referenceDensity/referenceDensity;
410 res.helmholtzTDD *= referenceHelmholtz/referenceDensity/referenceDensity/referenceTemperature;
411 res.helmholtzTTD *= referenceHelmholtz/referenceDensity/referenceTemperature/referenceTemperature;
412 res.helmholtzTTT *= referenceHelmholtz/referenceTemperature/referenceTemperature/referenceTemperature;
413
414 return res;
415 }
416
417 } // namespace Reaktoro
418