Cantera  2.4.0
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 http://www.cantera.org/license.txt for license and copyright information.
10 
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 
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 
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 
386 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 
400 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 
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 
454 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 
512 doublereal WaterPropsIAPWSphi::phi0_d() const
513 {
514  doublereal delta = DELTAsave;
515  return 1.0/delta;
516 }
517 
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 
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 
533 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 
616 doublereal WaterPropsIAPWSphi::phi0_dd() const
617 {
618  doublereal delta = DELTAsave;
619  return -1.0/(delta*delta);
620 }
621 
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 
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 
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 
646 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 
658 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 
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 
716 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 
729 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 
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 
797 doublereal WaterPropsIAPWSphi::phi0_dt() const
798 {
799  return 0.0;
800 }
801 
802 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 
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 
946 {
947  doublereal delta = DELTAsave;
948  doublereal rd = phiR_d();
949  return 1.0 + phi0() + phiR() + delta * rd;
950 }
951 
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 
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 
973 {
974  doublereal tau = TAUsave;
975  doublereal nt = phi0_t();
976  doublereal rt = phiR_t();
977  return tau * (nt + rt);
978 }
979 
980 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 
988 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
doublereal phiR_d() const
Calculate d_phiR_d(delta), the first derivative of phiR wrt delta.
Header for Lowest level of the classes which support a real water model (see class WaterPropsIAPWS an...
doublereal phi_dd(doublereal tau, doublereal delta)
2nd derivative of phi wrt delta
doublereal phiR_t() const
Calculate Equation 6.6 for dphiRdtau, the derivative residual part of the dimensionless Helmholtz fre...
This file contains definitions for utility functions and text for modules, inputfiles, logs, textlogs, (see Input File Handling, Diagnostic Output, and Writing messages to the screen).
doublereal dfind(doublereal p_red, doublereal tau, doublereal deltaGuess)
This function computes the reduced density, given the reduced pressure and the reduced temperature...
doublereal TAUsave
Last tau that was used to calculate polynomials.
doublereal phi0_dt() const
Calculate the mixed derivative d2_phi0/(dtau ddelta)
doublereal DELTAp[16]
Value of internally calculated polynomials of powers of delta.
doublereal gibbs_RT() const
Calculate the dimensionless Gibbs free energy.
doublereal DELTAsave
Last delta that was used to calculate polynomials.
doublereal pressureM_rhoRT(doublereal tau, doublereal delta)
Calculate the dimensionless pressure at tau and delta;.
doublereal cp_R() const
Calculate the dimensionless constant pressure heat capacity, Cv/R.
doublereal phi0_tt() const
Calculate d2_phi0/dtau2.
doublereal phi0_t() const
Calculate d_phi0/d(tau)
doublereal dimdpdrho(doublereal tau, doublereal delta)
Dimensionless derivative of p wrt rho at constant T.
doublereal phi_tt(doublereal tau, doublereal delta)
Second derivative of phi wrt tau.
doublereal phi0_d() const
Calculate d_phi0_d(delta), the first derivative of phi0 wrt delta.
doublereal dimdpdT(doublereal tau, doublereal delta)
Dimensionless derivative of p wrt T at constant rho.
doublereal cv_R() const
Calculate the dimensionless constant volume heat capacity, Cv/R.
doublereal phiR_dd() const
Calculate d2_phiR_dd(delta), the second derivative of phiR wrt delta.
WaterPropsIAPWSphi()
Base constructor.
void tdpolycalc(doublereal tau, doublereal delta)
Calculates internal polynomials in tau and delta.
doublereal phi_t(doublereal tau, doublereal delta)
First derivative of phi wrt tau.
doublereal phi0_dd() const
Calculate d2_phi0_dd(delta), the second derivative of phi0 wrt delta.
doublereal phiR_tt() const
Calculate Equation 6.6 for dphiRdtau, the second derivative residual part of the dimensionless Helmho...
doublereal TAUsqrt
sqrt of TAU
doublereal intEnergy_RT() const
Calculate the dimensionless internal energy, u/RT.
doublereal phi0() const
Calculate Equation 6.5 for phi0, the ideal gas part of the dimensionless Helmholtz free energy...
doublereal phi_d(doublereal tau, doublereal delta)
Calculate derivative of phi wrt delta.
doublereal phi(doublereal tau, doublereal delta)
Calculate the Phi function, which is the base function.
doublereal phiR_dt() const
Calculate the mixed derivative d2_phiR/(dtau ddelta)
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
doublereal entropy_R() const
Calculate the dimensionless entropy, s/R.
doublereal TAUp[52]
Value of internally calculated polynomials of powers of TAU.
doublereal enthalpy_RT() const
Calculate the dimensionless enthalpy, h/RT.