www.mooseframework.org
Water97FluidProperties.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
9 
10 template <>
11 InputParameters
13 {
14  InputParameters params = validParams<SinglePhaseFluidPropertiesPT>();
15  params.addClassDescription("Fluid properties for water and steam (H2O) using IAPWS-IF97");
16  return params;
17 }
18 
19 Water97FluidProperties::Water97FluidProperties(const InputParameters & parameters)
20  : SinglePhaseFluidPropertiesPT(parameters),
21  _Mh2o(18.015e-3),
22  _Rw(461.526),
23  _p_critical(22.064e6),
24  _T_critical(647.096),
25  _rho_critical(322.0),
26  _p_triple(611.657),
27  _T_triple(273.16)
28 {
29 }
30 
32 
33 std::string
35 {
36  return "water";
37 }
38 
39 Real
41 {
42  return _Mh2o;
43 }
44 
45 Real
47 {
48  return _p_critical;
49 }
50 
51 Real
53 {
54  return _T_critical;
55 }
56 
57 Real
59 {
60  return _rho_critical;
61 }
62 
63 Real
65 {
66  return _p_triple;
67 }
68 
69 Real
71 {
72  return _T_triple;
73 }
74 
75 Real
77 {
78  Real density, pi, tau;
79 
80  // Determine which region the point is in
81  unsigned int region = inRegion(pressure, temperature);
82 
83  switch (region)
84  {
85  case 1:
86  pi = pressure / _p_star[0];
87  tau = _T_star[0] / temperature;
88  density = pressure / (pi * _Rw * temperature * dgamma1_dpi(pi, tau));
89  break;
90 
91  case 2:
92  pi = pressure / _p_star[1];
93  tau = _T_star[1] / temperature;
94  density = pressure / (pi * _Rw * temperature * dgamma2_dpi(pi, tau));
95  break;
96 
97  case 3:
98  density = densityRegion3(pressure, temperature);
99  break;
100 
101  case 5:
102  pi = pressure / _p_star[4];
103  tau = _T_star[4] / temperature;
104  density = pressure / (pi * _Rw * temperature * dgamma5_dpi(pi, tau));
105  break;
106 
107  default:
108  mooseError(name(), ": inRegion() has given an incorrect region");
109  }
110  return density;
111 }
112 
113 void
115  Real pressure, Real temperature, Real & rho, Real & drho_dp, Real & drho_dT) const
116 {
117  Real pi, tau, ddensity_dp, ddensity_dT;
118 
119  // Determine which region the point is in
120  unsigned int region = inRegion(pressure, temperature);
121 
122  switch (region)
123  {
124  case 1:
125  {
126  pi = pressure / _p_star[0];
127  tau = _T_star[0] / temperature;
128  Real dgdp = dgamma1_dpi(pi, tau);
129  ddensity_dp = -d2gamma1_dpi2(pi, tau) / (_Rw * temperature * dgdp * dgdp);
130  ddensity_dT = -pressure * (dgdp - tau * d2gamma1_dpitau(pi, tau)) /
131  (_Rw * pi * temperature * temperature * dgdp * dgdp);
132  break;
133  }
134 
135  case 2:
136  {
137  pi = pressure / _p_star[1];
138  tau = _T_star[1] / temperature;
139  Real dgdp = dgamma2_dpi(pi, tau);
140  ddensity_dp = -d2gamma2_dpi2(pi, tau) / (_Rw * temperature * dgdp * dgdp);
141  ddensity_dT = -pressure * (dgdp - tau * d2gamma2_dpitau(pi, tau)) /
142  (_Rw * pi * temperature * temperature * dgdp * dgdp);
143  break;
144  }
145 
146  case 3:
147  {
148  // Calculate density first, then use that in Helmholtz free energy
149  Real density = densityRegion3(pressure, temperature);
150  Real delta = density / _rho_critical;
151  tau = _T_star[2] / temperature;
152  Real dpdd = dphi3_ddelta(delta, tau);
153  Real d2pdd2 = d2phi3_ddelta2(delta, tau);
154  ddensity_dp = 1.0 / (_Rw * temperature * delta * (2.0 * dpdd + delta * d2pdd2));
155  ddensity_dT = density * (tau * d2phi3_ddeltatau(delta, tau) - dpdd) / temperature /
156  (2.0 * dpdd + delta * d2pdd2);
157  break;
158  }
159 
160  case 5:
161  {
162  pi = pressure / _p_star[4];
163  tau = _T_star[4] / temperature;
164  Real dgdp = dgamma5_dpi(pi, tau);
165  ddensity_dp = -d2gamma5_dpi2(pi, tau) / (_Rw * temperature * dgdp * dgdp);
166  ddensity_dT = -pressure * (dgdp - tau * d2gamma5_dpitau(pi, tau)) /
167  (_Rw * pi * temperature * temperature * dgdp * dgdp);
168  break;
169  }
170 
171  default:
172  mooseError(name(), ": inRegion() has given an incorrect region");
173  }
174 
175  rho = this->rho(pressure, temperature);
176  drho_dp = ddensity_dp;
177  drho_dT = ddensity_dT;
178 }
179 
180 Real
182 {
183  Real internal_energy, pi, tau;
184 
185  // Determine which region the point is in
186  unsigned int region = inRegion(pressure, temperature);
187  switch (region)
188  {
189  case 1:
190  pi = pressure / _p_star[0];
191  tau = _T_star[0] / temperature;
192  internal_energy =
193  _Rw * temperature * (tau * dgamma1_dtau(pi, tau) - pi * dgamma1_dpi(pi, tau));
194  break;
195 
196  case 2:
197  pi = pressure / _p_star[1];
198  tau = _T_star[1] / temperature;
199  internal_energy =
200  _Rw * temperature * (tau * dgamma2_dtau(pi, tau) - pi * dgamma2_dpi(pi, tau));
201  break;
202 
203  case 3:
204  {
205  // Calculate density first, then use that in Helmholtz free energy
206  Real density3 = densityRegion3(pressure, temperature);
207  Real delta = density3 / _rho_critical;
208  tau = _T_star[2] / temperature;
209  internal_energy = _Rw * temperature * tau * dphi3_dtau(delta, tau);
210  break;
211  }
212 
213  case 5:
214  pi = pressure / _p_star[4];
215  tau = _T_star[4] / temperature;
216  internal_energy =
217  _Rw * temperature * (tau * dgamma5_dtau(pi, tau) - pi * dgamma5_dpi(pi, tau));
218  break;
219 
220  default:
221  mooseError(name(), ": inRegion() has given an incorrect region");
222  }
223  // Output in J/kg
224  return internal_energy;
225 }
226 
227 void
229  Real pressure, Real temperature, Real & e, Real & de_dp, Real & de_dT) const
230 {
231  Real pi, tau, dinternal_energy_dp, dinternal_energy_dT;
232 
233  // Determine which region the point is in
234  unsigned int region = inRegion(pressure, temperature);
235  switch (region)
236  {
237  case 1:
238  {
239  pi = pressure / _p_star[0];
240  tau = _T_star[0] / temperature;
241  Real dgdp = dgamma1_dpi(pi, tau);
242  Real d2gdpt = d2gamma1_dpitau(pi, tau);
243  dinternal_energy_dp =
244  _Rw * temperature * (tau * d2gdpt - dgdp - pi * d2gamma1_dpi2(pi, tau)) / _p_star[0];
245  dinternal_energy_dT =
246  _Rw * (pi * tau * d2gdpt - tau * tau * d2gamma1_dtau2(pi, tau) - pi * dgdp);
247  break;
248  }
249 
250  case 2:
251  {
252  pi = pressure / _p_star[1];
253  tau = _T_star[1] / temperature;
254  Real dgdp = dgamma2_dpi(pi, tau);
255  Real d2gdpt = d2gamma2_dpitau(pi, tau);
256  dinternal_energy_dp =
257  _Rw * temperature * (tau * d2gdpt - dgdp - pi * d2gamma2_dpi2(pi, tau)) / _p_star[1];
258  dinternal_energy_dT =
259  _Rw * (pi * tau * d2gdpt - tau * tau * d2gamma2_dtau2(pi, tau) - pi * dgdp);
260  break;
261  }
262 
263  case 3:
264  {
265  // Calculate density first, then use that in Helmholtz free energy
266  Real density3 = densityRegion3(pressure, temperature);
267  Real delta = density3 / _rho_critical;
268  tau = _T_star[2] / temperature;
269  Real dpdd = dphi3_ddelta(delta, tau);
270  Real d2pddt = d2phi3_ddeltatau(delta, tau);
271  Real d2pdd2 = d2phi3_ddelta2(delta, tau);
272  dinternal_energy_dp =
273  _T_star[2] * d2pddt / _rho_critical /
274  (2.0 * temperature * delta * dpdd + temperature * delta * delta * d2pdd2);
275  dinternal_energy_dT =
276  -_Rw * (delta * tau * d2pddt * (dpdd - tau * d2pddt) / (2.0 * dpdd + delta * d2pdd2) +
277  tau * tau * d2phi3_dtau2(delta, tau));
278  break;
279  }
280 
281  case 5:
282  {
283  pi = pressure / _p_star[4];
284  tau = _T_star[4] / temperature;
285  Real dgdp = dgamma5_dpi(pi, tau);
286  Real d2gdpt = d2gamma5_dpitau(pi, tau);
287  dinternal_energy_dp =
288  _Rw * temperature * (tau * d2gdpt - dgdp - pi * d2gamma5_dpi2(pi, tau)) / _p_star[4];
289  dinternal_energy_dT =
290  _Rw * (pi * tau * d2gdpt - tau * tau * d2gamma5_dtau2(pi, tau) - pi * dgdp);
291  break;
292  }
293 
294  default:
295  mooseError(name(), ": inRegion has given an incorrect region");
296  }
297 
298  e = this->e(pressure, temperature);
299  de_dp = dinternal_energy_dp;
300  de_dT = dinternal_energy_dT;
301 }
302 
303 void
305  Real temperature,
306  Real & rho,
307  Real & drho_dp,
308  Real & drho_dT,
309  Real & e,
310  Real & de_dp,
311  Real & de_dT) const
312 {
313  rho_dpT(pressure, temperature, rho, drho_dp, drho_dT);
314  e_dpT(pressure, temperature, e, de_dp, de_dT);
315 }
316 
317 Real
319 {
320  Real speed2, pi, tau, delta;
321 
322  // Determine which region the point is in
323  unsigned int region = inRegion(pressure, temperature);
324  switch (region)
325  {
326  case 1:
327  pi = pressure / _p_star[0];
328  tau = _T_star[0] / temperature;
329  speed2 = _Rw * temperature * std::pow(dgamma1_dpi(pi, tau), 2.0) /
330  (std::pow(dgamma1_dpi(pi, tau) - tau * d2gamma1_dpitau(pi, tau), 2.0) /
331  (tau * tau * d2gamma1_dtau2(pi, tau)) -
332  d2gamma1_dpi2(pi, tau));
333  break;
334 
335  case 2:
336  pi = pressure / _p_star[1];
337  tau = _T_star[1] / temperature;
338  speed2 = _Rw * temperature * std::pow(pi * dgamma2_dpi(pi, tau), 2.0) /
339  ((-pi * pi * d2gamma2_dpi2(pi, tau)) +
340  std::pow(pi * dgamma2_dpi(pi, tau) - tau * pi * d2gamma2_dpitau(pi, tau), 2.0) /
341  (tau * tau * d2gamma2_dtau2(pi, tau)));
342  break;
343 
344  case 3:
345  {
346  // Calculate density first, then use that in Helmholtz free energy
347  Real density3 = densityRegion3(pressure, temperature);
348  delta = density3 / _rho_critical;
349  tau = _T_star[2] / temperature;
350  speed2 =
351  _Rw * temperature *
352  (2.0 * delta * dphi3_ddelta(delta, tau) + delta * delta * d2phi3_ddelta2(delta, tau) -
353  std::pow(delta * dphi3_ddelta(delta, tau) - delta * tau * d2phi3_ddeltatau(delta, tau),
354  2.0) /
355  (tau * tau * d2phi3_dtau2(delta, tau)));
356  break;
357  }
358 
359  case 5:
360  pi = pressure / _p_star[4];
361  tau = _T_star[4] / temperature;
362  speed2 = _Rw * temperature * std::pow(pi * dgamma5_dpi(pi, tau), 2.0) /
363  ((-pi * pi * d2gamma5_dpi2(pi, tau)) +
364  std::pow(pi * dgamma5_dpi(pi, tau) - tau * pi * d2gamma5_dpitau(pi, tau), 2.0) /
365  (tau * tau * d2gamma5_dtau2(pi, tau)));
366  break;
367 
368  default:
369  mooseError(name(), ": inRegion() has given an incorrect region");
370  }
371 
372  return std::sqrt(speed2);
373 }
374 
375 Real
377 {
378  Real specific_heat, pi, tau, delta;
379 
380  // Determine which region the point is in
381  unsigned int region = inRegion(pressure, temperature);
382  switch (region)
383  {
384  case 1:
385  pi = pressure / _p_star[0];
386  tau = _T_star[0] / temperature;
387  specific_heat = -_Rw * tau * tau * d2gamma1_dtau2(pi, tau);
388  break;
389 
390  case 2:
391  pi = pressure / _p_star[1];
392  tau = _T_star[1] / temperature;
393  specific_heat = -_Rw * tau * tau * d2gamma2_dtau2(pi, tau);
394  break;
395 
396  case 3:
397  {
398  // Calculate density first, then use that in Helmholtz free energy
399  Real density3 = densityRegion3(pressure, temperature);
400  delta = density3 / _rho_critical;
401  tau = _T_star[2] / temperature;
402  specific_heat =
403  _Rw *
404  (-tau * tau * d2phi3_dtau2(delta, tau) +
405  (delta * dphi3_ddelta(delta, tau) - delta * tau * d2phi3_ddeltatau(delta, tau)) *
406  (delta * dphi3_ddelta(delta, tau) - delta * tau * d2phi3_ddeltatau(delta, tau)) /
407  (2.0 * delta * dphi3_ddelta(delta, tau) +
408  delta * delta * d2phi3_ddelta2(delta, tau)));
409  break;
410  }
411 
412  case 5:
413  pi = pressure / _p_star[4];
414  tau = _T_star[4] / temperature;
415  specific_heat = -_Rw * tau * tau * d2gamma5_dtau2(pi, tau);
416  break;
417 
418  default:
419  mooseError(name(), ": inRegion() has given an incorrect region");
420  }
421  return specific_heat;
422 }
423 
424 Real
426 {
427  Real specific_heat, pi, tau, delta;
428 
429  // Determine which region the point is in
430  unsigned int region = inRegion(pressure, temperature);
431  switch (region)
432  {
433  case 1:
434  pi = pressure / _p_star[0];
435  tau = _T_star[0] / temperature;
436  specific_heat = _Rw * (-tau * tau * d2gamma1_dtau2(pi, tau) +
437  std::pow(dgamma1_dpi(pi, tau) - tau * d2gamma1_dpitau(pi, tau), 2) /
438  d2gamma1_dpi2(pi, tau));
439  break;
440 
441  case 2:
442  pi = pressure / _p_star[1];
443  tau = _T_star[1] / temperature;
444  specific_heat = _Rw * (-tau * tau * d2gamma2_dtau2(pi, tau) +
445  std::pow(dgamma2_dpi(pi, tau) - tau * d2gamma2_dpitau(pi, tau), 2) /
446  d2gamma2_dpi2(pi, tau));
447  break;
448 
449  case 3:
450  {
451  // Calculate density first, then use that in Helmholtz free energy
452  Real density3 = densityRegion3(pressure, temperature);
453  delta = density3 / _rho_critical;
454  tau = _T_star[2] / temperature;
455  specific_heat = _Rw * (-tau * tau * d2phi3_dtau2(delta, tau));
456  break;
457  }
458 
459  case 5:
460  pi = pressure / _p_star[4];
461  tau = _T_star[4] / temperature;
462  specific_heat = _Rw * (-tau * tau * d2gamma5_dtau2(pi, tau) +
463  std::pow(dgamma5_dpi(pi, tau) - tau * d2gamma5_dpitau(pi, tau), 2) /
464  d2gamma5_dpi2(pi, tau));
465  break;
466 
467  default:
468  mooseError(name(), ": inRegion() has given an incorrect region");
469  }
470  return specific_heat;
471 }
472 
473 Real
475 {
476  Real rho = this->rho(pressure, temperature);
477  return this->mu_from_rho_T(rho, temperature);
478 }
479 
480 void
482  Real pressure, Real temperature, Real & mu, Real & dmu_dp, Real & dmu_dT) const
483 {
484  Real rho, drho_dp, drho_dT;
485  this->rho_dpT(pressure, temperature, rho, drho_dp, drho_dT);
486  Real dmu_drho;
487  this->mu_drhoT_from_rho_T(rho, temperature, drho_dT, mu, dmu_drho, dmu_dT);
488  dmu_dp = dmu_drho * drho_dp;
489 }
490 
491 Real
493 {
494  // Constants from Release on the IAPWS Formulation 2008 for the Viscosity of
495  // Ordinary Water Substance
496  const std::vector<int> I{0, 1, 2, 3, 0, 1, 2, 3, 5, 0, 1, 2, 3, 4, 0, 1, 0, 3, 4, 3, 5};
497  const std::vector<int> J{0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 4, 4, 5, 6, 6};
498  const std::vector<Real> H0{1.67752, 2.20462, 0.6366564, -0.241605};
499  const std::vector<Real> H1{
500  5.20094e-1, 8.50895e-2, -1.08374, -2.89555e-1, 2.22531e-1, 9.99115e-1, 1.88797,
501  1.26613, 1.20573e-1, -2.81378e-1, -9.06851e-1, -7.72479e-1, -4.89837e-1, -2.57040e-1,
502  1.61913e-1, 2.57399e-1, -3.25372e-2, 6.98452e-2, 8.72102e-3, -4.35673e-3, -5.93264e-4};
503 
504  Real mu_star = 1.e-6;
505  Real rhobar = density / _rho_critical;
506  Real Tbar = temperature / _T_critical;
507 
508  // Viscosity in limit of zero density
509  Real sum0 = 0.0;
510  for (std::size_t i = 0; i < H0.size(); ++i)
511  sum0 += H0[i] / std::pow(Tbar, i);
512 
513  Real mu0 = 100.0 * std::sqrt(Tbar) / sum0;
514 
515  // Residual component due to finite density
516  Real sum1 = 0.0;
517  for (std::size_t i = 0; i < H1.size(); ++i)
518  sum1 += std::pow(1.0 / Tbar - 1.0, I[i]) * H1[i] * std::pow(rhobar - 1.0, J[i]);
519 
520  Real mu1 = std::exp(rhobar * sum1);
521 
522  // The water viscosity (in Pa.s) is then given by
523  return mu_star * mu0 * mu1;
524 }
525 
526 void
528  Real temperature,
529  Real ddensity_dT,
530  Real & mu,
531  Real & dmu_drho,
532  Real & dmu_dT) const
533 {
534  // Constants from Release on the IAPWS Formulation 2008 for the Viscosity of
535  // Ordinary Water Substance
536  const std::vector<int> I{0, 1, 2, 3, 0, 1, 2, 3, 5, 0, 1, 2, 3, 4, 0, 1, 0, 3, 4, 3, 5};
537  const std::vector<int> J{0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 4, 4, 5, 6, 6};
538  const std::vector<Real> H0{1.67752, 2.20462, 0.6366564, -0.241605};
539  const std::vector<Real> H1{
540  5.20094e-1, 8.50895e-2, -1.08374, -2.89555e-1, 2.22531e-1, 9.99115e-1, 1.88797,
541  1.26613, 1.20573e-1, -2.81378e-1, -9.06851e-1, -7.72479e-1, -4.89837e-1, -2.57040e-1,
542  1.61913e-1, 2.57399e-1, -3.25372e-2, 6.98452e-2, 8.72102e-3, -4.35673e-3, -5.93264e-4};
543 
544  Real mu_star = 1.0e-6;
545  Real rhobar = density / _rho_critical;
546  Real Tbar = temperature / _T_critical;
547  Real drhobar_drho = 1.0 / _rho_critical;
548  Real dTbar_dT = 1.0 / _T_critical;
549 
550  // Limit of zero density. Derivative wrt rho is 0
551  Real sum0 = 0.0, dsum0_dTbar = 0.0;
552  for (std::size_t i = 0; i < H0.size(); ++i)
553  {
554  sum0 += H0[i] / std::pow(Tbar, i);
555  dsum0_dTbar -= i * H0[i] / std::pow(Tbar, i + 1);
556  }
557 
558  Real mu0 = 100.0 * std::sqrt(Tbar) / sum0;
559  Real dmu0_dTbar =
560  50.0 / std::sqrt(Tbar) / sum0 - 100.0 * std::sqrt(Tbar) * dsum0_dTbar / sum0 / sum0;
561 
562  // Residual component due to finite density
563  Real sum1 = 0.0, dsum1_drho = 0.0, dsum1_dTbar = 0.0;
564  for (std::size_t i = 0; i < H1.size(); ++i)
565  {
566  sum1 += std::pow(1.0 / Tbar - 1.0, I[i]) * H1[i] * std::pow(rhobar - 1.0, J[i]);
567  dsum1_drho +=
568  std::pow(1.0 / Tbar - 1.0, I[i]) * H1[i] * J[i] * std::pow(rhobar - 1.0, J[i] - 1);
569  dsum1_dTbar -= I[i] * std::pow(1.0 / Tbar - 1.0, I[i] - 1) * H1[i] *
570  std::pow(rhobar - 1.0, J[i]) / Tbar / Tbar;
571  }
572 
573  Real mu1 = std::exp(rhobar * sum1);
574  Real dmu1_drho = (sum1 + rhobar * dsum1_drho) * mu1;
575  Real dmu1_dTbar = (rhobar * dsum1_dTbar) * mu1;
576 
577  // Viscosity and its derivatives are then
578  mu = mu_star * mu0 * mu1;
579  dmu_drho = mu_star * mu0 * dmu1_drho * drhobar_drho;
580  dmu_dT = mu_star * (dmu0_dTbar * mu1 + mu0 * dmu1_dTbar) * dTbar_dT + dmu_drho * ddensity_dT;
581 }
582 
583 Real
585 {
586  Real rho = this->rho(pressure, temperature);
587  return this->k_from_rho_T(rho, temperature);
588 }
589 
590 void
592  Real /*pressure*/, Real /*temperature*/, Real & /*k*/, Real & /*dk_dp*/, Real & /*dk_dT*/) const
593 {
594  mooseError(name(), "k_dpT() is not implemented");
595 }
596 
597 Real
599 {
600  // Scale the density and temperature. Note that the scales are slightly
601  // different to the critical values used in IAPWS-IF97
602  Real Tbar = temperature / 647.26;
603  Real rhobar = density / 317.7;
604  std::vector<Real> a{0.0102811, 0.0299621, 0.0156146, -0.00422464};
605 
606  // Ideal gas component
607  Real sum0 = 0.0;
608 
609  for (std::size_t i = 0; i < a.size(); ++i)
610  sum0 += a[i] * std::pow(Tbar, i);
611 
612  Real lambda0 = std::sqrt(Tbar) * sum0;
613 
614  // The contribution due to finite density
615  Real lambda1 =
616  -0.39707 + 0.400302 * rhobar + 1.06 * std::exp(-0.171587 * std::pow(rhobar + 2.392190, 2.0));
617 
618  // Critical enhancement
619  Real DeltaT = std::abs(Tbar - 1.0) + 0.00308976;
620  Real Q = 2.0 + 0.0822994 / std::pow(DeltaT, 0.6);
621  Real S = (Tbar >= 1.0 ? 1.0 / DeltaT : 10.0932 / std::pow(DeltaT, 0.6));
622 
623  Real lambda2 = (0.0701309 / std::pow(Tbar, 10.0) + 0.011852) * std::pow(rhobar, 1.8) *
624  std::exp(0.642857 * (1.0 - std::pow(rhobar, 2.8))) +
625  0.00169937 * S * std::pow(rhobar, Q) *
626  std::exp((Q / (1.0 + Q)) * (1.0 - std::pow(rhobar, 1.0 + Q))) -
627  1.02 * std::exp(-4.11717 * std::pow(Tbar, 1.5) - 6.17937 / std::pow(rhobar, 5.0));
628 
629  return lambda0 + lambda1 + lambda2;
630 }
631 
632 Real
634 {
635  Real entropy, pi, tau, density3, delta;
636 
637  // Determine which region the point is in
638  unsigned int region = inRegion(pressure, temperature);
639  switch (region)
640  {
641  case 1:
642  pi = pressure / _p_star[0];
643  tau = _T_star[0] / temperature;
644  entropy = _Rw * (tau * dgamma1_dtau(pi, tau) - gamma1(pi, tau));
645  break;
646 
647  case 2:
648  pi = pressure / _p_star[1];
649  tau = _T_star[1] / temperature;
650  entropy = _Rw * (tau * dgamma2_dtau(pi, tau) - gamma2(pi, tau));
651  break;
652 
653  case 3:
654  // Calculate density first, then use that in Helmholtz free energy
655  density3 = densityRegion3(pressure, temperature);
656  delta = density3 / _rho_critical;
657  tau = _T_star[2] / temperature;
658  entropy = _Rw * (tau * dphi3_dtau(delta, tau) - phi3(delta, tau));
659  break;
660 
661  case 5:
662  pi = pressure / _p_star[4];
663  tau = _T_star[4] / temperature;
664  entropy = _Rw * (tau * dgamma5_dtau(pi, tau) - gamma5(pi, tau));
665  break;
666 
667  default:
668  mooseError(name(), ": inRegion() has given an incorrect region");
669  }
670  return entropy;
671 }
672 
673 Real
675 {
676  Real enthalpy, pi, tau, delta;
677 
678  // Determine which region the point is in
679  unsigned int region = inRegion(pressure, temperature);
680  switch (region)
681  {
682  case 1:
683  pi = pressure / _p_star[0];
684  tau = _T_star[0] / temperature;
685  enthalpy = _Rw * _T_star[0] * dgamma1_dtau(pi, tau);
686  break;
687 
688  case 2:
689  pi = pressure / _p_star[1];
690  tau = _T_star[1] / temperature;
691  enthalpy = _Rw * _T_star[1] * dgamma2_dtau(pi, tau);
692  break;
693 
694  case 3:
695  {
696  // Calculate density first, then use that in Helmholtz free energy
697  Real density3 = densityRegion3(pressure, temperature);
698  delta = density3 / _rho_critical;
699  tau = _T_star[2] / temperature;
700  enthalpy =
701  _Rw * temperature * (tau * dphi3_dtau(delta, tau) + delta * dphi3_ddelta(delta, tau));
702  break;
703  }
704 
705  case 5:
706  pi = pressure / _p_star[4];
707  tau = _T_star[4] / temperature;
708  enthalpy = _Rw * _T_star[4] * dgamma5_dtau(pi, tau);
709  break;
710 
711  default:
712  mooseError("Water97FluidProperties::inRegion has given an incorrect region");
713  }
714  return enthalpy;
715 }
716 
717 void
719  Real pressure, Real temperature, Real & h, Real & dh_dp, Real & dh_dT) const
720 {
721  Real enthalpy, pi, tau, delta, denthalpy_dp, denthalpy_dT;
722 
723  // Determine which region the point is in
724  unsigned int region = inRegion(pressure, temperature);
725  switch (region)
726  {
727  case 1:
728  pi = pressure / _p_star[0];
729  tau = _T_star[0] / temperature;
730  enthalpy = _Rw * _T_star[0] * dgamma1_dtau(pi, tau);
731  denthalpy_dp = _Rw * _T_star[0] * d2gamma1_dpitau(pi, tau) / _p_star[0];
732  denthalpy_dT = -_Rw * tau * tau * d2gamma1_dtau2(pi, tau);
733  break;
734 
735  case 2:
736  pi = pressure / _p_star[1];
737  tau = _T_star[1] / temperature;
738  enthalpy = _Rw * _T_star[1] * dgamma2_dtau(pi, tau);
739  denthalpy_dp = _Rw * _T_star[1] * d2gamma2_dpitau(pi, tau) / _p_star[1];
740  denthalpy_dT = -_Rw * tau * tau * d2gamma2_dtau2(pi, tau);
741  break;
742 
743  case 3:
744  {
745  // Calculate density first, then use that in Helmholtz free energy
746  Real density3 = densityRegion3(pressure, temperature);
747  delta = density3 / _rho_critical;
748  tau = _T_star[2] / temperature;
749  Real dpdd = dphi3_ddelta(delta, tau);
750  Real d2pddt = d2phi3_ddeltatau(delta, tau);
751  Real d2pdd2 = d2phi3_ddelta2(delta, tau);
752  enthalpy = _Rw * temperature * (tau * dphi3_dtau(delta, tau) + delta * dpdd);
753  denthalpy_dp = (d2pddt + dpdd + delta * d2pdd2) / _rho_critical /
754  (2.0 * delta * dpdd + delta * delta * d2pdd2);
755  denthalpy_dT = _Rw * delta * dpdd * (1.0 - tau * d2pddt / dpdd) *
756  (1.0 - tau * d2pddt / dpdd) / (2.0 + delta * d2pdd2 / dpdd) -
757  _Rw * tau * tau * d2phi3_dtau2(delta, tau);
758  break;
759  }
760 
761  case 5:
762  pi = pressure / _p_star[4];
763  tau = _T_star[4] / temperature;
764  enthalpy = _Rw * _T_star[4] * dgamma5_dtau(pi, tau);
765  denthalpy_dp = _Rw * _T_star[4] * d2gamma5_dpitau(pi, tau) / _p_star[4];
766  denthalpy_dT = -_Rw * tau * tau * d2gamma5_dtau2(pi, tau);
767  break;
768 
769  default:
770  mooseError("Water97FluidProperties::inRegion has given an incorrect region");
771  }
772  h = enthalpy;
773  dh_dp = denthalpy_dp;
774  dh_dT = denthalpy_dT;
775 }
776 
777 Real Water97FluidProperties::beta(Real /*pressure*/, Real /*temperature*/) const
778 {
779  mooseError(name(), ": beta() not implemented yet");
780 }
781 
782 Real
784 {
785  // Check whether the input temperature is within the region of validity of this equation.
786  // Valid for 273.15 K <= t <= 647.096 K
787  if (temperature < 273.15 || temperature > _T_critical)
788  mooseError(name(),
789  ": vaporPressure(): Temperature is outside range 273.15 K <= T "
790  "<= 647.096 K");
791 
792  // Constants for region 4 (the saturation curve up to the critical point)
793  const std::vector<Real> n4{0.11670521452767e4,
794  -0.72421316703206e6,
795  -0.17073846940092e2,
796  0.12020824702470e5,
797  -0.32325550322333e7,
798  0.14915108613530e2,
799  -0.48232657361591e4,
800  0.40511340542057e6,
801  -0.238555575678490,
802  0.65017534844798e3};
803 
804  Real theta, theta2, a, b, c, p;
805  theta = temperature + n4[8] / (temperature - n4[9]);
806  theta2 = theta * theta;
807  a = theta2 + n4[0] * theta + n4[1];
808  b = n4[2] * theta2 + n4[3] * theta + n4[4];
809  c = n4[5] * theta2 + n4[6] * theta + n4[7];
810  p = std::pow(2.0 * c / (-b + std::sqrt(b * b - 4.0 * a * c)), 4.0);
811 
812  return p * 1.e6;
813 }
814 
815 void
816 Water97FluidProperties::vaporPressure_dT(Real temperature, Real & psat, Real & dpsat_dT) const
817 {
818  // Check whether the input temperature is within the region of validity of this equation.
819  // Valid for 273.15 K <= t <= 647.096 K
820  if (temperature < 273.15 || temperature > _T_critical)
821  mooseError(name(),
822  ": vaporPressure_dT(): Temperature is outside range 273.15 K <= "
823  "T <= 647.096 K");
824 
825  // Constants for region 4 (the saturation curve up to the critical point)
826  const std::vector<Real> n4{0.11670521452767e4,
827  -0.72421316703206e6,
828  -0.17073846940092e2,
829  0.12020824702470e5,
830  -0.32325550322333e7,
831  0.14915108613530e2,
832  -0.48232657361591e4,
833  0.40511340542057e6,
834  -0.238555575678490,
835  0.65017534844798e3};
836 
837  Real theta, dtheta_dT, theta2, a, b, c, da_dtheta, db_dtheta, dc_dtheta;
838  theta = temperature + n4[8] / (temperature - n4[9]);
839  dtheta_dT = 1.0 - n4[8] / (temperature - n4[9]) / (temperature - n4[9]);
840  theta2 = theta * theta;
841 
842  a = theta2 + n4[0] * theta + n4[1];
843  b = n4[2] * theta2 + n4[3] * theta + n4[4];
844  c = n4[5] * theta2 + n4[6] * theta + n4[7];
845 
846  da_dtheta = 2.0 * theta + n4[0];
847  db_dtheta = 2.0 * n4[2] * theta + n4[3];
848  dc_dtheta = 2.0 * n4[5] * theta + n4[6];
849 
850  Real denominator = -b + std::sqrt(b * b - 4.0 * a * c);
851 
852  psat = std::pow(2.0 * c / denominator, 4.0) * 1.0e6;
853 
854  // The derivative wrt temperature is given by the chain rule
855  Real dpsat = 4.0 * std::pow(2.0 * c / denominator, 3.0);
856  dpsat *= (2.0 * dc_dtheta / denominator -
857  2.0 * c / denominator / denominator *
858  (-db_dtheta +
859  std::pow(b * b - 4.0 * a * c, -0.5) *
860  (b * db_dtheta - 2.0 * da_dtheta * c - 2.0 * a * dc_dtheta)));
861  dpsat_dT = dpsat * dtheta_dT * 1.0e6;
862 }
863 
864 Real
866 {
867  // Check whether the input pressure is within the region of validity of this equation.
868  // Valid for 611.213 Pa <= p <= 22.064 MPa
869  if (pressure < 611.23 || pressure > _p_critical)
870  mooseError(name(),
871  ": vaporTemperature(): Pressure is outside range 611.213 Pa <= "
872  "p <= 22.064 MPa");
873 
874  // Constants for region 4 (the saturation curve up to the critical point)
875  const std::vector<Real> n4{0.11670521452767e4,
876  -0.72421316703206e6,
877  -0.17073846940092e2,
878  0.12020824702470e5,
879  -0.32325550322333e7,
880  0.14915108613530e2,
881  -0.48232657361591e4,
882  0.40511340542057e6,
883  -0.238555575678490,
884  0.65017534844798e3};
885 
886  Real beta, beta2, e, f, g, d;
887  beta = std::pow(pressure / 1.e6, 0.25);
888  beta2 = beta * beta;
889  e = beta2 + n4[2] * beta + n4[5];
890  f = n4[0] * beta2 + n4[3] * beta + n4[6];
891  g = n4[1] * beta2 + n4[4] * beta + n4[7];
892  d = 2.0 * g / (-f - std::sqrt(f * f - 4.0 * e * g));
893 
894  return (n4[9] + d - std::sqrt((n4[9] + d) * (n4[9] + d) - 4.0 * (n4[8] + n4[9] * d))) / 2.0;
895 }
896 
897 Real
899 {
900  // Check whether the input temperature is within the region of validity of this equation.
901  // Valid for 623.15 K <= t <= 863.15 K
902  if (temperature < 623.15 || temperature > 863.15)
903  mooseError(name(), ": b23p(): Temperature is outside range of 623.15 K <= T <= 863.15 K");
904 
905  // Constants for the boundary between regions 2 and 3
906  const std::vector<Real> n23{0.34805185628969e3,
907  -0.11671859879975e1,
908  0.10192970039326e-2,
909  0.57254459862746e3,
910  0.13918839778870e2};
911 
912  return (n23[0] + n23[1] * temperature + n23[2] * temperature * temperature) * 1.e6;
913 }
914 
915 Real
917 {
918  // Check whether the input pressure is within the region of validity of this equation.
919  // Valid for 16.529 MPa <= p <= 100 MPa
920  if (pressure < 16.529e6 || pressure > 100.0e6)
921  mooseError(name(), ": b23T(): Pressure is outside range 16.529 MPa <= p <= 100 MPa");
922 
923  // Constants for the boundary between regions 2 and 3
924  const std::vector<Real> n23{0.34805185628969e3,
925  -0.11671859879975e1,
926  0.10192970039326e-2,
927  0.57254459862746e3,
928  0.13918839778870e2};
929 
930  return n23[3] + std::sqrt((pressure / 1.e6 - n23[4]) / n23[2]);
931 }
932 
933 unsigned int
935 {
936  // Valid for 273.15 K <= T <= 1073.15 K, p <= 100 MPa
937  // 1073.15 K <= T <= 2273.15 K, p <= 50 Mpa
938  if (temperature >= 273.15 && temperature <= 1073.15)
939  {
940  if (pressure < vaporPressure(273.15) || pressure > 100.0e6)
941  mooseError("Pressure ", pressure, " is out of range in ", name(), ": inRegion()");
942  }
943  else if (temperature > 1073.15 && temperature <= 2273.15)
944  {
945  if (pressure < 0.0 || pressure > 50.0e6)
946  mooseError("Pressure ", pressure, " is out of range in ", name(), ": inRegion()");
947  }
948  else
949  mooseError("Temperature ", temperature, " is out of range in ", name(), ": inRegion()");
950 
951  // Determine the phase region that the (P, T) point lies in
952  unsigned int region;
953 
954  if (temperature >= 273.15 && temperature <= 623.15)
955  {
956  if (pressure > vaporPressure(temperature) && pressure <= 100.0e6)
957  region = 1;
958  else
959  region = 2;
960  }
961  else if (temperature > 623.15 && temperature <= 863.15)
962  {
963  if (pressure <= b23p(temperature))
964  region = 2;
965  else
966  region = 3;
967  }
968  else if (temperature > 863.15 && temperature <= 1073.15)
969  region = 2;
970  else
971  region = 5;
972 
973  return region;
974 }
975 
976 Real
977 Water97FluidProperties::gamma1(Real pi, Real tau) const
978 {
979  Real sum = 0.0;
980  for (std::size_t i = 0; i < _n1.size(); ++i)
981  sum += _n1[i] * std::pow(7.1 - pi, _I1[i]) * std::pow(tau - 1.222, _J1[i]);
982 
983  return sum;
984 }
985 
986 Real
987 Water97FluidProperties::dgamma1_dpi(Real pi, Real tau) const
988 {
989  Real sum = 0.0;
990  for (std::size_t i = 0; i < _n1.size(); ++i)
991  sum += -_n1[i] * _I1[i] * std::pow(7.1 - pi, _I1[i] - 1) * std::pow(tau - 1.222, _J1[i]);
992 
993  return sum;
994 }
995 
996 Real
998 {
999  Real sum = 0.0;
1000  for (std::size_t i = 0; i < _n1.size(); ++i)
1001  sum += _n1[i] * _I1[i] * (_I1[i] - 1) * std::pow(7.1 - pi, _I1[i] - 2) *
1002  std::pow(tau - 1.222, _J1[i]);
1003 
1004  return sum;
1005 }
1006 
1007 Real
1009 {
1010  Real g = 0.0;
1011  for (std::size_t i = 0; i < _n1.size(); ++i)
1012  g += _n1[i] * _J1[i] * std::pow(7.1 - pi, _I1[i]) * std::pow(tau - 1.222, _J1[i] - 1);
1013 
1014  return g;
1015 }
1016 
1017 Real
1019 {
1020  Real dg = 0.0;
1021  for (std::size_t i = 0; i < _n1.size(); ++i)
1022  dg += _n1[i] * _J1[i] * (_J1[i] - 1) * std::pow(7.1 - pi, _I1[i]) *
1023  std::pow(tau - 1.222, _J1[i] - 2);
1024 
1025  return dg;
1026 }
1027 
1028 Real
1030 {
1031  Real dg = 0.0;
1032  for (std::size_t i = 0; i < _n1.size(); ++i)
1033  dg += -_n1[i] * _I1[i] * _J1[i] * std::pow(7.1 - pi, _I1[i] - 1) *
1034  std::pow(tau - 1.222, _J1[i] - 1);
1035 
1036  return dg;
1037 }
1038 
1039 Real
1040 Water97FluidProperties::gamma2(Real pi, Real tau) const
1041 {
1042  // Ideal gas part of the Gibbs free energy
1043  Real sum0 = 0.0;
1044  for (std::size_t i = 0; i < _n02.size(); ++i)
1045  sum0 += _n02[i] * std::pow(tau, _J02[i]);
1046 
1047  Real g0 = std::log(pi) + sum0;
1048 
1049  // Residual part of the Gibbs free energy
1050  Real gr = 0.0;
1051  for (std::size_t i = 0; i < _n2.size(); ++i)
1052  gr += _n2[i] * std::pow(pi, _I2[i]) * std::pow(tau - 0.5, _J2[i]);
1053 
1054  return g0 + gr;
1055 }
1056 
1057 Real
1058 Water97FluidProperties::dgamma2_dpi(Real pi, Real tau) const
1059 {
1060  // Ideal gas part of the Gibbs free energy
1061  Real dg0 = 1.0 / pi;
1062 
1063  // Residual part of the Gibbs free energy
1064  Real dgr = 0.0;
1065  for (std::size_t i = 0; i < _n2.size(); ++i)
1066  dgr += _n2[i] * _I2[i] * std::pow(pi, _I2[i] - 1) * std::pow(tau - 0.5, _J2[i]);
1067 
1068  return dg0 + dgr;
1069 }
1070 
1071 Real
1073 {
1074  // Ideal gas part of the Gibbs free energy
1075  Real dg0 = -1.0 / pi / pi;
1076 
1077  // Residual part of the Gibbs free energy
1078  Real dgr = 0.0;
1079  for (std::size_t i = 0; i < _n2.size(); ++i)
1080  dgr += _n2[i] * _I2[i] * (_I2[i] - 1) * std::pow(pi, _I2[i] - 2) * std::pow(tau - 0.5, _J2[i]);
1081 
1082  return dg0 + dgr;
1083 }
1084 
1085 Real
1087 {
1088  // Ideal gas part of the Gibbs free energy
1089  Real dg0 = 0.0;
1090  for (std::size_t i = 0; i < _n02.size(); ++i)
1091  dg0 += _n02[i] * _J02[i] * std::pow(tau, _J02[i] - 1);
1092 
1093  // Residual part of the Gibbs free energy
1094  Real dgr = 0.0;
1095  for (std::size_t i = 0; i < _n2.size(); ++i)
1096  dgr += _n2[i] * _J2[i] * std::pow(pi, _I2[i]) * std::pow(tau - 0.5, _J2[i] - 1);
1097 
1098  return dg0 + dgr;
1099 }
1100 
1101 Real
1103 {
1104  // Ideal gas part of the Gibbs free energy
1105  Real dg0 = 0.0;
1106  for (std::size_t i = 0; i < _n02.size(); ++i)
1107  dg0 += _n02[i] * _J02[i] * (_J02[i] - 1) * std::pow(tau, _J02[i] - 2);
1108 
1109  // Residual part of the Gibbs free energy
1110  Real dgr = 0.0;
1111  for (std::size_t i = 0; i < _n2.size(); ++i)
1112  dgr += _n2[i] * _J2[i] * (_J2[i] - 1) * std::pow(pi, _I2[i]) * std::pow(tau - 0.5, _J2[i] - 2);
1113 
1114  return dg0 + dgr;
1115 }
1116 
1117 Real
1119 {
1120  // Ideal gas part of the Gibbs free energy
1121  Real dg0 = 0.0;
1122 
1123  // Residual part of the Gibbs free energy
1124  Real dgr = 0.0;
1125  for (std::size_t i = 0; i < _n2.size(); ++i)
1126  dgr += _n2[i] * _I2[i] * _J2[i] * std::pow(pi, _I2[i] - 1) * std::pow(tau - 0.5, _J2[i] - 1);
1127 
1128  return dg0 + dgr;
1129 }
1130 
1131 Real
1132 Water97FluidProperties::phi3(Real delta, Real tau) const
1133 {
1134  Real sum = 0.0;
1135  for (std::size_t i = 1; i < _n3.size(); ++i)
1136  sum += _n3[i] * std::pow(delta, _I3[i]) * std::pow(tau, _J3[i]);
1137 
1138  return _n3[0] * std::log(delta) + sum;
1139 }
1140 
1141 Real
1142 Water97FluidProperties::dphi3_ddelta(Real delta, Real tau) const
1143 {
1144  Real sum = 0.0;
1145  for (std::size_t i = 1; i < _n3.size(); ++i)
1146  sum += _n3[i] * _I3[i] * std::pow(delta, _I3[i] - 1) * std::pow(tau, _J3[i]);
1147 
1148  return _n3[0] / delta + sum;
1149 }
1150 
1151 Real
1152 Water97FluidProperties::d2phi3_ddelta2(Real delta, Real tau) const
1153 {
1154  Real sum = 0.0;
1155  for (std::size_t i = 1; i < _n3.size(); ++i)
1156  sum += _n3[i] * _I3[i] * (_I3[i] - 1) * std::pow(delta, _I3[i] - 2) * std::pow(tau, _J3[i]);
1157 
1158  return -_n3[0] / delta / delta + sum;
1159 }
1160 
1161 Real
1162 Water97FluidProperties::dphi3_dtau(Real delta, Real tau) const
1163 {
1164  Real sum = 0.0;
1165  for (std::size_t i = 1; i < _n3.size(); ++i)
1166  sum += _n3[i] * _J3[i] * std::pow(delta, _I3[i]) * std::pow(tau, _J3[i] - 1);
1167 
1168  return sum;
1169 }
1170 
1171 Real
1172 Water97FluidProperties::d2phi3_dtau2(Real delta, Real tau) const
1173 {
1174  Real sum = 0.0;
1175  for (std::size_t i = 1; i < _n3.size(); ++i)
1176  sum += _n3[i] * _J3[i] * (_J3[i] - 1) * std::pow(delta, _I3[i]) * std::pow(tau, _J3[i] - 2);
1177 
1178  return sum;
1179 }
1180 
1181 Real
1182 Water97FluidProperties::d2phi3_ddeltatau(Real delta, Real tau) const
1183 {
1184  Real sum = 0.0;
1185  for (std::size_t i = 1; i < _n3.size(); ++i)
1186  sum += _n3[i] * _I3[i] * _J3[i] * std::pow(delta, _I3[i] - 1) * std::pow(tau, _J3[i] - 1);
1187 
1188  return sum;
1189 }
1190 
1191 Real
1192 Water97FluidProperties::gamma5(Real pi, Real tau) const
1193 {
1194  // Ideal gas part of the Gibbs free energy
1195  Real sum0 = 0.0;
1196  for (std::size_t i = 0; i < _n05.size(); ++i)
1197  sum0 += _n05[i] * std::pow(tau, _J05[i]);
1198 
1199  Real g0 = std::log(pi) + sum0;
1200 
1201  // Residual part of the Gibbs free energy
1202  Real gr = 0.0;
1203  for (std::size_t i = 0; i < _n5.size(); ++i)
1204  gr += _n5[i] * std::pow(pi, _I5[i]) * std::pow(tau, _J5[i]);
1205 
1206  return g0 + gr;
1207 }
1208 
1209 Real
1210 Water97FluidProperties::dgamma5_dpi(Real pi, Real tau) const
1211 {
1212  // Ideal gas part of the Gibbs free energy
1213  Real dg0 = 1.0 / pi;
1214 
1215  // Residual part of the Gibbs free energy
1216  Real dgr = 0.0;
1217  for (std::size_t i = 0; i < _n5.size(); ++i)
1218  dgr += _n5[i] * _I5[i] * std::pow(pi, _I5[i] - 1) * std::pow(tau, _J5[i]);
1219 
1220  return dg0 + dgr;
1221 }
1222 
1223 Real
1225 {
1226  // Ideal gas part of the Gibbs free energy
1227  Real dg0 = -1.0 / pi / pi;
1228 
1229  // Residual part of the Gibbs free energy
1230  Real dgr = 0.0;
1231  for (std::size_t i = 0; i < _n5.size(); ++i)
1232  dgr += _n5[i] * _I5[i] * (_I5[i] - 1) * std::pow(pi, _I5[i] - 2) * std::pow(tau, _J5[i]);
1233 
1234  return dg0 + dgr;
1235 }
1236 
1237 Real
1239 {
1240  // Ideal gas part of the Gibbs free energy
1241  Real dg0 = 0.0;
1242  for (std::size_t i = 0; i < _n05.size(); ++i)
1243  dg0 += _n05[i] * _J05[i] * std::pow(tau, _J05[i] - 1);
1244 
1245  // Residual part of the Gibbs free energy
1246  Real dgr = 0.0;
1247  for (std::size_t i = 0; i < _n5.size(); ++i)
1248  dgr += _n5[i] * _J5[i] * std::pow(pi, _I5[i]) * std::pow(tau, _J5[i] - 1);
1249 
1250  return dg0 + dgr;
1251 }
1252 
1253 Real
1255 {
1256  // Ideal gas part of the Gibbs free energy
1257  Real dg0 = 0.0;
1258  for (std::size_t i = 0; i < _n05.size(); ++i)
1259  dg0 += _n05[i] * _J05[i] * (_J05[i] - 1) * std::pow(tau, _J05[i] - 2);
1260 
1261  // Residual part of the Gibbs free energy
1262  Real dgr = 0.0;
1263  for (std::size_t i = 0; i < _n5.size(); ++i)
1264  dgr += _n5[i] * _J5[i] * (_J5[i] - 1) * std::pow(pi, _I5[i]) * std::pow(tau, _J5[i] - 2);
1265 
1266  return dg0 + dgr;
1267 }
1268 
1269 Real
1271 {
1272  // Ideal gas part of the Gibbs free energy
1273  Real dg0 = 0.0;
1274 
1275  // Residual part of the Gibbs free energy
1276  Real dgr = 0.0;
1277  for (std::size_t i = 0; i < _n5.size(); ++i)
1278  dgr += _n5[i] * _I5[i] * _J5[i] * std::pow(pi, _I5[i] - 1) * std::pow(tau, _J5[i] - 1);
1279 
1280  return dg0 + dgr;
1281 }
1282 
1283 unsigned int
1285 {
1286  Real pMPa = pressure / 1.0e6;
1287  const Real P3cd = 19.00881189173929;
1288  unsigned int subregion = 0;
1289 
1290  if (pMPa > 40.0 && pMPa <= 100.0)
1291  {
1292  if (temperature <= tempXY(pressure, AB))
1293  subregion = 0;
1294  else // (temperature > tempXY(pressure, AB))
1295  subregion = 1;
1296  }
1297  else if (pMPa > 25.0 && pMPa <= 40.0)
1298  {
1299  if (temperature <= tempXY(pressure, CD))
1300  subregion = 2;
1301  else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, AB))
1302  subregion = 3;
1303  else if (temperature > tempXY(pressure, AB) && temperature <= tempXY(pressure, EF))
1304  subregion = 4;
1305  else // (temperature > tempXY(pressure, EF))
1306  subregion = 5;
1307  }
1308  else if (pMPa > 23.5 && pMPa <= 25.0)
1309  {
1310  if (temperature <= tempXY(pressure, CD))
1311  subregion = 2;
1312  else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, GH))
1313  subregion = 6;
1314  else if (temperature > tempXY(pressure, GH) && temperature <= tempXY(pressure, EF))
1315  subregion = 7;
1316  else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, IJ))
1317  subregion = 8;
1318  else if (temperature > tempXY(pressure, IJ) && temperature <= tempXY(pressure, JK))
1319  subregion = 9;
1320  else // (temperature > tempXY(pressure, JK))
1321  subregion = 10;
1322  }
1323  else if (pMPa > 23.0 && pMPa <= 23.5)
1324  {
1325  if (temperature <= tempXY(pressure, CD))
1326  subregion = 2;
1327  else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, GH))
1328  subregion = 11;
1329  else if (temperature > tempXY(pressure, GH) && temperature <= tempXY(pressure, EF))
1330  subregion = 7;
1331  else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, IJ))
1332  subregion = 8;
1333  else if (temperature > tempXY(pressure, IJ) && temperature <= tempXY(pressure, JK))
1334  subregion = 9;
1335  else // (temperature > tempXY(pressure, JK))
1336  subregion = 10;
1337  }
1338  else if (pMPa > 22.5 && pMPa <= 23.0)
1339  {
1340  if (temperature <= tempXY(pressure, CD))
1341  subregion = 2;
1342  else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, GH))
1343  subregion = 11;
1344  else if (temperature > tempXY(pressure, GH) && temperature <= tempXY(pressure, MN))
1345  subregion = 12;
1346  else if (temperature > tempXY(pressure, MN) && temperature <= tempXY(pressure, EF))
1347  subregion = 13;
1348  else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, OP))
1349  subregion = 14;
1350  else if (temperature > tempXY(pressure, OP) && temperature <= tempXY(pressure, IJ))
1351  subregion = 15;
1352  else if (temperature > tempXY(pressure, IJ) && temperature <= tempXY(pressure, JK))
1353  subregion = 9;
1354  else // (temperature > tempXY(pressure, JK))
1355  subregion = 10;
1356  }
1357  else if (pMPa > vaporPressure(643.15) * 1.0e-6 &&
1358  pMPa <= 22.5) // vaporPressure(643.15) = 21.04 MPa
1359  {
1360  if (temperature <= tempXY(pressure, CD))
1361  subregion = 2;
1362  else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, QU))
1363  subregion = 16;
1364  else if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, RX))
1365  {
1366  if (pMPa > 22.11 && pMPa <= 22.5)
1367  {
1368  if (temperature <= tempXY(pressure, UV))
1369  subregion = 20;
1370  else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1371  subregion = 21;
1372  else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1373  subregion = 22;
1374  else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1375  subregion = 23;
1376  }
1377  else if (pMPa > 22.064 && pMPa <= 22.11)
1378  {
1379  if (temperature <= tempXY(pressure, UV))
1380  subregion = 20;
1381  else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1382  subregion = 24;
1383  else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1384  subregion = 25;
1385  else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1386  subregion = 23;
1387  }
1388  else if (temperature <= vaporTemperature(pressure))
1389  {
1390  if (pMPa > 21.93161551 && pMPa <= 22.064)
1391  if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, UV))
1392  subregion = 20;
1393  else
1394  subregion = 24;
1395  else // (pMPa > vaporPressure(643.15) * 1.0e-6 && pMPa <= 21.93161551)
1396  subregion = 20;
1397  }
1398  else if (temperature > vaporTemperature(pressure))
1399  {
1400  if (pMPa > 21.90096265 && pMPa <= 22.064)
1401  {
1402  if (temperature <= tempXY(pressure, WX))
1403  subregion = 25;
1404  else
1405  subregion = 23;
1406  }
1407  else
1408  subregion = 23;
1409  }
1410  }
1411  else if (temperature > tempXY(pressure, RX) && temperature <= tempXY(pressure, JK))
1412  subregion = 17;
1413  else
1414  subregion = 10;
1415  }
1416  else if (pMPa > 20.5 &&
1417  pMPa <= vaporPressure(643.15) * 1.0e-6) // vaporPressure(643.15) = 21.04 MPa
1418  {
1419  if (temperature <= tempXY(pressure, CD))
1420  subregion = 2;
1421  else if (temperature > tempXY(pressure, CD) && temperature <= vaporTemperature(pressure))
1422  subregion = 16;
1423  else if (temperature > vaporTemperature(pressure) && temperature <= tempXY(pressure, JK))
1424  subregion = 17;
1425  else // (temperature > tempXY(pressure, JK))
1426  subregion = 10;
1427  }
1428  else if (pMPa > P3cd && pMPa <= 20.5) // P3cd = 19.00881189173929
1429  {
1430  if (temperature <= tempXY(pressure, CD))
1431  subregion = 2;
1432  else if (temperature > tempXY(pressure, CD) && temperature <= vaporTemperature(pressure))
1433  subregion = 18;
1434  else
1435  subregion = 19;
1436  }
1437  else if (pMPa > vaporPressure(623.15) * 1.0e-6 && pMPa <= P3cd)
1438  {
1439  if (temperature < vaporTemperature(pressure))
1440  subregion = 2;
1441  else
1442  subregion = 19;
1443  }
1444  else if (pMPa > 22.11 && pMPa <= 22.5)
1445  {
1446  if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, UV))
1447  subregion = 20;
1448  else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1449  subregion = 21;
1450  else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1451  subregion = 22;
1452  else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1453  subregion = 23;
1454  }
1455  else if (pMPa > 22.064 && pMPa <= 22.11)
1456  {
1457  if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, UV))
1458  subregion = 20;
1459  else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1460  subregion = 24;
1461  else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1462  subregion = 25;
1463  else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1464  subregion = 23;
1465  }
1466  else
1467  mooseError(name(), ": subregion3(): Shouldn't have got here!");
1468 
1469  return subregion;
1470 }
1471 
1472 Real
1474 {
1475  Real pi = pressure / 1.0e6;
1476 
1477  const std::vector<std::vector<int>> I{{0, 1, 2, -1, -2},
1478  {0, 1, 2, 3},
1479  {0, 1, 2, 3, 4},
1480  {0, 1, 2, 3, 4},
1481  {0, 1, 2, 3, 4},
1482  {0, 1, 2, 3},
1483  {0, 1, 2, -1, -2},
1484  {0, 1, 2, 3},
1485  {0, 1, 2, 3},
1486  {0, 1, 2, 3, 4},
1487  {0, 1, 2, -1, -2}};
1488 
1489  const std::vector<std::vector<Real>> n{
1490  {0.154793642129415e4,
1491  -0.187661219490113e3,
1492  0.213144632222113e2,
1493  -0.191887498864292e4,
1494  0.918419702359447e3},
1495  {0.585276966696349e3, 0.278233532206915e1, -0.127283549295878e-1, 0.159090746562729e-3},
1496  {-0.249284240900418e5,
1497  0.428143584791546e4,
1498  -0.269029173140130e3,
1499  0.751608051114157e1,
1500  -0.787105249910383e-1},
1501  {0.584814781649163e3,
1502  -0.616179320924617,
1503  0.260763050899562,
1504  -0.587071076864459e-2,
1505  0.515308185433082e-4},
1506  {0.617229772068439e3,
1507  -0.770600270141675e1,
1508  0.697072596851896,
1509  -0.157391839848015e-1,
1510  0.137897492684194e-3},
1511  {0.535339483742384e3, 0.761978122720128e1, -0.158365725441648, 0.192871054508108e-2},
1512  {0.969461372400213e3,
1513  -0.332500170441278e3,
1514  0.642859598466067e2,
1515  0.773845935768222e3,
1516  -0.152313732937084e4},
1517  {0.565603648239126e3, 0.529062258221222e1, -0.102020639611016, 0.122240301070145e-2},
1518  {0.584561202520006e3, -0.102961025163669e1, 0.243293362700452, -0.294905044740799e-2},
1519  {0.528199646263062e3, 0.890579602135307e1, -0.222814134903755, 0.286791682263697e-2},
1520  {0.728052609145380e1,
1521  0.973505869861952e2,
1522  0.147370491183191e2,
1523  0.329196213998375e3,
1524  0.873371668682417e3}};
1525 
1526  // Choose the constants based on the string xy
1527  unsigned int row;
1528 
1529  switch (xy)
1530  {
1531  case AB:
1532  row = 0;
1533  break;
1534  case CD:
1535  row = 1;
1536  break;
1537  case GH:
1538  row = 2;
1539  break;
1540  case IJ:
1541  row = 3;
1542  break;
1543  case JK:
1544  row = 4;
1545  break;
1546  case MN:
1547  row = 5;
1548  break;
1549  case OP:
1550  row = 6;
1551  break;
1552  case QU:
1553  row = 7;
1554  break;
1555  case RX:
1556  row = 8;
1557  break;
1558  case UV:
1559  row = 9;
1560  break;
1561  case WX:
1562  row = 10;
1563  break;
1564  default:
1565  row = 0;
1566  }
1567 
1568  Real sum = 0.0;
1569 
1570  if (xy == AB || xy == OP || xy == WX)
1571  for (std::size_t i = 0; i < n[row].size(); ++i)
1572  sum += n[row][i] * std::pow(std::log(pi), I[row][i]);
1573  else if (xy == EF)
1574  sum += 3.727888004 * (pi - _p_critical / 1.0e6) + _T_critical;
1575  else
1576  for (std::size_t i = 0; i < n[row].size(); ++i)
1577  sum += n[row][i] * std::pow(pi, I[row][i]);
1578 
1579  return sum;
1580 }
1581 
1582 Real
1584  Real pi, Real theta, Real a, Real b, Real c, Real d, Real e, unsigned int sid) const
1585 {
1586  Real sum = 0.0;
1587 
1588  for (std::size_t i = 0; i < _n3s[sid].size(); ++i)
1589  sum += _n3s[sid][i] * std::pow(std::pow(pi - a, c), _I3s[sid][i]) *
1590  std::pow(std::pow(theta - b, d), _J3s[sid][i]);
1591 
1592  return std::pow(sum, e);
1593 }
1594 
1595 Real
1597 {
1598  // Region 3 is subdivided into 26 subregions, each with a given backwards equation
1599  // to directly calculate density from pressure and temperature without the need for
1600  // expensive iterations. Find the subregion id that the point is in:
1601  unsigned int sid = subregion3(pressure, temperature);
1602 
1603  Real vstar, pi, theta, a, b, c, d, e;
1604  unsigned int N;
1605 
1606  vstar = _par3[sid][0];
1607  pi = pressure / _par3[sid][1] / 1.0e6;
1608  theta = temperature / _par3[sid][2];
1609  a = _par3[sid][3];
1610  b = _par3[sid][4];
1611  c = _par3[sid][5];
1612  d = _par3[sid][6];
1613  e = _par3[sid][7];
1614  N = _par3N[sid];
1615 
1616  Real sum = 0.0;
1617  Real volume = 0.0;
1618 
1619  // Note that subregion 13 is the only different formulation
1620  if (sid == 13)
1621  {
1622  for (std::size_t i = 0; i < N; ++i)
1623  sum += _n3s[sid][i] * std::pow(pi - a, _I3s[sid][i]) * std::pow(theta - b, _J3s[sid][i]);
1624 
1625  volume = vstar * std::exp(sum);
1626  }
1627  else
1628  volume = vstar * subregionVolume(pi, theta, a, b, c, d, e, sid);
1629 
1630  // Density is the inverse of volume
1631  return 1.0 / volume;
1632 }
1633 
1634 Real Water97FluidProperties::henryConstant(Real /*temperature*/) const
1635 {
1636  mooseError(name(), ": henryConstant() not defined");
1637 }
1638 
1639 void
1641  Real & /* Kh */,
1642  Real & /* dKh_dT */) const
1643 {
1644  mooseError(name(), ": henryConstant_dT() not defined");
1645 }
1646 
1647 unsigned int
1649 {
1650  unsigned int region;
1651 
1652  // Need to calculate enthalpies at the boundaries to delineate regions
1653  Real p273 = vaporPressure(273.15);
1654  Real p623 = vaporPressure(623.15);
1655 
1656  if (pressure >= p273 && pressure <= p623)
1657  {
1658  if (enthalpy >= h(pressure, 273.15) && enthalpy <= h(pressure, vaporTemperature(pressure)))
1659  region = 1;
1660  else if (enthalpy > h(pressure, vaporTemperature(pressure)) && enthalpy <= h(pressure, 1073.15))
1661  region = 2;
1662  else if (enthalpy > h(pressure, 1073.15) && enthalpy <= h(pressure, 2273.15))
1663  region = 5;
1664  else
1665  mooseError("Enthalpy ", enthalpy, " is out of range in ", name(), ": inRegionPH()");
1666  }
1667  else if (pressure > p623 && pressure <= 50.0e6)
1668  {
1669  if (enthalpy >= h(pressure, 273.15) && enthalpy <= h(pressure, 623.15))
1670  region = 1;
1671  else if (enthalpy > h(pressure, 623.15) && enthalpy <= h(pressure, b23T(pressure)))
1672  region = 3;
1673  else if (enthalpy > h(pressure, b23T(pressure)) && enthalpy <= h(pressure, 1073.15))
1674  region = 2;
1675  else if (enthalpy > h(pressure, 1073.15) && enthalpy <= h(pressure, 2273.15))
1676  region = 5;
1677  else
1678  mooseError("Enthalpy ", enthalpy, " is out of range in ", name(), ": inRegionPH()");
1679  }
1680  else if (pressure > 50.0e6 && pressure <= 100.0e6)
1681  {
1682  if (enthalpy >= h(pressure, 273.15) && enthalpy <= h(pressure, 623.15))
1683  region = 1;
1684  else if (enthalpy > h(pressure, 623.15) && enthalpy <= h(pressure, b23T(pressure)))
1685  region = 3;
1686  else if (enthalpy > h(pressure, b23T(pressure)) && enthalpy <= h(pressure, 1073.15))
1687  region = 2;
1688  else
1689  mooseError("Enthalpy ", enthalpy, " is out of range in ", name(), ": inRegionPH()");
1690  }
1691  else
1692  mooseError("Pressure ", pressure, " is out of range in ", name(), ": inRegionPH()");
1693 
1694  return region;
1695 }
1696 
1697 unsigned int
1699 {
1700  unsigned int subregion;
1701 
1702  if (pressure <= 4.0e6)
1703  subregion = 1;
1704  else if (pressure > 4.0e6 && pressure < 6.5467e6)
1705  subregion = 2;
1706  else
1707  {
1708  if (enthalpy >= b2bc(pressure))
1709  subregion = 2;
1710  else
1711  subregion = 3;
1712  }
1713 
1714  return subregion;
1715 }
1716 
1717 unsigned int
1719 {
1720  unsigned int subregion;
1721 
1722  if (enthalpy <= b3ab(pressure))
1723  subregion = 1;
1724  else
1725  subregion = 2;
1726 
1727  return subregion;
1728 }
1729 
1730 Real
1732 {
1733  // Check whether the input pressure is within the region of validity of this equation.
1734  if (pressure < 6.5467e6 || pressure > 100.0e6)
1735  mooseError(name(), ": b2bc(): Pressure is outside range of 6.5467 MPa <= p <= 100 MPa");
1736 
1737  Real pi = pressure / 1.0e6;
1738 
1739  return (0.26526571908428e4 + std::sqrt((pi - 0.45257578905948e1) / 0.12809002730136e-3)) * 1.0e3;
1740 }
1741 
1742 Real
1744 {
1745  Real temperature = 0.0;
1746 
1747  // Determine which region the point is in
1748  unsigned int region = inRegionPH(pressure, enthalpy);
1749 
1750  switch (region)
1751  {
1752  case 1:
1753  temperature = temperature_from_ph1(pressure, enthalpy);
1754  break;
1755 
1756  case 2:
1757  {
1758  // First, determine which subregion the point is in:
1759  unsigned int subregion = subregion2ph(pressure, enthalpy);
1760 
1761  if (subregion == 1)
1762  temperature = temperature_from_ph2a(pressure, enthalpy);
1763  else if (subregion == 2)
1764  temperature = temperature_from_ph2b(pressure, enthalpy);
1765  else
1766  temperature = temperature_from_ph2c(pressure, enthalpy);
1767  break;
1768  }
1769 
1770  case 3:
1771  {
1772  // First, determine which subregion the point is in:
1773  unsigned int subregion = subregion3ph(pressure, enthalpy);
1774 
1775  if (subregion == 1)
1776  temperature = temperature_from_ph3a(pressure, enthalpy);
1777  else
1778  temperature = temperature_from_ph3b(pressure, enthalpy);
1779  break;
1780  }
1781 
1782  case 5:
1783  mooseError(name(), ": temperature_from_ph() not implemented for region 5");
1784  break;
1785 
1786  default:
1787  mooseError(name(), ": inRegionPH() has given an incorrect region");
1788  }
1789 
1790  return temperature;
1791 }
1792 
1793 Real
1795 {
1796  Real pi = pressure / 1.0e6;
1797  Real eta = enthalpy / 2500.0e3;
1798  Real sum = 0.0;
1799 
1800  for (std::size_t i = 0; i < _nph1.size(); ++i)
1801  sum += _nph1[i] * std::pow(pi, _Iph1[i]) * std::pow(eta + 1.0, _Jph1[i]);
1802 
1803  return sum;
1804 }
1805 
1806 Real
1808 {
1809  Real pi = pressure / 1.0e6;
1810  Real eta = enthalpy / 2000.0e3;
1811  Real sum = 0.0;
1812 
1813  for (std::size_t i = 0; i < _nph2a.size(); ++i)
1814  sum += _nph2a[i] * std::pow(pi, _Iph2a[i]) * std::pow(eta - 2.1, _Jph2a[i]);
1815 
1816  return sum;
1817 }
1818 
1819 Real
1821 {
1822  Real pi = pressure / 1.0e6;
1823  Real eta = enthalpy / 2000.0e3;
1824  Real sum = 0.0;
1825 
1826  for (std::size_t i = 0; i < _nph2b.size(); ++i)
1827  sum += _nph2b[i] * std::pow(pi - 2.0, _Iph2b[i]) * std::pow(eta - 2.6, _Jph2b[i]);
1828 
1829  return sum;
1830 }
1831 
1832 Real
1834 {
1835  Real pi = pressure / 1.0e6;
1836  Real eta = enthalpy / 2000.0e3;
1837  Real sum = 0.0;
1838 
1839  for (std::size_t i = 0; i < _nph2c.size(); ++i)
1840  sum += _nph2c[i] * std::pow(pi + 25.0, _Iph2c[i]) * std::pow(eta - 1.8, _Jph2c[i]);
1841 
1842  return sum;
1843 }
1844 
1845 Real
1847 {
1848  // Check whether the input pressure is within the region of validity of this equation.
1849  if (pressure < b23p(623.15) || pressure > 100.0e6)
1850  mooseError(name(), ": b3ab(): Pressure is outside range of 16.529 MPa <= p <= 100 MPa");
1851 
1852  Real pi = pressure / 1.0e6;
1853  Real eta = 0.201464004206875e4 + 0.374696550136983e1 * pi - 0.219921901054187e-1 * pi * pi +
1854  0.87513168600995e-4 * pi * pi * pi;
1855 
1856  return eta * 1.0e3;
1857 }
1858 
1859 Real
1861 {
1862  Real pi = pressure / 100.0e6;
1863  Real eta = enthalpy / 2300.0e3;
1864  Real sum = 0.0;
1865 
1866  for (std::size_t i = 0; i < _nph3a.size(); ++i)
1867  sum += _nph3a[i] * std::pow(pi + 0.24, _Iph3a[i]) * std::pow(eta - 0.615, _Jph3a[i]);
1868 
1869  return sum * 760.0;
1870 }
1871 
1872 Real
1874 {
1875  Real pi = pressure / 100.0e6;
1876  Real eta = enthalpy / 2800.0e3;
1877  Real sum = 0.0;
1878 
1879  for (std::size_t i = 0; i < _nph3b.size(); ++i)
1880  sum += _nph3b[i] * std::pow(pi + 0.298, _Iph3b[i]) * std::pow(eta - 0.72, _Jph3b[i]);
1881 
1882  return sum * 860.0;
1883 }
const std::vector< int > _Jph2c
unsigned int inRegionPH(Real pressure, Real enthalpy) const
Determines the phase region that the given pressure and enthaply values lie in.
const Real _T_critical
Critical temperature (K)
const std::vector< int > _I5
const std::vector< Real > _n05
const std::vector< std::vector< Real > > _par3
const std::vector< int > _Iph3a
virtual Real c(Real pressure, Real temperature) const override
Speed of sound.
unsigned int subregion3(Real pressure, Real temperature) const
Provides the correct subregion index for a (P,T) point in region 3.
virtual Real beta(Real pressure, Real temperature) const override
Thermal expansion coefficient.
const std::vector< int > _I2
virtual Real criticalPressure() const
Water critical pressure.
const Real _p_critical
Critical pressure (Pa)
Real dphi3_ddelta(Real delta, Real tau) const
Derivative of Helmholtz free energy in Region 3 wrt delta.
Real d2gamma5_dtau2(Real pi, Real tau) const
Second derivative of Gibbs free energy in Region 5 wrt tau.
virtual Real cp(Real pressure, Real temperature) const override
Isobaric specific heat capacity.
virtual std::string fluidName() const override
Fluid name.
Real temperature_from_ph1(Real pressure, Real enthalpy) const
Backwards equation T(p, h) in Region 1 Eq.
const std::vector< int > _Jph2b
const std::vector< int > _Iph2b
const std::vector< int > _I1
virtual void henryConstant_dT(Real temperature, Real &Kh, Real &dKh_dT) const override
Henry&#39;s law constant for dissolution in water and derivative wrt temperature.
virtual Real cv(Real pressure, Real temperature) const override
Isochoric specific heat.
void vaporPressure_dT(Real temperature, Real &psat, Real &dpsat_dT) const
Saturation pressure as a function of temperature and derivative wrt temperature.
const std::vector< Real > _T_star
Temperature scale for each region.
Water97FluidProperties(const InputParameters &parameters)
virtual void mu_dpT(Real pressure, Real temperature, Real &mu, Real &dmu_dp, Real &dmu_dT) const override
Real dgamma5_dtau(Real pi, Real tau) const
Derivative of Gibbs free energy in Region 5 wrt tau.
Real temperature_from_ph2a(Real pressure, Real enthalpy) const
Backwards equation T(p, h) in Region 2a Eq.
const std::vector< int > _J1
Real gamma5(Real pi, Real tau) const
Gibbs free energy in Region 5.
Real tempXY(Real pressure, subregionEnum xy) const
Boundaries between subregions in region 3.
const std::string density
Definition: NS.h:15
unsigned int subregion2ph(Real pressure, Real enthalpy) const
Provides the correct subregion index for a (P,h) point in region 2.
const std::vector< int > _Jph3b
const std::vector< std::vector< Real > > _n3s
Constants for all 26 subregions in region 3.
const std::vector< Real > _p_star
Pressure scale for each region.
virtual void k_dpT(Real pressure, Real temperature, Real &k, Real &dk_dp, Real &dk_dT) const override
Thermal conductivity and its derivatives wrt pressure and temperature.
virtual void rho_e_dpT(Real pressure, Real temperature, Real &rho, Real &drho_dp, Real &drho_dT, Real &e, Real &de_dp, Real &de_dT) const override
Density and internal energy and their derivatives wrt pressure and temperature.
Real temperature_from_ph2b(Real pressure, Real enthalpy) const
Backwards equation T(p, h) in Region 2b Eq.
unsigned int subregion3ph(Real pressure, Real enthalpy) const
Provides the correct subregion index for a (P,h) point in region 3.
const std::string temperature
Definition: NS.h:25
Real dgamma2_dpi(Real pi, Real tau) const
Derivative of Gibbs free energy in Region 2 wrt pi.
virtual void h_dpT(Real pressure, Real temperature, Real &h, Real &dh_dp, Real &dh_dT) const override
Enthalpy and its derivatives wrt pressure and temperature.
Real b3ab(Real pressure) const
Boundary between subregions a and b in region 3.
const std::vector< int > _J2
const Real _T_triple
Triple point temperature (K)
virtual Real triplePointPressure() const
Water triple point pressure.
const std::string enthalpy
Definition: NS.h:26
virtual Real rho(Real pressure, Real temperature) const override
Density.
Real d2phi3_ddelta2(Real delta, Real tau) const
Second derivative of Helmholtz free energy in Region 3 wrt delta.
const std::vector< int > _Jph2a
const std::vector< Real > _nph2a
Real b23p(Real temperature) const
Auxillary equation for the boundary between regions 2 and 3.
Common class for single phase fluid properties using a pressure and temperature formulation.
const std::vector< int > _Iph1
subregionEnum
Enum of subregion ids for region 3.
InputParameters validParams< SinglePhaseFluidPropertiesPT >()
const std::vector< int > _Jph3a
const std::vector< std::vector< int > > _I3s
Real d2phi3_dtau2(Real delta, Real tau) const
Second derivative of Helmholtz free energy in Region 3 wrt tau.
Real temperature_from_ph3b(Real pressure, Real enthalpy) const
Backwards equation T(p, h) in Region 3b Eq.
Real b23T(Real pressure) const
Auxillary equation for the boundary between regions 2 and 3.
Real densityRegion3(Real pressure, Real temperature) const
Density function for Region 3 - supercritical water and steam.
Real dphi3_dtau(Real delta, Real tau) const
Derivative of Helmholtz free energy in Region 3 wrt tau.
virtual void mu_drhoT_from_rho_T(Real density, Real temperature, Real ddensity_dT, Real &mu, Real &dmu_drho, Real &dmu_dT) const override
Dynamic viscosity and its derivatives wrt density and temperature.
InputParameters validParams< Water97FluidProperties >()
const std::vector< int > _J5
Real subregionVolume(Real pi, Real theta, Real a, Real b, Real c, Real d, Real e, unsigned int sid) const
Specific volume in all subregions of region 3 EXCEPT subregion n (13).
const std::vector< Real > _n3
Constants for region 3.
Real d2gamma2_dpitau(Real pi, Real tau) const
Second derivative of Gibbs free energy in Region 2 wrt pi and tau.
const std::vector< Real > _nph2b
const std::vector< Real > _nph1
virtual Real mu(Real pressure, Real temperature) const override
Real dgamma1_dpi(Real pi, Real tau) const
Derivative of Gibbs free energy in Region 1 wrt pi.
const std::vector< Real > _n1
Reference constants used in to calculate thermophysical properties of water.
unsigned int inRegion(Real pressure, Real temperature) const
Determines the phase region that the given pressure and temperature values lie in.
const std::vector< unsigned int > _par3N
virtual Real molarMass() const override
Water molar mass.
Real gamma2(Real pi, Real tau) const
Gibbs free energy in Region 2 - superheated steam.
Real vaporPressure(Real temperature) const
Saturation pressure as a function of temperature.
virtual Real h(Real pressure, Real temperature) const override
Specific enthalpy.
Real temperature_from_ph2c(Real pressure, Real enthalpy) const
Backwards equation T(p, h) in Region 2c Eq.
const std::vector< int > _Iph2c
Real d2gamma5_dpi2(Real pi, Real tau) const
Second derivative of Gibbs free energy in Region 5 wrt pi.
virtual Real mu_from_rho_T(Real density, Real temperature) const override
const std::vector< Real > _n5
const std::vector< int > _J3
const std::vector< int > _Iph3b
const std::vector< Real > _nph2c
Real dgamma2_dtau(Real pi, Real tau) const
Derivative of Gibbs free energy in Region 2 wrt tau.
Real d2phi3_ddeltatau(Real delta, Real tau) const
Second derivative of Helmholtz free energy in Region 3 wrt delta and tau.
Real temperature_from_ph3a(Real pressure, Real enthalpy) const
Backwards equation T(p, h) in Region 3a Eq.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
Real vaporTemperature(Real pressure) const
Saturation temperature as a function of pressure.
virtual void rho_dpT(Real pressure, Real temperature, Real &rho, Real &drho_dp, Real &drho_dT) const override
Density and its derivatives wrt pressure and temperature.
const std::vector< int > _I3
virtual Real criticalDensity() const
Water critical density.
Real d2gamma5_dpitau(Real pi, Real tau) const
Second derivative of Gibbs free energy in Region 5 wrt pi and tau.
Real gamma1(Real pi, Real tau) const
Gibbs free energy in Region 1 - single phase liquid region.
const Real _rho_critical
Critical density (kg/m^3)
Real phi3(Real delta, Real tau) const
Helmholtz free energy in Region 3.
Real d2gamma1_dpi2(Real pi, Real tau) const
Second derivative of Gibbs free energy in Region 1 wrt pi.
virtual Real s(Real pressure, Real temperature) const override
Specific entropy.
const std::vector< Real > _nph3b
virtual Real triplePointTemperature() const
Water triple point temperature.
Real d2gamma2_dpi2(Real pi, Real tau) const
Second derivative of Gibbs free energy in Region 2 wrt pi.
const std::string pressure
Definition: NS.h:24
const std::vector< int > _J05
Constants for region 5.
virtual void e_dpT(Real pressure, Real temperature, Real &e, Real &de_dp, Real &de_dT) const override
Internal energy and its derivatives wrt pressure and temperature.
Real dgamma1_dtau(Real pi, Real tau) const
Derivative of Gibbs free energy in Region 1 wrt tau.
const std::vector< int > _Jph1
virtual Real k(Real pressure, Real temperature) const override
Thermal conductivity.
virtual Real temperature_from_ph(Real pressure, Real enthalpy) const
Backwards equation T(p, h) From Revised Release on the IAPWS Industrial Formulation 1997 for the Ther...
Real d2gamma2_dtau2(Real pi, Real tau) const
Second derivative of Gibbs free energy in Region 2 wrt tau.
const std::vector< std::vector< int > > _J3s
const std::string internal_energy
Definition: NS.h:28
const std::vector< Real > _nph3a
const Real _Rw
Specific gas constant for H2O (universal gas constant / molar mass of water - J/kg/K) ...
const std::vector< int > _J02
virtual Real criticalTemperature() const
Water critical temperature.
Real d2gamma1_dpitau(Real pi, Real tau) const
Second derivative of Gibbs free energy in Region 1 wrt pi and tau.
const Real _p_triple
Triple point pressure (Pa)
Real b2bc(Real pressure) const
Boundary between subregions b and c in region 2.
Real d2gamma1_dtau2(Real pi, Real tau) const
Second derivative of Gibbs free energy in Region 1 wrt tau.
const std::vector< Real > _n02
Constants for region 2.
virtual Real henryConstant(Real temperature) const override
Henry&#39;s law constant for dissolution in water.
virtual Real k_from_rho_T(Real density, Real temperature) const override
Thermal conductivity as a function of density and temperature.
const Real _Mh2o
Water molar mass (kg/mol)
virtual Real e(Real pressure, Real temperature) const override
Internal energy.
const std::vector< int > _Iph2a
Real dgamma5_dpi(Real pi, Real tau) const
Derivative of Gibbs free energy in Region 5 wrt pi.
const std::vector< Real > _n2