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