1 /**
2 * @file WaterPropsIAPWSphi.cpp
3 * Definitions for Lowest level of the classes which support a real water
4 * model (see class \link Cantera::WaterPropsIAPWS WaterPropsIAPWS\endlink and
5 * class \link Cantera::WaterPropsIAPWSphi WaterPropsIAPWSphi \endlink).
6 */
7
8 // This file is part of Cantera. See License.txt in the top-level directory or
9 // at https://cantera.org/license.txt for license and copyright information.
10
11 #include "cantera/thermo/WaterPropsIAPWSphi.h"
12 #include "cantera/base/global.h"
13
14 #include <cmath>
15 #include <algorithm>
16
17 namespace Cantera
18 {
19
20 using std::sqrt;
21 using std::log;
22 using std::exp;
23 using std::pow;
24 using std::fabs;
25
26 /*
27 * The added constants were calculated so that u = s = 0
28 * for liquid at the triple point. These where determined
29 * by the program testPress. I'm not quite satisfied with
30 * the result, but will let it stand for the moment.
31 * H didn't turn out to be .611872 J/kg, but .611782 J/kg.
32 * There may be a slight error here somehow.
33 */
34 // \cond
35 static const doublereal ni0[9] = {
36 0.0,
37 -8.32044648201 - 0.000000001739715,
38 6.6832105268 + 0.000000000793232,
39 3.00632,
40 0.012436,
41 0.97315,
42 1.27950,
43 0.96956,
44 0.24873
45 };
46
47 static const doublereal gammi0[9] = {
48 0.0,
49 0.0,
50 0.0,
51 0.0,
52 1.28728967,
53 3.53734222,
54 7.74073708,
55 9.24437796,
56 27.5075105
57 };
58
59 static const int ciR[56] = {
60 0, // 0
61 0, // 1
62 0,
63 0,
64 0,
65 0, // 5
66 0,
67 0,
68 1,
69 1,
70 1, // 10
71 1,
72 1,
73 1,
74 1,
75 1, // 15
76 1,
77 1,
78 1,
79 1,
80 1, // 20
81 1,
82 1,
83 2,
84 2,
85 2, // 25
86 2,
87 2,
88 2,
89 2,
90 2, // 30
91 2,
92 2,
93 2,
94 2,
95 2, // 35
96 2,
97 2,
98 2,
99 2,
100 2, // 40
101 2,
102 2,
103 3,
104 3,
105 3, // 45
106 3,
107 4,
108 6,
109 6,
110 6, // 50
111 6,
112 0,
113 0,
114 0,
115 0 // 55
116 };
117
118 static const int diR[55] = {
119 0, // 0
120 1, // 1
121 1,
122 1,
123 2,
124 2, // 5
125 3,
126 4,
127 1,
128 1,
129 1, // 10
130 2,
131 2,
132 3,
133 4,
134 4, // 15
135 5,
136 7,
137 9,
138 10,
139 11, // 20
140 13,
141 15,
142 1,
143 2,
144 2, // 25
145 2,
146 3,
147 4,
148 4,
149 4, // 30
150 5,
151 6,
152 6,
153 7,
154 9, // 35
155 9,
156 9,
157 9,
158 9,
159 10, // 40
160 10,
161 12,
162 3,
163 4,
164 4, // 45
165 5,
166 14,
167 3,
168 6,
169 6, // 50
170 6,
171 3,
172 3,
173 3 // 54
174 };
175
176 static const int tiR[55] = {
177 0, // 0
178 0, // 1
179 0,
180 0,
181 0,
182 0, // 5
183 0,
184 0,
185 4, // 8
186 6,
187 12, // 10
188 1,
189 5,
190 4,
191 2,
192 13, // 15
193 9,
194 3,
195 4,
196 11,
197 4, // 20
198 13,
199 1,
200 7,
201 1,
202 9, // 25
203 10,
204 10,
205 3,
206 7,
207 10, // 30
208 10,
209 6,
210 10,
211 10,
212 1, // 35
213 2,
214 3,
215 4,
216 8,
217 6, // 40
218 9,
219 8,
220 16,
221 22,
222 23, // 45
223 23,
224 10,
225 50,
226 44,
227 46, // 50
228 50,
229 0,
230 1,
231 4 // 54
232 };
233
234 static const doublereal ni[57] = {
235 +0.0,
236 +0.12533547935523E-1, // 1
237 +0.78957634722828E1, // 2
238 -0.87803203303561E1, // 3
239 +0.31802509345418E0, // 4
240 -0.26145533859358E0, // 5
241 -0.78199751687981E-2, // 6
242 +0.88089493102134E-2, // 7
243 -0.66856572307965E0, // 8
244 +0.20433810950965, // 9
245 -0.66212605039687E-4, // 10
246 -0.19232721156002E0, // 11
247 -0.25709043003438E0, // 12
248 +0.16074868486251E0, // 13
249 -0.40092828925807E-1, // 14
250 +0.39343422603254E-6, // 15
251 -0.75941377088144E-5, // 16
252 +0.56250979351888E-3, // 17
253 -0.15608652257135E-4, // 18
254 +0.11537996422951E-8, // 19
255 +0.36582165144204E-6, // 20
256 -0.13251180074668E-11,// 21
257 -0.62639586912454E-9, // 22
258 -0.10793600908932E0, // 23
259 +0.17611491008752E-1, // 24
260 +0.22132295167546E0, // 25
261 -0.40247669763528E0, // 26
262 +0.58083399985759E0, // 27
263 +0.49969146990806E-2, // 28
264 -0.31358700712549E-1, // 29
265 -0.74315929710341E0, // 30
266 +0.47807329915480E0, // 31
267 +0.20527940895948E-1, // 32
268 -0.13636435110343E0, // 33
269 +0.14180634400617E-1, // 34
270 +0.83326504880713E-2, // 35
271 -0.29052336009585E-1, // 36
272 +0.38615085574206E-1, // 37
273 -0.20393486513704E-1, // 38
274 -0.16554050063734E-2, // 39
275 +0.19955571979541E-2, // 40
276 +0.15870308324157E-3, // 41
277 -0.16388568342530E-4, // 42
278 +0.43613615723811E-1, // 43
279 +0.34994005463765E-1, // 44
280 -0.76788197844621E-1, // 45
281 +0.22446277332006E-1, // 46
282 -0.62689710414685E-4, // 47
283 -0.55711118565645E-9, // 48
284 -0.19905718354408E0, // 49
285 +0.31777497330738E0, // 50
286 -0.11841182425981E0, // 51
287 -0.31306260323435E2, // 52
288 +0.31546140237781E2, // 53
289 -0.25213154341695E4, // 54
290 -0.14874640856724E0, // 55
291 +0.31806110878444E0 // 56
292 };
293
294 static const doublereal alphai[3] = {
295 +20.,
296 +20.,
297 +20.
298 };
299
300 static const doublereal betai[3] = {
301 +150.,
302 +150.,
303 +250.
304 };
305
306 static const doublereal gammai[3] = {
307 +1.21,
308 +1.21,
309 +1.25
310 };
311
312 static const doublereal epsi[3] = {
313 +1.0,
314 +1.0,
315 +1.0
316 };
317
318 static const doublereal ai[2] = {
319 +3.5,
320 +3.5
321 };
322
323 static const doublereal bi[2] = {
324 +0.85,
325 +0.95
326 };
327
328 static const doublereal Bi[2] = {
329 +0.2,
330 +0.2
331 };
332
333 static const doublereal Ci[2] = {
334 +28.0,
335 +32.0
336 };
337
338 static const doublereal Di[2] = {
339 +700.,
340 +800.
341 };
342
343 static const doublereal Ai[2] = {
344 +0.32,
345 +0.32
346 };
347
348 static const doublereal Bbetai[2] = {
349 +0.3,
350 +0.3
351 };
352 // \endcond
353
WaterPropsIAPWSphi()354 WaterPropsIAPWSphi::WaterPropsIAPWSphi() :
355 TAUsave(-1.0),
356 TAUsqrt(-1.0),
357 DELTAsave(-1.0)
358 {
359 for (int i = 0; i < 52; i++) {
360 TAUp[i] = 1.0;
361 }
362 for (int i = 0; i < 16; i++) {
363 DELTAp[i] = 1.0;
364 }
365 }
366
tdpolycalc(doublereal tau,doublereal delta)367 void WaterPropsIAPWSphi::tdpolycalc(doublereal tau, doublereal delta)
368 {
369 if ((tau != TAUsave) || 1) {
370 TAUsave = tau;
371 TAUsqrt = sqrt(tau);
372 TAUp[0] = 1.0;
373 for (int i = 1; i < 51; i++) {
374 TAUp[i] = TAUp[i-1] * tau;
375 }
376 }
377 if ((delta != DELTAsave) || 1) {
378 DELTAsave = delta;
379 DELTAp[0] = 1.0;
380 for (int i = 1; i <= 15; i++) {
381 DELTAp[i] = DELTAp[i-1] * delta;
382 }
383 }
384 }
385
phi0() const386 doublereal WaterPropsIAPWSphi::phi0() const
387 {
388 doublereal tau = TAUsave;
389 doublereal delta = DELTAsave;
390 doublereal retn = log(delta) + ni0[1] + ni0[2]*tau + ni0[3]*log(tau);
391
392 retn += ni0[4] * log(1.0 - exp(-gammi0[4]*tau));
393 retn += ni0[5] * log(1.0 - exp(-gammi0[5]*tau));
394 retn += ni0[6] * log(1.0 - exp(-gammi0[6]*tau));
395 retn += ni0[7] * log(1.0 - exp(-gammi0[7]*tau));
396 retn += ni0[8] * log(1.0 - exp(-gammi0[8]*tau));
397 return retn;
398 }
399
phiR() const400 doublereal WaterPropsIAPWSphi::phiR() const
401 {
402 doublereal tau = TAUsave;
403 doublereal delta = DELTAsave;
404 int i, j;
405
406 // Write out the first seven polynomials in the expression
407 doublereal T375 = pow(tau, 0.375);
408 doublereal val = (ni[1] * delta / TAUsqrt +
409 ni[2] * delta * TAUsqrt * T375 +
410 ni[3] * delta * tau +
411 ni[4] * DELTAp[2] * TAUsqrt +
412 ni[5] * DELTAp[2] * T375 * T375 +
413 ni[6] * DELTAp[3] * T375 +
414 ni[7] * DELTAp[4] * tau);
415 // Next, do polynomial contributions 8 to 51
416 for (i = 8; i <= 51; i++) {
417 val += (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * exp(-DELTAp[ciR[i]]));
418 }
419
420 // Next do contributions 52 to 54
421 for (j = 0; j < 3; j++) {
422 i = 52 + j;
423 doublereal dtmp = delta - epsi[j];
424 doublereal ttmp = tau - gammai[j];
425 val += (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] *
426 exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
427 }
428
429 // Next do contributions 55 and 56
430 for (j = 0; j < 2; j++) {
431 i = 55 + j;
432 doublereal deltam1 = delta - 1.0;
433 doublereal dtmp2 = deltam1 * deltam1;
434 doublereal atmp = 0.5 / Bbetai[j];
435 doublereal theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
436 doublereal triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
437 doublereal ttmp = tau - 1.0;
438 doublereal triagtmp = pow(triag, bi[j]);
439 doublereal phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
440 val += (ni[i] * triagtmp * delta * phi);
441 }
442
443 return val;
444 }
445
phi(doublereal tau,doublereal delta)446 doublereal WaterPropsIAPWSphi::phi(doublereal tau, doublereal delta)
447 {
448 tdpolycalc(tau, delta);
449 doublereal nau = phi0();
450 doublereal res = phiR();
451 return nau + res;
452 }
453
phiR_d() const454 doublereal WaterPropsIAPWSphi::phiR_d() const
455 {
456 doublereal tau = TAUsave;
457 doublereal delta = DELTAsave;
458 int i, j;
459
460 // Write out the first seven polynomials in the expression
461 doublereal T375 = pow(tau, 0.375);
462 doublereal val = (ni[1] / TAUsqrt +
463 ni[2] * TAUsqrt * T375 +
464 ni[3] * tau +
465 ni[4] * 2.0 * delta * TAUsqrt +
466 ni[5] * 2.0 * delta * T375 * T375 +
467 ni[6] * 3.0 * DELTAp[2] * T375 +
468 ni[7] * 4.0 * DELTAp[3] * tau);
469 // Next, do polynomial contributions 8 to 51
470 for (i = 8; i <= 51; i++) {
471 val += ((ni[i] * exp(-DELTAp[ciR[i]]) * DELTAp[diR[i] - 1] *
472 TAUp[tiR[i]]) * (diR[i] - ciR[i]* DELTAp[ciR[i]]));
473 }
474
475 // Next do contributions 52 to 54
476 for (j = 0; j < 3; j++) {
477 i = 52 + j;
478 doublereal dtmp = delta - epsi[j];
479 doublereal ttmp = tau - gammai[j];
480 doublereal tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] *
481 exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
482 val += tmp * (diR[i]/delta - 2.0 * alphai[j] * dtmp);
483 }
484
485 // Next do contributions 55 and 56
486 for (j = 0; j < 2; j++) {
487 i = 55 + j;
488 doublereal deltam1 = delta - 1.0;
489 doublereal dtmp2 = deltam1 * deltam1;
490 doublereal atmp = 0.5 / Bbetai[j];
491 doublereal theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
492 doublereal triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
493 doublereal ttmp = tau - 1.0;
494 doublereal triagtmp = pow(triag, bi[j]);
495 doublereal triagtmpm1 = pow(triag, bi[j]-1.0);
496 doublereal atmpM1 = atmp - 1.0;
497 doublereal ptmp = pow(dtmp2,atmpM1);
498 doublereal p2tmp = pow(dtmp2, ai[j]-1.0);
499 doublereal dtriagddelta =
500 deltam1 *(Ai[j] * theta * 2.0 / Bbetai[j] * ptmp +
501 2.0*Bi[j]*ai[j]*p2tmp);
502 doublereal phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
503 doublereal dphiddelta = -2.0*Ci[j]*deltam1*phi;
504 doublereal dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta;
505 doublereal tmp = ni[i] * (triagtmp * (phi + delta*dphiddelta) +
506 dtriagtmpddelta * delta * phi);
507 val += tmp;
508 }
509 return val;
510 }
511
phi0_d() const512 doublereal WaterPropsIAPWSphi::phi0_d() const
513 {
514 doublereal delta = DELTAsave;
515 return 1.0/delta;
516 }
517
phi_d(doublereal tau,doublereal delta)518 doublereal WaterPropsIAPWSphi::phi_d(doublereal tau, doublereal delta)
519 {
520 tdpolycalc(tau, delta);
521 doublereal nau = phi0_d();
522 doublereal res = phiR_d();
523 return nau + res;
524 }
525
pressureM_rhoRT(doublereal tau,doublereal delta)526 doublereal WaterPropsIAPWSphi::pressureM_rhoRT(doublereal tau, doublereal delta)
527 {
528 tdpolycalc(tau, delta);
529 doublereal res = phiR_d();
530 return 1.0 + delta * res;
531 }
532
phiR_dd() const533 doublereal WaterPropsIAPWSphi::phiR_dd() const
534 {
535 doublereal tau = TAUsave;
536 doublereal delta = DELTAsave;
537 int i, j;
538 doublereal atmp;
539
540 // Write out the first seven polynomials in the expression
541 doublereal T375 = pow(tau, 0.375);
542 doublereal val = (ni[4] * 2.0 * TAUsqrt +
543 ni[5] * 2.0 * T375 * T375 +
544 ni[6] * 6.0 * delta * T375 +
545 ni[7] * 12.0 * DELTAp[2] * tau);
546 // Next, do polynomial contributions 8 to 51
547 for (i = 8; i <= 51; i++) {
548 doublereal dtmp = DELTAp[ciR[i]];
549 doublereal tmp = ni[i] * exp(-dtmp) * TAUp[tiR[i]];
550 if (diR[i] == 1) {
551 atmp = 1.0/delta;
552 } else {
553 atmp = DELTAp[diR[i] - 2];
554 }
555 tmp *= atmp *((diR[i] - ciR[i]*dtmp)*(diR[i]-1.0-ciR[i]*dtmp) -
556 ciR[i]*ciR[i]*dtmp);
557 val += tmp;
558 }
559
560 // Next do contributions 52 to 54
561 for (j = 0; j < 3; j++) {
562 i = 52 + j;
563 doublereal dtmp = delta - epsi[j];
564 doublereal ttmp = tau - gammai[j];
565 doublereal tmp = (ni[i] * TAUp[tiR[i]] *
566 exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
567 doublereal deltmp = DELTAp[diR[i]];
568 doublereal deltmpM1 = deltmp/delta;
569 doublereal deltmpM2 = deltmpM1 / delta;
570 doublereal d2tmp = dtmp * dtmp;
571
572 val += tmp * (-2.0*alphai[j]*deltmp +
573 4.0 * alphai[j] * alphai[j] * deltmp * d2tmp -
574 4.0 * diR[i] * alphai[j] * deltmpM1 * dtmp +
575 diR[i] * (diR[i] - 1.0) * deltmpM2);
576 }
577
578 // Next do contributions 55 and 56
579 for (j = 0; j < 2; j++) {
580 i = 55 + j;
581 doublereal deltam1 = delta - 1.0;
582 doublereal dtmp2 = deltam1 * deltam1;
583 atmp = 0.5 / Bbetai[j];
584 doublereal theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
585 doublereal triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
586 doublereal ttmp = tau - 1.0;
587 doublereal triagtmp = pow(triag, bi[j]);
588 doublereal triagtmpm1 = pow(triag, bi[j]-1.0);
589 doublereal atmpM1 = atmp - 1.0;
590 doublereal ptmp = pow(dtmp2,atmpM1);
591 doublereal p2tmp = pow(dtmp2, ai[j]-1.0);
592 doublereal dtriagddelta =
593 deltam1 *(Ai[j] * theta * 2.0 / Bbetai[j] * ptmp +
594 2.0*Bi[j]*ai[j]*p2tmp);
595 doublereal phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
596 doublereal dphiddelta = -2.0*Ci[j]*deltam1*phi;
597 doublereal dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta;
598 doublereal d2phiddelta2 = 2.0 * Ci[j] * phi * (2.0*Ci[j]*dtmp2 - 1.0);
599 doublereal pptmp = ptmp / dtmp2;
600 doublereal d2triagddelta2 = dtriagddelta / deltam1;
601 d2triagddelta2 +=
602 dtmp2 *(4.0*Bi[j]*ai[j]*(ai[j]-1.0)*pow(dtmp2,ai[j]-2.0) +
603 2.0*Ai[j]*Ai[j]/(Bbetai[j]*Bbetai[j])*ptmp*ptmp +
604 Ai[j]*theta*4.0/Bbetai[j]*(atmp-1.0)*pptmp);
605 doublereal d2triagtmpd2delta =
606 bi[j] * (triagtmpm1 * d2triagddelta2 +
607 (bi[j]-1.0)*triagtmpm1/triag*dtriagddelta*dtriagddelta);
608 doublereal ctmp = (triagtmp * (2.0*dphiddelta + delta*d2phiddelta2) +
609 2.0*dtriagtmpddelta*(phi + delta * dphiddelta) +
610 d2triagtmpd2delta * delta * phi);
611 val += ni[i] * ctmp;
612 }
613 return val;
614 }
615
phi0_dd() const616 doublereal WaterPropsIAPWSphi::phi0_dd() const
617 {
618 doublereal delta = DELTAsave;
619 return -1.0/(delta*delta);
620 }
621
phi_dd(doublereal tau,doublereal delta)622 doublereal WaterPropsIAPWSphi::phi_dd(doublereal tau, doublereal delta)
623 {
624 tdpolycalc(tau, delta);
625 doublereal nau = phi0_dd();
626 doublereal res = phiR_dd();
627 return nau + res;
628 }
629
dimdpdrho(doublereal tau,doublereal delta)630 doublereal WaterPropsIAPWSphi::dimdpdrho(doublereal tau, doublereal delta)
631 {
632 tdpolycalc(tau, delta);
633 doublereal res1 = phiR_d();
634 doublereal res2 = phiR_dd();
635 return 1.0 + delta * (2.0*res1 + delta*res2);
636 }
637
dimdpdT(doublereal tau,doublereal delta)638 doublereal WaterPropsIAPWSphi::dimdpdT(doublereal tau, doublereal delta)
639 {
640 tdpolycalc(tau, delta);
641 doublereal res1 = phiR_d();
642 doublereal res2 = phiR_dt();
643 return (1.0 + delta * res1) - tau * delta * (res2);
644 }
645
phi0_t() const646 doublereal WaterPropsIAPWSphi::phi0_t() const
647 {
648 doublereal tau = TAUsave;
649 doublereal retn = ni0[2] + ni0[3]/tau;
650 retn += (ni0[4] * gammi0[4] * (1.0/(1.0 - exp(-gammi0[4]*tau)) - 1.0));
651 retn += (ni0[5] * gammi0[5] * (1.0/(1.0 - exp(-gammi0[5]*tau)) - 1.0));
652 retn += (ni0[6] * gammi0[6] * (1.0/(1.0 - exp(-gammi0[6]*tau)) - 1.0));
653 retn += (ni0[7] * gammi0[7] * (1.0/(1.0 - exp(-gammi0[7]*tau)) - 1.0));
654 retn += (ni0[8] * gammi0[8] * (1.0/(1.0 - exp(-gammi0[8]*tau)) - 1.0));
655 return retn;
656 }
657
phiR_t() const658 doublereal WaterPropsIAPWSphi::phiR_t() const
659 {
660 doublereal tau = TAUsave;
661 doublereal delta = DELTAsave;
662 int i, j;
663 doublereal atmp, tmp;
664
665 // Write out the first seven polynomials in the expression
666 doublereal T375 = pow(tau, 0.375);
667 doublereal val = ((-0.5) *ni[1] * delta / TAUsqrt / tau +
668 ni[2] * delta * 0.875 / TAUsqrt * T375 +
669 ni[3] * delta +
670 ni[4] * DELTAp[2] * 0.5 / TAUsqrt +
671 ni[5] * DELTAp[2] * 0.75 * T375 * T375 / tau +
672 ni[6] * DELTAp[3] * 0.375 * T375 / tau +
673 ni[7] * DELTAp[4]);
674 // Next, do polynomial contributions 8 to 51
675 for (i = 8; i <= 51; i++) {
676 tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]-1] * exp(-DELTAp[ciR[i]]));
677 val += tiR[i] * tmp;
678 }
679
680 // Next do contributions 52 to 54
681 for (j = 0; j < 3; j++) {
682 i = 52 + j;
683 doublereal dtmp = delta - epsi[j];
684 doublereal ttmp = tau - gammai[j];
685 tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] *
686 exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
687 val += tmp *(tiR[i]/tau - 2.0 * betai[j]*ttmp);
688 }
689
690 // Next do contributions 55 and 56
691 for (j = 0; j < 2; j++) {
692 i = 55 + j;
693 doublereal deltam1 = delta - 1.0;
694 doublereal dtmp2 = deltam1 * deltam1;
695 atmp = 0.5 / Bbetai[j];
696 doublereal theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
697 doublereal triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
698 doublereal ttmp = tau - 1.0;
699 doublereal triagtmp = pow(triag, bi[j]);
700 doublereal phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
701 doublereal dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag;
702 doublereal dphidtau = - 2.0 * Di[j] * ttmp * phi;
703 val += ni[i] * delta * (dtriagtmpdtau * phi + triagtmp * dphidtau);
704 }
705 return val;
706 }
707
phi_t(doublereal tau,doublereal delta)708 doublereal WaterPropsIAPWSphi::phi_t(doublereal tau, doublereal delta)
709 {
710 tdpolycalc(tau, delta);
711 doublereal nau = phi0_t();
712 doublereal res = phiR_t();
713 return nau + res;
714 }
715
phi0_tt() const716 doublereal WaterPropsIAPWSphi::phi0_tt() const
717 {
718 doublereal tau = TAUsave;
719 doublereal tmp, itmp;
720 doublereal retn = - ni0[3]/(tau * tau);
721 for (int i = 4; i <= 8; i++) {
722 tmp = exp(-gammi0[i]*tau);
723 itmp = 1.0 - tmp;
724 retn -= (ni0[i] * gammi0[i] * gammi0[i] * tmp / (itmp * itmp));
725 }
726 return retn;
727 }
728
phiR_tt() const729 doublereal WaterPropsIAPWSphi::phiR_tt() const
730 {
731 doublereal tau = TAUsave;
732 doublereal delta = DELTAsave;
733 int i, j;
734 doublereal atmp, tmp;
735
736 // Write out the first seven polynomials in the expression
737 doublereal T375 = pow(tau, 0.375);
738 doublereal val = ((-0.5) * (-1.5) * ni[1] * delta / (TAUsqrt * tau * tau) +
739 ni[2] * delta * 0.875 * (-0.125) * T375 / (TAUsqrt * tau) +
740 ni[4] * DELTAp[2] * 0.5 * (-0.5)/ (TAUsqrt * tau) +
741 ni[5] * DELTAp[2] * 0.75 *(-0.25) * T375 * T375 / (tau * tau) +
742 ni[6] * DELTAp[3] * 0.375 *(-0.625) * T375 / (tau * tau));
743 // Next, do polynomial contributions 8 to 51
744 for (i = 8; i <= 51; i++) {
745 if (tiR[i] > 1) {
746 tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]-2] * exp(-DELTAp[ciR[i]]));
747 val += tiR[i] * (tiR[i] - 1.0) * tmp;
748 }
749 }
750
751 // Next do contributions 52 to 54
752 for (j = 0; j < 3; j++) {
753 i = 52 + j;
754 doublereal dtmp = delta - epsi[j];
755 doublereal ttmp = tau - gammai[j];
756 tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] *
757 exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
758 atmp = tiR[i]/tau - 2.0 * betai[j]*ttmp;
759 val += tmp *(atmp * atmp - tiR[i]/(tau*tau) - 2.0*betai[j]);
760 }
761
762 // Next do contributions 55 and 56
763 for (j = 0; j < 2; j++) {
764 i = 55 + j;
765 doublereal deltam1 = delta - 1.0;
766 doublereal dtmp2 = deltam1 * deltam1;
767 atmp = 0.5 / Bbetai[j];
768 doublereal theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
769 doublereal triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
770 doublereal ttmp = tau - 1.0;
771 doublereal triagtmp = pow(triag, bi[j]);
772 doublereal triagtmpM1 = triagtmp / triag;
773 doublereal phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
774 doublereal dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag;
775 doublereal dphidtau = - 2.0 * Di[j] * ttmp * phi;
776 doublereal d2triagtmpdtau2 =
777 (2 * bi[j] * triagtmpM1 +
778 4 * theta * theta * bi[j] * (bi[j]-1.0) * triagtmpM1 / triag);
779 doublereal d2phidtau2 = 2.0*Di[j]*phi *(2.0*Di[j]*ttmp*ttmp - 1.0);
780 tmp = (d2triagtmpdtau2 * phi +
781 2 * dtriagtmpdtau * dphidtau +
782 triagtmp * d2phidtau2);
783 val += ni[i] * delta * tmp;
784 }
785
786 return val;
787 }
788
phi_tt(doublereal tau,doublereal delta)789 doublereal WaterPropsIAPWSphi::phi_tt(doublereal tau, doublereal delta)
790 {
791 tdpolycalc(tau, delta);
792 doublereal nau = phi0_tt();
793 doublereal res = phiR_tt();
794 return nau + res;
795 }
796
phi0_dt() const797 doublereal WaterPropsIAPWSphi::phi0_dt() const
798 {
799 return 0.0;
800 }
801
phiR_dt() const802 doublereal WaterPropsIAPWSphi::phiR_dt() const
803 {
804 doublereal tau = TAUsave;
805 doublereal delta = DELTAsave;
806 int i, j;
807 doublereal tmp;
808 // Write out the first seven polynomials in the expression
809 doublereal T375 = pow(tau, 0.375);
810 doublereal val = (ni[1] * (-0.5) / (TAUsqrt * tau) +
811 ni[2] * (0.875) * T375 / TAUsqrt +
812 ni[3] +
813 ni[4] * 2.0 * delta * (0.5) / TAUsqrt +
814 ni[5] * 2.0 * delta * (0.75) * T375 * T375 / tau +
815 ni[6] * 3.0 * DELTAp[2] * 0.375 * T375 / tau +
816 ni[7] * 4.0 * DELTAp[3]);
817 // Next, do polynomial contributions 8 to 51
818 for (i = 8; i <= 51; i++) {
819 tmp = (ni[i] * tiR[i] * exp(-DELTAp[ciR[i]]) * DELTAp[diR[i] - 1] *
820 TAUp[tiR[i] - 1]);
821 val += tmp * (diR[i] - ciR[i] * DELTAp[ciR[i]]);
822 }
823
824 // Next do contributions 52 to 54
825 for (j = 0; j < 3; j++) {
826 i = 52 + j;
827 doublereal dtmp = delta - epsi[j];
828 doublereal ttmp = tau - gammai[j];
829 tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] *
830 exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
831 val += tmp * ((diR[i]/delta - 2.0 * alphai[j] * dtmp) *
832 (tiR[i]/tau - 2.0 * betai[j] * ttmp));
833 }
834
835 // Next do contributions 55 and 56
836 for (j = 0; j < 2; j++) {
837 i = 55 + j;
838 doublereal deltam1 = delta - 1.0;
839 doublereal dtmp2 = deltam1 * deltam1;
840 doublereal atmp = 0.5 / Bbetai[j];
841 doublereal theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
842 doublereal triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
843 doublereal ttmp = tau - 1.0;
844 doublereal triagtmp = pow(triag, bi[j]);
845 doublereal triagtmpm1 = pow(triag, bi[j]-1.0);
846 doublereal atmpM1 = atmp - 1.0;
847 doublereal ptmp = pow(dtmp2,atmpM1);
848 doublereal p2tmp = pow(dtmp2, ai[j]-1.0);
849 doublereal dtriagddelta =
850 deltam1 *(Ai[j] * theta * 2.0 / Bbetai[j] * ptmp +
851 2.0*Bi[j]*ai[j]*p2tmp);
852 doublereal phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
853 doublereal dphiddelta = -2.0*Ci[j]*deltam1*phi;
854 doublereal dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta;
855 doublereal dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag;
856 doublereal dphidtau = - 2.0 * Di[j] * ttmp * phi;
857 doublereal d2phiddeltadtau = 4.0 * Ci[j] * Di[j] * deltam1 * ttmp * phi;
858 doublereal d2triagtmpddeltadtau =
859 (-Ai[j] * bi[j] * 2.0 / Bbetai[j] * triagtmpm1 * deltam1 * ptmp
860 -2.0 * theta * bi[j] * (bi[j] - 1.0) * triagtmpm1 / triag * dtriagddelta);
861 doublereal tmp = ni[i] * (triagtmp * (dphidtau + delta*d2phiddeltadtau) +
862 delta * dtriagtmpddelta * dphidtau +
863 dtriagtmpdtau * (phi + delta * dphiddelta) +
864 d2triagtmpddeltadtau * delta * phi);
865 val += tmp;
866 }
867 return val;
868 }
869
dfind(doublereal p_red,doublereal tau,doublereal deltaGuess)870 doublereal WaterPropsIAPWSphi::dfind(doublereal p_red, doublereal tau, doublereal deltaGuess)
871 {
872 doublereal dd = deltaGuess;
873 bool conv = false;
874 doublereal deldd = dd;
875 doublereal pcheck = 1.0E-30 + 1.0E-8 * p_red;
876 for (int n = 0; n < 200; n++) {
877
878 // Calculate the internal polynomials, and then calculate the phi deriv
879 // functions needed by this routine.
880 tdpolycalc(tau, dd);
881 doublereal q1 = phiR_d();
882 doublereal q2 = phiR_dd();
883
884 // Calculate the predicted reduced pressure, pred0, based on the current
885 // tau and dd.
886 doublereal pred0 = dd + dd * dd * q1;
887
888 // Calculate the derivative of the predicted reduced pressure wrt the
889 // reduced density, dd, This is dpddelta
890 doublereal dpddelta = 1.0 + 2.0 * dd * q1 + dd * dd * q2;
891
892 // If dpddelta is negative, then we are in the middle of the 2 phase
893 // region, beyond the stability curve. We need to adjust the initial
894 // guess outwards and start a new iteration.
895 if (dpddelta <= 0.0) {
896 if (deltaGuess > 1.0) {
897 dd = dd * 1.05;
898 }
899 if (deltaGuess < 1.0) {
900 dd = dd * 0.95;
901 }
902 continue;
903 }
904
905 // Check for convergence
906 if (fabs(pred0-p_red) < pcheck) {
907 conv = true;
908 break;
909 }
910
911 // Dampen and crop the update
912 doublereal dpdx = dpddelta;
913 if (n < 10) {
914 dpdx = dpddelta * 1.1;
915 }
916 dpdx = std::max(dpdx, 0.001);
917
918 // Formulate the update to reduced density using Newton's method. Then,
919 // crop it to a max value of 0.02
920 deldd = - (pred0 - p_red) / dpdx;
921 if (fabs(deldd) > 0.05) {
922 deldd = deldd * 0.05 / fabs(deldd);
923 }
924
925 // updated the reduced density value
926 dd += deldd;
927 if (fabs(deldd/dd) < 1.0E-14) {
928 conv = true;
929 break;
930 }
931
932 // Check for negative densities
933 if (dd <= 0.0) {
934 dd = 1.0E-24;
935 }
936 }
937
938 // Check for convergence, and return 0.0 if it wasn't achieved.
939 if (! conv) {
940 dd = 0.0;
941 }
942 return dd;
943 }
944
gibbs_RT() const945 doublereal WaterPropsIAPWSphi::gibbs_RT() const
946 {
947 doublereal delta = DELTAsave;
948 doublereal rd = phiR_d();
949 return 1.0 + phi0() + phiR() + delta * rd;
950 }
951
enthalpy_RT() const952 doublereal WaterPropsIAPWSphi::enthalpy_RT() const
953 {
954 doublereal delta = DELTAsave;
955 doublereal tau = TAUsave;
956 doublereal rd = phiR_d();
957 doublereal nt = phi0_t();
958 doublereal rt = phiR_t();
959 return 1.0 + tau * (nt + rt) + delta * rd;
960 }
961
entropy_R() const962 doublereal WaterPropsIAPWSphi::entropy_R() const
963 {
964 doublereal tau = TAUsave;
965 doublereal nt = phi0_t();
966 doublereal rt = phiR_t();
967 doublereal p0 = phi0();
968 doublereal pR = phiR();
969 return tau * (nt + rt) - p0 - pR;
970 }
971
intEnergy_RT() const972 doublereal WaterPropsIAPWSphi::intEnergy_RT() const
973 {
974 doublereal tau = TAUsave;
975 doublereal nt = phi0_t();
976 doublereal rt = phiR_t();
977 return tau * (nt + rt);
978 }
979
cv_R() const980 doublereal WaterPropsIAPWSphi::cv_R() const
981 {
982 doublereal tau = TAUsave;
983 doublereal ntt = phi0_tt();
984 doublereal rtt = phiR_tt();
985 return - tau * tau * (ntt + rtt);
986 }
987
cp_R() const988 doublereal WaterPropsIAPWSphi::cp_R() const
989 {
990 doublereal tau = TAUsave;
991 doublereal delta = DELTAsave;
992 doublereal cvR = cv_R();
993 doublereal rd = phiR_d();
994 doublereal rdd = phiR_dd();
995 doublereal rdt = phiR_dt();
996 doublereal num = (1.0 + delta * rd - delta * tau * rdt);
997 doublereal cpR = cvR + (num * num /
998 (1.0 + 2.0 * delta * rd + delta * delta * rdd));
999 return cpR;
1000 }
1001
1002 } // namespace Cantera
1003