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