Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 : InputParameters
17 254 : Water97FluidProperties::validParams()
18 : {
19 254 : InputParameters params = SinglePhaseFluidProperties::validParams();
20 254 : params.addClassDescription("Fluid properties for water and steam (H2O) using IAPWS-IF97");
21 254 : return params;
22 0 : }
23 :
24 130 : Water97FluidProperties::Water97FluidProperties(const InputParameters & parameters)
25 : : SinglePhaseFluidProperties(parameters),
26 130 : _Mh2o(18.015e-3),
27 130 : _Rw(461.526),
28 130 : _p_critical(22.064e6),
29 130 : _T_critical(647.096),
30 130 : _rho_critical(322.0),
31 130 : _p_triple(611.657),
32 13780 : _T_triple(273.16)
33 : {
34 780 : }
35 :
36 260 : Water97FluidProperties::~Water97FluidProperties() {}
37 :
38 : std::string
39 24 : Water97FluidProperties::fluidName() const
40 : {
41 24 : return "water";
42 : }
43 :
44 : Real
45 53 : Water97FluidProperties::molarMass() const
46 : {
47 53 : return _Mh2o;
48 : }
49 :
50 : Real
51 1 : Water97FluidProperties::criticalPressure() const
52 : {
53 1 : return _p_critical;
54 : }
55 :
56 : Real
57 1 : Water97FluidProperties::criticalTemperature() const
58 : {
59 1 : return _T_critical;
60 : }
61 :
62 : Real
63 1 : Water97FluidProperties::criticalDensity() const
64 : {
65 1 : return _rho_critical;
66 : }
67 :
68 : Real
69 1 : Water97FluidProperties::triplePointPressure() const
70 : {
71 1 : return _p_triple;
72 : }
73 :
74 : Real
75 1 : Water97FluidProperties::triplePointTemperature() const
76 : {
77 1 : return _T_triple;
78 : }
79 :
80 : Real
81 220942 : Water97FluidProperties::rho_from_p_T(Real pressure, Real temperature) const
82 : {
83 220942 : return rho_from_p_T_template(pressure, temperature);
84 : }
85 :
86 : ADReal
87 9 : Water97FluidProperties::rho_from_p_T(const ADReal & pressure, const ADReal & temperature) const
88 : {
89 9 : return rho_from_p_T_template(pressure, temperature);
90 : }
91 :
92 : void
93 20 : Water97FluidProperties::rho_from_p_T(
94 : Real pressure, Real temperature, Real & rho, Real & drho_dp, Real & drho_dT) const
95 : {
96 : rho_from_p_T_template(pressure, temperature, rho, drho_dp, drho_dT);
97 20 : }
98 :
99 : void
100 0 : Water97FluidProperties::rho_from_p_T(const ADReal & pressure,
101 : const ADReal & temperature,
102 : ADReal & rho,
103 : ADReal & drho_dp,
104 : ADReal & drho_dT) const
105 : {
106 0 : rho_from_p_T_template(pressure, temperature, rho, drho_dp, drho_dT);
107 0 : }
108 :
109 : Real
110 0 : Water97FluidProperties::v_from_p_T(Real pressure, Real temperature) const
111 : {
112 0 : return v_from_p_T_template(pressure, temperature);
113 : }
114 :
115 : ADReal
116 0 : Water97FluidProperties::v_from_p_T(const ADReal & pressure, const ADReal & temperature) const
117 : {
118 0 : return v_from_p_T_template(pressure, temperature);
119 : }
120 :
121 : void
122 4 : Water97FluidProperties::v_from_p_T(
123 : Real pressure, Real temperature, Real & v, Real & dv_dp, Real & dv_dT) const
124 : {
125 : v_from_p_T_template(pressure, temperature, v, dv_dp, dv_dT);
126 4 : }
127 :
128 : void
129 12 : Water97FluidProperties::v_from_p_T(const ADReal & pressure,
130 : const ADReal & temperature,
131 : ADReal & v,
132 : ADReal & dv_dp,
133 : ADReal & dv_dT) const
134 : {
135 12 : v_from_p_T_template(pressure, temperature, v, dv_dp, dv_dT);
136 12 : }
137 :
138 : Real
139 7 : Water97FluidProperties::p_from_v_e(const Real v, const Real e) const
140 : {
141 7 : return p_from_v_e_template(v, e);
142 : }
143 :
144 : ADReal
145 15 : Water97FluidProperties::p_from_v_e(const ADReal & v, const ADReal & e) const
146 : {
147 15 : return p_from_v_e_template(v, e);
148 : }
149 :
150 : Real
151 3 : Water97FluidProperties::e_from_p_rho(Real p, Real rho) const
152 : {
153 3 : return e_from_p_rho_template(p, rho);
154 : }
155 :
156 : ADReal
157 3 : Water97FluidProperties::e_from_p_rho(const ADReal & p, const ADReal & rho) const
158 : {
159 3 : return e_from_p_rho_template(p, rho);
160 : }
161 :
162 : void
163 15 : Water97FluidProperties::e_from_p_rho(
164 : const Real p, const Real rho, Real & e, Real & de_dp, Real & de_drho) const
165 : {
166 : e_from_p_rho_template(p, rho, e, de_dp, de_drho);
167 15 : }
168 :
169 : void
170 31 : Water97FluidProperties::e_from_p_rho(
171 : const ADReal & p, const ADReal & rho, ADReal & e, ADReal & de_dp, ADReal & de_drho) const
172 : {
173 31 : e_from_p_rho_template(p, rho, e, de_dp, de_drho);
174 31 : }
175 :
176 : Real
177 220258 : Water97FluidProperties::e_from_p_T(Real pressure, Real temperature) const
178 : {
179 220258 : return e_from_p_T_template(pressure, temperature);
180 : }
181 :
182 : ADReal
183 6 : Water97FluidProperties::e_from_p_T(const ADReal & pressure, const ADReal & temperature) const
184 : {
185 6 : return e_from_p_T_template(pressure, temperature);
186 : }
187 :
188 : void
189 7 : Water97FluidProperties::e_from_p_T(
190 : Real pressure, Real temperature, Real & e, Real & de_dp, Real & de_dT) const
191 : {
192 7 : e_from_p_T_template(pressure, temperature, e, de_dp, de_dT);
193 7 : }
194 :
195 : void
196 1 : Water97FluidProperties::e_from_p_T(const ADReal & pressure,
197 : const ADReal & temperature,
198 : ADReal & e,
199 : ADReal & de_dp,
200 : ADReal & de_dT) const
201 : {
202 1 : e_from_p_T_template(pressure, temperature, e, de_dp, de_dT);
203 1 : }
204 :
205 : ADReal
206 2 : Water97FluidProperties::e_from_v_h(const ADReal & v, const ADReal & h) const
207 : {
208 2 : const auto [p, T] = p_T_from_v_h(v, h);
209 2 : return e_from_p_T(p, T);
210 : }
211 :
212 : Real
213 1 : Water97FluidProperties::T_from_v_e(const Real v, const Real e) const
214 : {
215 1 : return p_T_from_v_e(v, e).second;
216 : }
217 :
218 : ADReal
219 1 : Water97FluidProperties::T_from_v_e(const ADReal & v, const ADReal & e) const
220 : {
221 1 : return p_T_from_v_e(v, e).second;
222 : }
223 :
224 : ADReal
225 2 : Water97FluidProperties::c_from_v_e(const ADReal & v, const ADReal & e) const
226 : {
227 2 : const auto [p, T] = p_T_from_v_e(v, e);
228 2 : return c_from_p_T(p, T);
229 : }
230 :
231 : Real
232 231 : Water97FluidProperties::c_from_p_T(const Real p, const Real T) const
233 : {
234 231 : return c_from_p_T_template(p, T);
235 : }
236 :
237 : ADReal
238 5 : Water97FluidProperties::c_from_p_T(const ADReal & p, const ADReal & T) const
239 : {
240 5 : return c_from_p_T_template(p, T);
241 : }
242 :
243 : Real
244 666 : Water97FluidProperties::cp_from_p_T(Real pressure, Real temperature) const
245 : {
246 666 : return cp_from_p_T_template(pressure, temperature);
247 : }
248 :
249 : ADReal
250 5 : Water97FluidProperties::cp_from_p_T(const ADReal & pressure, const ADReal & temperature) const
251 : {
252 5 : return cp_from_p_T_template(pressure, temperature);
253 : }
254 :
255 : ADReal
256 2 : Water97FluidProperties::cp_from_v_e(const ADReal & v, const ADReal & e) const
257 : {
258 2 : const auto [p, T] = p_T_from_v_e(v, e);
259 2 : return cp_from_p_T(p, T);
260 : }
261 :
262 : Real
263 223 : Water97FluidProperties::cv_from_p_T(Real pressure, Real temperature) const
264 : {
265 223 : return cv_from_p_T_template(pressure, temperature);
266 : }
267 :
268 : ADReal
269 5 : Water97FluidProperties::cv_from_p_T(const ADReal & pressure, const ADReal & temperature) const
270 : {
271 5 : return cv_from_p_T_template(pressure, temperature);
272 : }
273 :
274 : ADReal
275 2 : Water97FluidProperties::cv_from_v_e(const ADReal & v, const ADReal & e) const
276 : {
277 2 : const auto [p, T] = p_T_from_v_e(v, e);
278 2 : return cv_from_p_T(p, T);
279 : }
280 :
281 : Real
282 220244 : Water97FluidProperties::mu_from_p_T(Real pressure, Real temperature) const
283 : {
284 220244 : return mu_from_p_T_template(pressure, temperature);
285 : }
286 :
287 : ADReal
288 3 : Water97FluidProperties::mu_from_p_T(const ADReal & pressure, const ADReal & temperature) const
289 : {
290 3 : return mu_from_p_T_template(pressure, temperature);
291 : }
292 :
293 : void
294 6 : Water97FluidProperties::mu_from_p_T(
295 : Real pressure, Real temperature, Real & mu, Real & dmu_dp, Real & dmu_dT) const
296 : {
297 : Real rho, drho_dp, drho_dT;
298 6 : this->rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
299 : Real dmu_drho;
300 6 : this->mu_from_rho_T(rho, temperature, drho_dT, mu, dmu_drho, dmu_dT);
301 6 : dmu_dp = dmu_drho * drho_dp;
302 6 : }
303 :
304 : Real
305 17 : Water97FluidProperties::mu_from_rho_T(Real density, Real temperature) const
306 : {
307 17 : return mu_from_rho_T_template(density, temperature);
308 : }
309 :
310 : void
311 9 : Water97FluidProperties::mu_from_rho_T(Real density,
312 : Real temperature,
313 : Real ddensity_dT,
314 : Real & mu,
315 : Real & dmu_drho,
316 : Real & dmu_dT) const
317 : {
318 : const Real mu_star = 1.0e-6;
319 9 : const Real rhobar = density / _rho_critical;
320 9 : const Real Tbar = temperature / _T_critical;
321 9 : const Real drhobar_drho = 1.0 / _rho_critical;
322 9 : const Real dTbar_dT = 1.0 / _T_critical;
323 :
324 : // Limit of zero density. Derivative wrt rho is 0
325 : Real sum0 = 0.0, dsum0_dTbar = 0.0;
326 45 : for (std::size_t i = 0; i < _mu_H0.size(); ++i)
327 : {
328 36 : sum0 += _mu_H0[i] / MathUtils::pow(Tbar, i);
329 36 : dsum0_dTbar -= i * _mu_H0[i] / MathUtils::pow(Tbar, i + 1);
330 : }
331 :
332 9 : const Real mu0 = 100.0 * std::sqrt(Tbar) / sum0;
333 : const Real dmu0_dTbar =
334 9 : 50.0 / std::sqrt(Tbar) / sum0 - 100.0 * std::sqrt(Tbar) * dsum0_dTbar / sum0 / sum0;
335 :
336 : // Residual component due to finite density
337 : Real sum1 = 0.0, dsum1_drhobar = 0.0, dsum1_dTbar = 0.0;
338 63 : for (unsigned int i = 0; i < 6; ++i)
339 : {
340 54 : const Real fact = MathUtils::pow(1.0 / Tbar - 1.0, i);
341 54 : const Real dfact_dTbar = i * MathUtils::pow(1.0 / Tbar - 1.0, i - 1) / Tbar / Tbar;
342 :
343 432 : for (unsigned int j = 0; j < 7; ++j)
344 : {
345 378 : sum1 += fact * _mu_Hij[i][j] * MathUtils::pow(rhobar - 1.0, j);
346 378 : dsum1_dTbar -= dfact_dTbar * _mu_Hij[i][j] * MathUtils::pow(rhobar - 1.0, j);
347 378 : dsum1_drhobar += j * fact * _mu_Hij[i][j] * MathUtils::pow(rhobar - 1.0, j - 1);
348 : }
349 : }
350 :
351 9 : const Real mu1 = std::exp(rhobar * sum1);
352 9 : const Real dmu1_drhobar = (sum1 + rhobar * dsum1_drhobar) * mu1;
353 9 : const Real dmu1_dTbar = (rhobar * dsum1_dTbar) * mu1;
354 :
355 : // Viscosity and its derivatives are then
356 9 : mu = mu_star * mu0 * mu1;
357 9 : dmu_drho = mu_star * mu0 * dmu1_drhobar * drhobar_drho;
358 9 : dmu_dT = mu_star * (dmu0_dTbar * mu1 + mu0 * dmu1_dTbar) * dTbar_dT + dmu_drho * ddensity_dT;
359 9 : }
360 :
361 : ADReal
362 2 : Water97FluidProperties::mu_from_v_e(const ADReal & v, const ADReal & e) const
363 : {
364 2 : const auto [rho, T] = rho_T_from_v_e(v, e);
365 2 : return mu_from_rho_T_template(rho, T);
366 : }
367 :
368 : void
369 1 : Water97FluidProperties::rho_mu_from_p_T(Real pressure,
370 : Real temperature,
371 : Real & rho,
372 : Real & mu) const
373 : {
374 1 : rho = this->rho_from_p_T(pressure, temperature);
375 1 : mu = this->mu_from_rho_T(rho, temperature);
376 1 : }
377 :
378 : void
379 1 : Water97FluidProperties::rho_mu_from_p_T(Real pressure,
380 : Real temperature,
381 : Real & rho,
382 : Real & drho_dp,
383 : Real & drho_dT,
384 : Real & mu,
385 : Real & dmu_dp,
386 : Real & dmu_dT) const
387 : {
388 1 : this->rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
389 : Real dmu_drho;
390 1 : this->mu_from_rho_T(rho, temperature, drho_dT, mu, dmu_drho, dmu_dT);
391 1 : dmu_dp = dmu_drho * drho_dp;
392 1 : }
393 :
394 : Real
395 229 : Water97FluidProperties::k_from_p_T(Real pressure, Real temperature) const
396 : {
397 229 : return k_from_p_T_template(pressure, temperature);
398 : }
399 :
400 : ADReal
401 4 : Water97FluidProperties::k_from_p_T(const ADReal & pressure, const ADReal & temperature) const
402 : {
403 4 : return k_from_p_T_template(pressure, temperature);
404 : }
405 :
406 : void
407 0 : Water97FluidProperties::k_from_p_T(Real, Real, Real &, Real &, Real &) const
408 : {
409 0 : mooseError("k_from_p_T() is not implemented.");
410 : }
411 :
412 : Real
413 0 : Water97FluidProperties::k_from_rho_T(Real density, Real temperature) const
414 : {
415 0 : return k_from_rho_T_template(density, temperature);
416 : }
417 :
418 : Real
419 4 : Water97FluidProperties::k_from_v_e(Real v, Real e) const
420 : {
421 4 : return k_from_v_e_template(v, e);
422 : }
423 :
424 : ADReal
425 4 : Water97FluidProperties::k_from_v_e(const ADReal & v, const ADReal & e) const
426 : {
427 4 : return k_from_v_e_template(v, e);
428 : }
429 :
430 : Real
431 263 : Water97FluidProperties::s_from_p_T(Real pressure, Real temperature) const
432 : {
433 263 : return s_from_p_T_template(pressure, temperature);
434 : }
435 :
436 : ADReal
437 1 : Water97FluidProperties::s_from_p_T(const ADReal & pressure, const ADReal & temperature) const
438 : {
439 1 : return s_from_p_T_template(pressure, temperature);
440 : }
441 :
442 : void
443 7 : Water97FluidProperties::s_from_p_T(
444 : const Real pressure, const Real temperature, Real & s, Real & ds_dp, Real & ds_dT) const
445 : {
446 : s_from_p_T_template(pressure, temperature, s, ds_dp, ds_dT);
447 7 : }
448 :
449 : Real
450 14 : Water97FluidProperties::s_from_h_p(Real enthalpy, Real pressure) const
451 : {
452 14 : Real T = T_from_p_h(pressure, enthalpy);
453 14 : return s_from_p_T(pressure, T);
454 : }
455 :
456 : ADReal
457 0 : Water97FluidProperties::s_from_h_p(const ADReal & enthalpy, const ADReal & pressure) const
458 : {
459 0 : ADReal temperature = T_from_p_h_ad(pressure, enthalpy);
460 0 : return s_from_p_T_template(pressure, temperature);
461 : }
462 :
463 : void
464 1 : Water97FluidProperties::s_from_h_p(
465 : const Real enthalpy, const Real pressure, Real & s, Real & ds_dh, Real & ds_dp) const
466 : {
467 1 : ADReal p = pressure;
468 : Moose::derivInsert(p.derivatives(), 0, 1.0);
469 :
470 1 : ADReal h = enthalpy;
471 : Moose::derivInsert(h.derivatives(), 1, 1.0);
472 :
473 1 : ADReal T = T_from_p_h_ad(p, h);
474 1 : ADReal entropy = s_from_p_T_template(p, T);
475 :
476 1 : ds_dh = entropy.derivatives()[1];
477 1 : ds_dp = entropy.derivatives()[0];
478 1 : s = entropy.value();
479 1 : }
480 :
481 : void
482 3 : Water97FluidProperties::s_from_p_T(const ADReal & pressure,
483 : const ADReal & temperature,
484 : ADReal & s,
485 : ADReal & ds_dp,
486 : ADReal & ds_dT) const
487 : {
488 3 : s_from_p_T_template(pressure, temperature, s, ds_dp, ds_dT);
489 3 : }
490 :
491 : Real
492 220852 : Water97FluidProperties::h_from_p_T(Real pressure, Real temperature) const
493 : {
494 220852 : return h_from_p_T_template(pressure, temperature);
495 : }
496 :
497 : ADReal
498 1 : Water97FluidProperties::h_from_p_T(const ADReal & pressure, const ADReal & temperature) const
499 : {
500 1 : return h_from_p_T_template(pressure, temperature);
501 : }
502 :
503 : void
504 15 : Water97FluidProperties::h_from_p_T(
505 : Real pressure, Real temperature, Real & h, Real & dh_dp, Real & dh_dT) const
506 : {
507 : h_from_p_T_template(pressure, temperature, h, dh_dp, dh_dT);
508 15 : }
509 :
510 : void
511 15 : Water97FluidProperties::h_from_p_T(const ADReal & pressure,
512 : const ADReal & temperature,
513 : ADReal & h,
514 : ADReal & dh_dp,
515 : ADReal & dh_dT) const
516 : {
517 15 : h_from_p_T_template(pressure, temperature, h, dh_dp, dh_dT);
518 15 : }
519 :
520 : Real
521 1768248 : Water97FluidProperties::vaporPressure(Real temperature) const
522 : {
523 : // Check whether the input temperature is within the region of validity of this equation.
524 : // Valid for 273.15 K <= t <= 647.096 K
525 1768248 : if (temperature < 273.15 || temperature > _T_critical)
526 0 : mooseException(name(),
527 : ": vaporPressure(): Temperature ",
528 : temperature,
529 : " is outside range 273.15 K <= T "
530 : "<= 647.096 K");
531 :
532 : Real theta, theta2, a, b, c, p;
533 1768248 : theta = temperature + _n4[8] / (temperature - _n4[9]);
534 1768248 : theta2 = theta * theta;
535 1768248 : a = theta2 + _n4[0] * theta + _n4[1];
536 1768248 : b = _n4[2] * theta2 + _n4[3] * theta + _n4[4];
537 1768248 : c = _n4[5] * theta2 + _n4[6] * theta + _n4[7];
538 1768248 : p = Utility::pow<4>(2.0 * c / (-b + std::sqrt(b * b - 4.0 * a * c)));
539 :
540 1768248 : return p * 1.e6;
541 : }
542 :
543 : ADReal
544 60 : Water97FluidProperties::vaporTemperature_ad(const ADReal & pressure) const
545 : {
546 : // Check whether the input pressure is within the region of validity of this equation.
547 : // Valid for 611.213 Pa <= p <= 22.064 MPa
548 60 : if (pressure.value() < 611.23 || pressure.value() > _p_critical)
549 0 : mooseException(name() + ": vaporTemperature(): Pressure ",
550 : pressure.value(),
551 : " is outside range 611.213 Pa <= p <= 22.064 MPa");
552 :
553 120 : const ADReal beta = std::pow(pressure / 1.e6, 0.25);
554 : const ADReal beta2 = beta * beta;
555 60 : const ADReal e = beta2 + _n4[2] * beta + _n4[5];
556 60 : const ADReal f = _n4[0] * beta2 + _n4[3] * beta + _n4[6];
557 60 : const ADReal g = _n4[1] * beta2 + _n4[4] * beta + _n4[7];
558 180 : const ADReal d = 2.0 * g / (-f - std::sqrt(f * f - 4.0 * e * g));
559 :
560 240 : return (_n4[9] + d - std::sqrt((_n4[9] + d) * (_n4[9] + d) - 4.0 * (_n4[8] + _n4[9] * d))) / 2.0;
561 : }
562 :
563 : Real
564 60 : Water97FluidProperties::vaporTemperature(Real pressure) const
565 : {
566 60 : const ADReal p = pressure;
567 :
568 60 : return vaporTemperature_ad(p).value();
569 : }
570 :
571 : void
572 0 : Water97FluidProperties::vaporTemperature(Real pressure, Real & Tsat, Real & dTsat_dp) const
573 : {
574 0 : ADReal p = pressure;
575 : Moose::derivInsert(p.derivatives(), 0, 1.0);
576 :
577 0 : const ADReal T = vaporTemperature_ad(p);
578 :
579 0 : Tsat = T.value();
580 0 : dTsat_dp = T.derivatives()[0];
581 0 : }
582 :
583 : Real
584 103 : Water97FluidProperties::b23p(Real temperature) const
585 : {
586 : // Check whether the input temperature is within the region of validity of this equation.
587 : // Valid for 623.15 K <= t <= 863.15 K
588 103 : if (temperature < 623.15 || temperature > 863.15)
589 0 : mooseException(name(),
590 : ": b23p(): Temperature ",
591 : temperature,
592 : " is outside range of 623.15 K <= T <= 863.15 K");
593 :
594 103 : return (_n23[0] + _n23[1] * temperature + _n23[2] * temperature * temperature) * 1.e6;
595 : }
596 :
597 : Real
598 20 : Water97FluidProperties::b23T(Real pressure) const
599 : {
600 : // Check whether the input pressure is within the region of validity of this equation.
601 : // Valid for 16.529 MPa <= p <= 100 MPa
602 20 : if (pressure < 16.529e6 || pressure > 100.0e6)
603 0 : mooseException(
604 : name(), ": b23T(): Pressure ", pressure, "is outside range 16.529 MPa <= p <= 100 MPa");
605 :
606 20 : return _n23[3] + std::sqrt((pressure / 1.e6 - _n23[4]) / _n23[2]);
607 : }
608 :
609 : unsigned int
610 39 : Water97FluidProperties::inRegionPH(Real pressure, Real enthalpy) const
611 : {
612 : unsigned int region;
613 :
614 : // Need to calculate enthalpies at the boundaries to delineate regions
615 39 : Real p273 = vaporPressure(273.15);
616 39 : Real p623 = vaporPressure(623.15);
617 :
618 39 : if (pressure >= p273 && pressure <= p623)
619 : {
620 42 : if (enthalpy >= h_from_p_T(pressure, 273.15) &&
621 21 : enthalpy <= h_from_p_T(pressure, vaporTemperature(pressure)))
622 : region = 1;
623 24 : else if (enthalpy > h_from_p_T(pressure, vaporTemperature(pressure)) &&
624 12 : enthalpy <= h_from_p_T(pressure, 1073.15))
625 : region = 2;
626 0 : else if (enthalpy > h_from_p_T(pressure, 1073.15) && enthalpy <= h_from_p_T(pressure, 2273.15))
627 : region = 5;
628 : else
629 0 : mooseException("Enthalpy ", enthalpy, " is out of range in ", name(), ": inRegionPH()");
630 : }
631 18 : else if (pressure > p623 && pressure <= 50.0e6)
632 : {
633 9 : if (enthalpy >= h_from_p_T(pressure, 273.15) && enthalpy <= h_from_p_T(pressure, 623.15))
634 : region = 1;
635 18 : else if (enthalpy > h_from_p_T(pressure, 623.15) &&
636 9 : enthalpy <= h_from_p_T(pressure, b23T(pressure)))
637 : region = 3;
638 6 : else if (enthalpy > h_from_p_T(pressure, b23T(pressure)) &&
639 3 : enthalpy <= h_from_p_T(pressure, 1073.15))
640 : region = 2;
641 0 : else if (enthalpy > h_from_p_T(pressure, 1073.15) && enthalpy <= h_from_p_T(pressure, 2273.15))
642 : region = 5;
643 : else
644 0 : mooseException("Enthalpy ", enthalpy, " is out of range in ", name(), ": inRegionPH()");
645 : }
646 9 : else if (pressure > 50.0e6 && pressure <= 100.0e6)
647 : {
648 9 : if (enthalpy >= h_from_p_T(pressure, 273.15) && enthalpy <= h_from_p_T(pressure, 623.15))
649 : region = 1;
650 10 : else if (enthalpy > h_from_p_T(pressure, 623.15) &&
651 5 : enthalpy <= h_from_p_T(pressure, b23T(pressure)))
652 : region = 3;
653 4 : else if (enthalpy > h_from_p_T(pressure, b23T(pressure)) &&
654 2 : enthalpy <= h_from_p_T(pressure, 1073.15))
655 : region = 2;
656 : else
657 0 : mooseException("Enthalpy ", enthalpy, " is out of range in ", name(), ": inRegionPH()");
658 : }
659 : else
660 0 : mooseException("Pressure ", pressure, " is out of range in ", name(), ": inRegionPH()");
661 :
662 39 : return region;
663 : }
664 :
665 : unsigned int
666 17 : Water97FluidProperties::subregion2ph(Real pressure, Real enthalpy) const
667 : {
668 : unsigned int subregion;
669 :
670 17 : if (pressure <= 4.0e6)
671 : subregion = 1;
672 7 : else if (pressure > 4.0e6 && pressure < 6.5467e6)
673 : subregion = 2;
674 : else
675 : {
676 5 : if (enthalpy >= b2bc(pressure))
677 : subregion = 2;
678 : else
679 : subregion = 3;
680 : }
681 :
682 17 : return subregion;
683 : }
684 :
685 : unsigned int
686 9 : Water97FluidProperties::subregion3ph(Real pressure, Real enthalpy) const
687 : {
688 : unsigned int subregion;
689 :
690 9 : if (enthalpy <= b3ab(pressure))
691 : subregion = 1;
692 : else
693 : subregion = 2;
694 :
695 9 : return subregion;
696 : }
697 :
698 : Real
699 6 : Water97FluidProperties::b2bc(Real pressure) const
700 : {
701 : // Check whether the input pressure is within the region of validity of this equation.
702 6 : if (pressure < 6.5467e6 || pressure > 100.0e6)
703 0 : mooseException(
704 : name(), ": b2bc(): Pressure ", pressure, " is outside range of 6.5467 MPa <= p <= 100 MPa");
705 :
706 6 : Real pi = pressure / 1.0e6;
707 :
708 6 : return (0.26526571908428e4 + std::sqrt((pi - 0.45257578905948e1) / 0.12809002730136e-3)) * 1.0e3;
709 : }
710 :
711 : Real
712 36 : Water97FluidProperties::T_from_p_h(Real pressure, Real enthalpy) const
713 : {
714 36 : const ADReal p = pressure;
715 36 : const ADReal h = enthalpy;
716 :
717 36 : return T_from_p_h_ad(p, h).value();
718 : }
719 :
720 : void
721 0 : Water97FluidProperties::T_from_p_h(
722 : Real pressure, Real enthalpy, Real & temperature, Real & dT_dp, Real & dT_dh) const
723 : {
724 0 : ADReal p = pressure;
725 : Moose::derivInsert(p.derivatives(), 0, 1.0);
726 0 : ADReal h = enthalpy;
727 : Moose::derivInsert(h.derivatives(), 1, 1.0);
728 :
729 0 : const ADReal T = T_from_p_h_ad(p, h);
730 :
731 0 : temperature = T.value();
732 0 : dT_dp = T.derivatives()[0];
733 0 : dT_dh = T.derivatives()[1];
734 0 : }
735 :
736 : ADReal
737 2 : Water97FluidProperties::T_from_p_h(const ADReal & pressure, const ADReal & enthalpy) const
738 : {
739 2 : return T_from_p_h_ad(pressure, enthalpy);
740 : }
741 :
742 : ADReal
743 39 : Water97FluidProperties::T_from_p_h_ad(const ADReal & pressure, const ADReal & enthalpy) const
744 : {
745 39 : ADReal temperature = 0.0;
746 :
747 : // Determine which region the point is in
748 39 : const unsigned int region = inRegionPH(pressure.value(), enthalpy.value());
749 :
750 39 : switch (region)
751 : {
752 13 : case 1:
753 13 : temperature = temperature_from_ph1(pressure, enthalpy);
754 13 : break;
755 :
756 : case 2:
757 : {
758 : // First, determine which subregion the point is in:
759 17 : const unsigned int subregion = subregion2ph(pressure.value(), enthalpy.value());
760 :
761 17 : if (subregion == 1)
762 10 : temperature = temperature_from_ph2a(pressure, enthalpy);
763 7 : else if (subregion == 2)
764 3 : temperature = temperature_from_ph2b(pressure, enthalpy);
765 : else
766 4 : temperature = temperature_from_ph2c(pressure, enthalpy);
767 : break;
768 : }
769 :
770 : case 3:
771 : {
772 : // First, determine which subregion the point is in:
773 9 : const unsigned int subregion = subregion3ph(pressure.value(), enthalpy.value());
774 :
775 9 : if (subregion == 1)
776 4 : temperature = temperature_from_ph3a(pressure, enthalpy);
777 : else
778 5 : temperature = temperature_from_ph3b(pressure, enthalpy);
779 : break;
780 : }
781 :
782 0 : case 5:
783 0 : mooseError("temperature_from_ph() not implemented for region 5");
784 : break;
785 :
786 0 : default:
787 0 : mooseError("inRegionPH() has given an incorrect region");
788 : }
789 :
790 39 : return temperature;
791 : }
792 :
793 : ADReal
794 13 : Water97FluidProperties::temperature_from_ph1(const ADReal & pressure, const ADReal & enthalpy) const
795 : {
796 13 : const ADReal pi = pressure / 1.0e6;
797 13 : const ADReal eta = enthalpy / 2500.0e3;
798 13 : ADReal sum = 0.0;
799 :
800 273 : for (std::size_t i = 0; i < _nph1.size(); ++i)
801 520 : sum += _nph1[i] * std::pow(pi, _Iph1[i]) * std::pow(eta + 1.0, _Jph1[i]);
802 :
803 13 : return sum;
804 : }
805 :
806 : ADReal
807 10 : Water97FluidProperties::temperature_from_ph2a(const ADReal & pressure,
808 : const ADReal & enthalpy) const
809 : {
810 10 : const ADReal pi = pressure / 1.0e6;
811 10 : const ADReal eta = enthalpy / 2000.0e3;
812 10 : ADReal sum = 0.0;
813 :
814 : // Factor out the negative in std::pow(eta - 2.1, _Jph2a[i]) to avoid fpe in dbg (see #13163)
815 10 : const Real sgn = MathUtils::sign(eta.value() - 2.1);
816 :
817 370 : for (std::size_t i = 0; i < _nph2a.size(); ++i)
818 360 : sum += _nph2a[i] * std::pow(pi, _Iph2a[i]) * std::pow(std::abs(eta - 2.1), _Jph2a[i]) *
819 720 : std::pow(sgn, _Jph2a[i]);
820 :
821 10 : return sum;
822 : }
823 :
824 : ADReal
825 3 : Water97FluidProperties::temperature_from_ph2b(const ADReal & pressure,
826 : const ADReal & enthalpy) const
827 : {
828 3 : const ADReal pi = pressure / 1.0e6;
829 3 : const ADReal eta = enthalpy / 2000.0e3;
830 3 : ADReal sum = 0.0;
831 :
832 : // Factor out the negatives in std::pow(pi - 2.0, _Iph2b[i])* std::pow(eta - 2.6, _Jph2b[i])
833 : // to avoid fpe in dbg (see #13163)
834 3 : const Real sgn0 = MathUtils::sign(pi.value() - 2.0);
835 3 : const Real sgn1 = MathUtils::sign(eta.value() - 2.6);
836 :
837 117 : for (std::size_t i = 0; i < _nph2b.size(); ++i)
838 114 : sum += _nph2b[i] * std::pow(std::abs(pi - 2.0), _Iph2b[i]) * std::pow(sgn0, _Iph2b[i]) *
839 228 : std::pow(std::abs(eta - 2.6), _Jph2b[i]) * std::pow(sgn1, _Jph2b[i]);
840 :
841 3 : return sum;
842 : }
843 :
844 : ADReal
845 4 : Water97FluidProperties::temperature_from_ph2c(const ADReal & pressure,
846 : const ADReal & enthalpy) const
847 : {
848 4 : const ADReal pi = pressure / 1.0e6;
849 4 : const ADReal eta = enthalpy / 2000.0e3;
850 4 : ADReal sum = 0.0;
851 :
852 : // Factor out the negative in std::pow(eta - 1.8, _Jph2c[i]) to avoid fpe in dbg (see #13163)
853 4 : const Real sgn = MathUtils::sign(eta.value() - 1.8);
854 :
855 96 : for (std::size_t i = 0; i < _nph2c.size(); ++i)
856 184 : sum += _nph2c[i] * std::pow(pi + 25.0, _Iph2c[i]) * std::pow(std::abs(eta - 1.8), _Jph2c[i]) *
857 184 : std::pow(sgn, _Jph2c[i]);
858 :
859 4 : return sum;
860 : }
861 :
862 : Real
863 10 : Water97FluidProperties::b3ab(Real pressure) const
864 : {
865 : // Check whether the input pressure is within the region of validity of this equation.
866 10 : if (pressure < b23p(623.15) || pressure > 100.0e6)
867 0 : mooseException(
868 : name(), ": b3ab(): Pressure ", pressure, "is outside range of 16.529 MPa <= p <= 100 MPa");
869 :
870 10 : Real pi = pressure / 1.0e6;
871 10 : Real eta = 0.201464004206875e4 + 0.374696550136983e1 * pi - 0.219921901054187e-1 * pi * pi +
872 10 : 0.87513168600995e-4 * pi * pi * pi;
873 :
874 10 : return eta * 1.0e3;
875 : }
876 :
877 : ADReal
878 4 : Water97FluidProperties::temperature_from_ph3a(const ADReal & pressure,
879 : const ADReal & enthalpy) const
880 : {
881 4 : const ADReal pi = pressure / 100.0e6;
882 4 : const ADReal eta = enthalpy / 2300.0e3;
883 4 : ADReal sum = 0.0;
884 :
885 128 : for (std::size_t i = 0; i < _nph3a.size(); ++i)
886 372 : sum += _nph3a[i] * std::pow(pi + 0.24, _Iph3a[i]) * std::pow(eta - 0.615, _Jph3a[i]);
887 :
888 8 : return sum * 760.0;
889 : }
890 :
891 : ADReal
892 5 : Water97FluidProperties::temperature_from_ph3b(const ADReal & pressure,
893 : const ADReal & enthalpy) const
894 : {
895 5 : const ADReal pi = pressure / 100.0e6;
896 5 : const ADReal eta = enthalpy / 2800.0e3;
897 5 : ADReal sum = 0.0;
898 :
899 170 : for (std::size_t i = 0; i < _nph3b.size(); ++i)
900 495 : sum += _nph3b[i] * std::pow(pi + 0.298, _Iph3b[i]) * std::pow(eta - 0.72, _Jph3b[i]);
901 :
902 10 : return sum * 860.0;
903 : }
904 :
905 : Real
906 9 : Water97FluidProperties::henryConstant(Real temperature, const std::vector<Real> & coeffs) const
907 : {
908 9 : const Real A = coeffs[0];
909 9 : const Real B = coeffs[1];
910 9 : const Real C = coeffs[2];
911 :
912 9 : const Real Tr = temperature / 647.096;
913 9 : const Real tau = 1.0 - Tr;
914 :
915 : const Real lnkh =
916 9 : A / Tr + B * std::pow(tau, 0.355) / Tr + C * std::pow(Tr, -0.41) * std::exp(tau);
917 :
918 : // The vapor pressure used in this formulation
919 : const std::vector<Real> a{
920 9 : -7.85951783, 1.84408259, -11.7866497, 22.6807411, -15.9618719, 1.80122502};
921 9 : const std::vector<Real> b{1.0, 1.5, 3.0, 3.5, 4.0, 7.5};
922 : Real sum = 0.0;
923 :
924 63 : for (std::size_t i = 0; i < a.size(); ++i)
925 54 : sum += a[i] * std::pow(tau, b[i]);
926 :
927 9 : return 22.064e6 * std::exp(sum / Tr) * std::exp(lnkh);
928 9 : }
929 :
930 : void
931 1 : Water97FluidProperties::henryConstant(Real temperature,
932 : const std::vector<Real> & coeffs,
933 : Real & Kh,
934 : Real & dKh_dT) const
935 : {
936 1 : const Real A = coeffs[0];
937 1 : const Real B = coeffs[1];
938 1 : const Real C = coeffs[2];
939 :
940 : const Real pc = 22.064e6;
941 : const Real Tc = 647.096;
942 :
943 1 : const Real Tr = temperature / Tc;
944 1 : const Real tau = 1.0 - Tr;
945 :
946 : const Real lnkh =
947 1 : A / Tr + B * std::pow(tau, 0.355) / Tr + C * std::pow(Tr, -0.41) * std::exp(tau);
948 : const Real dlnkh_dT =
949 1 : (-A / Tr / Tr - B * std::pow(tau, 0.355) / Tr / Tr - 0.355 * B * std::pow(tau, -0.645) / Tr -
950 1 : 0.41 * C * std::pow(Tr, -1.41) * std::exp(tau) - C * std::pow(Tr, -0.41) * std::exp(tau)) /
951 1 : Tc;
952 :
953 : // The vapor pressure used in this formulation
954 : const std::vector<Real> a{
955 1 : -7.85951783, 1.84408259, -11.7866497, 22.6807411, -15.9618719, 1.80122502};
956 1 : const std::vector<Real> b{1.0, 1.5, 3.0, 3.5, 4.0, 7.5};
957 : Real sum = 0.0;
958 : Real dsum = 0.0;
959 :
960 7 : for (std::size_t i = 0; i < a.size(); ++i)
961 : {
962 6 : sum += a[i] * std::pow(tau, b[i]);
963 6 : dsum += a[i] * b[i] * std::pow(tau, b[i] - 1.0);
964 : }
965 :
966 1 : const Real p = pc * std::exp(sum / Tr);
967 1 : const Real dp_dT = -p / Tc / Tr * (sum / Tr + dsum);
968 :
969 : // Henry's constant and its derivative wrt temperature
970 1 : Kh = p * std::exp(lnkh);
971 1 : dKh_dT = (p * dlnkh_dT + dp_dT) * std::exp(lnkh);
972 1 : }
973 :
974 : ADReal
975 0 : Water97FluidProperties::henryConstant(const ADReal & temperature,
976 : const std::vector<Real> & coeffs) const
977 : {
978 0 : const Real T = temperature.value();
979 0 : Real Kh_real = 0.0;
980 0 : Real dKh_dT_real = 0.0;
981 0 : henryConstant(T, coeffs, Kh_real, dKh_dT_real);
982 :
983 0 : ADReal Kh = Kh_real;
984 0 : Kh.derivatives() = temperature.derivatives() * dKh_dT_real;
985 :
986 0 : return Kh;
987 : }
988 :
989 : unsigned int
990 154 : Water97FluidProperties::subregion3(const Real pressure, const Real temperature) const
991 : {
992 154 : Real pMPa = pressure / 1.0e6;
993 : const Real P3cd = 19.00881189173929;
994 : unsigned int subregion = 0;
995 :
996 154 : if (pMPa > 40.0 && pMPa <= 100.0)
997 : {
998 16 : if (temperature <= tempXY(pressure, AB))
999 : subregion = 0;
1000 : else // (temperature > tempXY(pressure, AB))
1001 : subregion = 1;
1002 : }
1003 138 : else if (pMPa > 25.0 && pMPa <= 40.0)
1004 : {
1005 48 : if (temperature <= tempXY(pressure, CD))
1006 : subregion = 2;
1007 12 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, AB))
1008 : subregion = 3;
1009 8 : else if (temperature > tempXY(pressure, AB) && temperature <= tempXY(pressure, EF))
1010 : subregion = 4;
1011 : else // (temperature > tempXY(pressure, EF))
1012 : subregion = 5;
1013 : }
1014 90 : else if (pMPa > 23.5 && pMPa <= 25.0)
1015 : {
1016 16 : if (temperature <= tempXY(pressure, CD))
1017 : subregion = 2;
1018 16 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, GH))
1019 : subregion = 6;
1020 12 : else if (temperature > tempXY(pressure, GH) && temperature <= tempXY(pressure, EF))
1021 : subregion = 7;
1022 8 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, IJ))
1023 : subregion = 8;
1024 4 : else if (temperature > tempXY(pressure, IJ) && temperature <= tempXY(pressure, JK))
1025 : subregion = 9;
1026 : else // (temperature > tempXY(pressure, JK))
1027 : subregion = 10;
1028 : }
1029 74 : else if (pMPa > 23.0 && pMPa <= 23.5)
1030 : {
1031 2 : if (temperature <= tempXY(pressure, CD))
1032 : subregion = 2;
1033 2 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, GH))
1034 : subregion = 11;
1035 2 : else if (temperature > tempXY(pressure, GH) && temperature <= tempXY(pressure, EF))
1036 : subregion = 7;
1037 2 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, IJ))
1038 : subregion = 8;
1039 2 : else if (temperature > tempXY(pressure, IJ) && temperature <= tempXY(pressure, JK))
1040 : subregion = 9;
1041 : else // (temperature > tempXY(pressure, JK))
1042 : subregion = 10;
1043 : }
1044 72 : else if (pMPa > 22.5 && pMPa <= 23.0)
1045 : {
1046 22 : if (temperature <= tempXY(pressure, CD))
1047 : subregion = 2;
1048 22 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, GH))
1049 : subregion = 11;
1050 18 : else if (temperature > tempXY(pressure, GH) && temperature <= tempXY(pressure, MN))
1051 : subregion = 12;
1052 14 : else if (temperature > tempXY(pressure, MN) && temperature <= tempXY(pressure, EF))
1053 : subregion = 13;
1054 10 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, OP))
1055 : subregion = 14;
1056 6 : else if (temperature > tempXY(pressure, OP) && temperature <= tempXY(pressure, IJ))
1057 : subregion = 15;
1058 2 : else if (temperature > tempXY(pressure, IJ) && temperature <= tempXY(pressure, JK))
1059 : subregion = 9;
1060 : else // (temperature > tempXY(pressure, JK))
1061 : subregion = 10;
1062 : }
1063 50 : else if (pMPa > vaporPressure(643.15) * 1.0e-6 &&
1064 : pMPa <= 22.5) // vaporPressure(643.15) = 21.04 MPa
1065 : {
1066 40 : if (temperature <= tempXY(pressure, CD))
1067 : subregion = 2;
1068 40 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, QU))
1069 : subregion = 16;
1070 36 : else if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, RX))
1071 : {
1072 24 : if (pMPa > 22.11 && pMPa <= 22.5)
1073 : {
1074 10 : if (temperature <= tempXY(pressure, UV))
1075 : subregion = 20;
1076 10 : else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1077 : subregion = 21;
1078 6 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1079 : subregion = 22;
1080 : else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1081 : subregion = 23;
1082 : }
1083 14 : else if (pMPa > 22.064 && pMPa <= 22.11)
1084 : {
1085 2 : if (temperature <= tempXY(pressure, UV))
1086 : subregion = 20;
1087 2 : else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1088 : subregion = 24;
1089 2 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1090 : subregion = 25;
1091 : else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1092 : subregion = 23;
1093 : }
1094 12 : else if (temperature <= vaporTemperature(pressure))
1095 : {
1096 8 : if (pMPa > 21.93161551 && pMPa <= 22.064)
1097 6 : if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, UV))
1098 : subregion = 20;
1099 : else
1100 : subregion = 24;
1101 : else // (pMPa > vaporPressure(643.15) * 1.0e-6 && pMPa <= 21.93161551)
1102 : subregion = 20;
1103 : }
1104 4 : else if (temperature > vaporTemperature(pressure))
1105 : {
1106 4 : if (pMPa > 21.90096265 && pMPa <= 22.064)
1107 : {
1108 4 : if (temperature <= tempXY(pressure, WX))
1109 : subregion = 25;
1110 : else
1111 : subregion = 23;
1112 : }
1113 : else
1114 : subregion = 23;
1115 : }
1116 : }
1117 12 : else if (temperature > tempXY(pressure, RX) && temperature <= tempXY(pressure, JK))
1118 : subregion = 17;
1119 : else
1120 : subregion = 10;
1121 : }
1122 10 : else if (pMPa > 20.5 &&
1123 0 : pMPa <= vaporPressure(643.15) * 1.0e-6) // vaporPressure(643.15) = 21.04 MPa
1124 : {
1125 0 : if (temperature <= tempXY(pressure, CD))
1126 : subregion = 2;
1127 0 : else if (temperature > tempXY(pressure, CD) && temperature <= vaporTemperature(pressure))
1128 : subregion = 16;
1129 0 : else if (temperature > vaporTemperature(pressure) && temperature <= tempXY(pressure, JK))
1130 : subregion = 17;
1131 : else // (temperature > tempXY(pressure, JK))
1132 : subregion = 10;
1133 : }
1134 10 : else if (pMPa > P3cd && pMPa <= 20.5) // P3cd = 19.00881189173929
1135 : {
1136 8 : if (temperature <= tempXY(pressure, CD))
1137 : subregion = 2;
1138 6 : else if (temperature > tempXY(pressure, CD) && temperature <= vaporTemperature(pressure))
1139 : subregion = 18;
1140 : else
1141 : subregion = 19;
1142 : }
1143 2 : else if (pMPa > vaporPressure(623.15) * 1.0e-6 && pMPa <= P3cd)
1144 : {
1145 2 : if (temperature < vaporTemperature(pressure))
1146 : subregion = 2;
1147 : else
1148 : subregion = 19;
1149 : }
1150 0 : else if (pMPa > 22.11 && pMPa <= 22.5)
1151 : {
1152 0 : if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, UV))
1153 : subregion = 20;
1154 0 : else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1155 : subregion = 21;
1156 0 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1157 : subregion = 22;
1158 : else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1159 : subregion = 23;
1160 : }
1161 0 : else if (pMPa > 22.064 && pMPa <= 22.11)
1162 : {
1163 0 : if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, UV))
1164 : subregion = 20;
1165 0 : else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1166 : subregion = 24;
1167 0 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1168 : subregion = 25;
1169 : else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1170 : subregion = 23;
1171 : }
1172 : else
1173 0 : mooseError("subregion3(): Shouldn't have got here!");
1174 :
1175 154 : return subregion;
1176 : }
1177 :
1178 : unsigned int
1179 884162 : Water97FluidProperties::inRegion(const Real pressure, const Real temperature) const
1180 : {
1181 : // Valid for 273.15 K <= T <= 1073.15 K, p <= 100 MPa
1182 : // 1073.15 K <= T <= 2273.15 K, p <= 50 Mpa
1183 884162 : if (temperature >= 273.15 && temperature <= 1073.15)
1184 : {
1185 884111 : if (pressure < vaporPressure(273.15) || pressure > 100.0e6)
1186 1 : mooseException("Pressure ", pressure, " is out of range in ", name(), ": inRegion()");
1187 : }
1188 51 : else if (temperature > 1073.15 && temperature <= 2273.15)
1189 : {
1190 51 : if (pressure < 0.0 || pressure > 50.0e6)
1191 1 : mooseException("Pressure ", pressure, " is out of range in ", name(), ": inRegion()");
1192 : }
1193 : else
1194 0 : mooseException("Temperature ", temperature, " is out of range in ", name(), ": inRegion()");
1195 :
1196 : // Determine the phase region that the (P, T) point lies in
1197 : unsigned int region;
1198 :
1199 884160 : if (temperature >= 273.15 && temperature <= 623.15)
1200 : {
1201 883999 : if (pressure > vaporPressure(temperature) && pressure <= 100.0e6)
1202 : region = 1;
1203 : else
1204 : region = 2;
1205 : }
1206 161 : else if (temperature > 623.15 && temperature <= 863.15)
1207 : {
1208 92 : if (pressure <= b23p(temperature))
1209 : region = 2;
1210 : else
1211 : region = 3;
1212 : }
1213 69 : else if (temperature > 863.15 && temperature <= 1073.15)
1214 : region = 2;
1215 : else
1216 : region = 5;
1217 :
1218 884160 : return region;
1219 : }
|
