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 130 : 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 : using std::pow, std::sqrt;
547 : // Check whether the input pressure is within the region of validity of this equation.
548 : // Valid for 611.213 Pa <= p <= 22.064 MPa
549 60 : if (pressure.value() < 611.23 || pressure.value() > _p_critical)
550 0 : mooseException(name() + ": vaporTemperature(): Pressure ",
551 : pressure.value(),
552 : " is outside range 611.213 Pa <= p <= 22.064 MPa");
553 :
554 120 : const ADReal beta = pow(pressure / 1.e6, 0.25);
555 : const ADReal beta2 = beta * beta;
556 60 : const ADReal e = beta2 + _n4[2] * beta + _n4[5];
557 60 : const ADReal f = _n4[0] * beta2 + _n4[3] * beta + _n4[6];
558 60 : const ADReal g = _n4[1] * beta2 + _n4[4] * beta + _n4[7];
559 180 : const ADReal d = 2.0 * g / (-f - sqrt(f * f - 4.0 * e * g));
560 :
561 240 : return (_n4[9] + d - sqrt((_n4[9] + d) * (_n4[9] + d) - 4.0 * (_n4[8] + _n4[9] * d))) / 2.0;
562 : }
563 :
564 : Real
565 60 : Water97FluidProperties::vaporTemperature(Real pressure) const
566 : {
567 60 : const ADReal p = pressure;
568 :
569 60 : return vaporTemperature_ad(p).value();
570 : }
571 :
572 : void
573 0 : Water97FluidProperties::vaporTemperature(Real pressure, Real & Tsat, Real & dTsat_dp) const
574 : {
575 0 : ADReal p = pressure;
576 : Moose::derivInsert(p.derivatives(), 0, 1.0);
577 :
578 0 : const ADReal T = vaporTemperature_ad(p);
579 :
580 0 : Tsat = T.value();
581 0 : dTsat_dp = T.derivatives()[0];
582 0 : }
583 :
584 : Real
585 103 : Water97FluidProperties::b23p(Real temperature) const
586 : {
587 : // Check whether the input temperature is within the region of validity of this equation.
588 : // Valid for 623.15 K <= t <= 863.15 K
589 103 : if (temperature < 623.15 || temperature > 863.15)
590 0 : mooseException(name(),
591 : ": b23p(): Temperature ",
592 : temperature,
593 : " is outside range of 623.15 K <= T <= 863.15 K");
594 :
595 103 : return (_n23[0] + _n23[1] * temperature + _n23[2] * temperature * temperature) * 1.e6;
596 : }
597 :
598 : Real
599 20 : Water97FluidProperties::b23T(Real pressure) const
600 : {
601 : // Check whether the input pressure is within the region of validity of this equation.
602 : // Valid for 16.529 MPa <= p <= 100 MPa
603 20 : if (pressure < 16.529e6 || pressure > 100.0e6)
604 0 : mooseException(
605 : name(), ": b23T(): Pressure ", pressure, "is outside range 16.529 MPa <= p <= 100 MPa");
606 :
607 20 : return _n23[3] + std::sqrt((pressure / 1.e6 - _n23[4]) / _n23[2]);
608 : }
609 :
610 : unsigned int
611 39 : Water97FluidProperties::inRegionPH(Real pressure, Real enthalpy) const
612 : {
613 : unsigned int region;
614 :
615 : // Need to calculate enthalpies at the boundaries to delineate regions
616 39 : Real p273 = vaporPressure(273.15);
617 39 : Real p623 = vaporPressure(623.15);
618 :
619 39 : if (pressure >= p273 && pressure <= p623)
620 : {
621 42 : if (enthalpy >= h_from_p_T(pressure, 273.15) &&
622 21 : enthalpy <= h_from_p_T(pressure, vaporTemperature(pressure)))
623 : region = 1;
624 24 : else if (enthalpy > h_from_p_T(pressure, vaporTemperature(pressure)) &&
625 12 : enthalpy <= h_from_p_T(pressure, 1073.15))
626 : region = 2;
627 0 : else if (enthalpy > h_from_p_T(pressure, 1073.15) && enthalpy <= h_from_p_T(pressure, 2273.15))
628 : region = 5;
629 : else
630 0 : mooseException("Enthalpy ", enthalpy, " is out of range in ", name(), ": inRegionPH()");
631 : }
632 18 : else if (pressure > p623 && pressure <= 50.0e6)
633 : {
634 9 : if (enthalpy >= h_from_p_T(pressure, 273.15) && enthalpy <= h_from_p_T(pressure, 623.15))
635 : region = 1;
636 18 : else if (enthalpy > h_from_p_T(pressure, 623.15) &&
637 9 : enthalpy <= h_from_p_T(pressure, b23T(pressure)))
638 : region = 3;
639 6 : else if (enthalpy > h_from_p_T(pressure, b23T(pressure)) &&
640 3 : enthalpy <= h_from_p_T(pressure, 1073.15))
641 : region = 2;
642 0 : else if (enthalpy > h_from_p_T(pressure, 1073.15) && enthalpy <= h_from_p_T(pressure, 2273.15))
643 : region = 5;
644 : else
645 0 : mooseException("Enthalpy ", enthalpy, " is out of range in ", name(), ": inRegionPH()");
646 : }
647 9 : else if (pressure > 50.0e6 && pressure <= 100.0e6)
648 : {
649 9 : if (enthalpy >= h_from_p_T(pressure, 273.15) && enthalpy <= h_from_p_T(pressure, 623.15))
650 : region = 1;
651 10 : else if (enthalpy > h_from_p_T(pressure, 623.15) &&
652 5 : enthalpy <= h_from_p_T(pressure, b23T(pressure)))
653 : region = 3;
654 4 : else if (enthalpy > h_from_p_T(pressure, b23T(pressure)) &&
655 2 : enthalpy <= h_from_p_T(pressure, 1073.15))
656 : region = 2;
657 : else
658 0 : mooseException("Enthalpy ", enthalpy, " is out of range in ", name(), ": inRegionPH()");
659 : }
660 : else
661 0 : mooseException("Pressure ", pressure, " is out of range in ", name(), ": inRegionPH()");
662 :
663 39 : return region;
664 : }
665 :
666 : unsigned int
667 17 : Water97FluidProperties::subregion2ph(Real pressure, Real enthalpy) const
668 : {
669 : unsigned int subregion;
670 :
671 17 : if (pressure <= 4.0e6)
672 : subregion = 1;
673 7 : else if (pressure > 4.0e6 && pressure < 6.5467e6)
674 : subregion = 2;
675 : else
676 : {
677 5 : if (enthalpy >= b2bc(pressure))
678 : subregion = 2;
679 : else
680 : subregion = 3;
681 : }
682 :
683 17 : return subregion;
684 : }
685 :
686 : unsigned int
687 9 : Water97FluidProperties::subregion3ph(Real pressure, Real enthalpy) const
688 : {
689 : unsigned int subregion;
690 :
691 9 : if (enthalpy <= b3ab(pressure))
692 : subregion = 1;
693 : else
694 : subregion = 2;
695 :
696 9 : return subregion;
697 : }
698 :
699 : Real
700 6 : Water97FluidProperties::b2bc(Real pressure) const
701 : {
702 : // Check whether the input pressure is within the region of validity of this equation.
703 6 : if (pressure < 6.5467e6 || pressure > 100.0e6)
704 0 : mooseException(
705 : name(), ": b2bc(): Pressure ", pressure, " is outside range of 6.5467 MPa <= p <= 100 MPa");
706 :
707 6 : Real pi = pressure / 1.0e6;
708 :
709 6 : return (0.26526571908428e4 + std::sqrt((pi - 0.45257578905948e1) / 0.12809002730136e-3)) * 1.0e3;
710 : }
711 :
712 : Real
713 36 : Water97FluidProperties::T_from_p_h(Real pressure, Real enthalpy) const
714 : {
715 36 : const ADReal p = pressure;
716 36 : const ADReal h = enthalpy;
717 :
718 36 : return T_from_p_h_ad(p, h).value();
719 : }
720 :
721 : void
722 0 : Water97FluidProperties::T_from_p_h(
723 : Real pressure, Real enthalpy, Real & temperature, Real & dT_dp, Real & dT_dh) const
724 : {
725 0 : ADReal p = pressure;
726 : Moose::derivInsert(p.derivatives(), 0, 1.0);
727 0 : ADReal h = enthalpy;
728 : Moose::derivInsert(h.derivatives(), 1, 1.0);
729 :
730 0 : const ADReal T = T_from_p_h_ad(p, h);
731 :
732 0 : temperature = T.value();
733 0 : dT_dp = T.derivatives()[0];
734 0 : dT_dh = T.derivatives()[1];
735 0 : }
736 :
737 : ADReal
738 2 : Water97FluidProperties::T_from_p_h(const ADReal & pressure, const ADReal & enthalpy) const
739 : {
740 2 : return T_from_p_h_ad(pressure, enthalpy);
741 : }
742 :
743 : ADReal
744 39 : Water97FluidProperties::T_from_p_h_ad(const ADReal & pressure, const ADReal & enthalpy) const
745 : {
746 39 : ADReal temperature = 0.0;
747 :
748 : // Determine which region the point is in
749 39 : const unsigned int region = inRegionPH(pressure.value(), enthalpy.value());
750 :
751 39 : switch (region)
752 : {
753 13 : case 1:
754 13 : temperature = temperature_from_ph1(pressure, enthalpy);
755 13 : break;
756 :
757 : case 2:
758 : {
759 : // First, determine which subregion the point is in:
760 17 : const unsigned int subregion = subregion2ph(pressure.value(), enthalpy.value());
761 :
762 17 : if (subregion == 1)
763 10 : temperature = temperature_from_ph2a(pressure, enthalpy);
764 7 : else if (subregion == 2)
765 3 : temperature = temperature_from_ph2b(pressure, enthalpy);
766 : else
767 4 : temperature = temperature_from_ph2c(pressure, enthalpy);
768 : break;
769 : }
770 :
771 : case 3:
772 : {
773 : // First, determine which subregion the point is in:
774 9 : const unsigned int subregion = subregion3ph(pressure.value(), enthalpy.value());
775 :
776 9 : if (subregion == 1)
777 4 : temperature = temperature_from_ph3a(pressure, enthalpy);
778 : else
779 5 : temperature = temperature_from_ph3b(pressure, enthalpy);
780 : break;
781 : }
782 :
783 0 : case 5:
784 0 : mooseError("temperature_from_ph() not implemented for region 5");
785 : break;
786 :
787 0 : default:
788 0 : mooseError("inRegionPH() has given an incorrect region");
789 : }
790 :
791 39 : return temperature;
792 : }
793 :
794 : ADReal
795 13 : Water97FluidProperties::temperature_from_ph1(const ADReal & pressure, const ADReal & enthalpy) const
796 : {
797 : using std::pow;
798 :
799 13 : const ADReal pi = pressure / 1.0e6;
800 13 : const ADReal eta = enthalpy / 2500.0e3;
801 13 : ADReal sum = 0.0;
802 :
803 273 : for (std::size_t i = 0; i < _nph1.size(); ++i)
804 520 : sum += _nph1[i] * pow(pi, _Iph1[i]) * pow(eta + 1.0, _Jph1[i]);
805 :
806 13 : return sum;
807 : }
808 :
809 : ADReal
810 10 : Water97FluidProperties::temperature_from_ph2a(const ADReal & pressure,
811 : const ADReal & enthalpy) const
812 : {
813 : using std::pow, std::abs;
814 :
815 10 : const ADReal pi = pressure / 1.0e6;
816 10 : const ADReal eta = enthalpy / 2000.0e3;
817 10 : ADReal sum = 0.0;
818 :
819 : // Factor out the negative in pow(eta - 2.1, _Jph2a[i]) to avoid fpe in dbg (see #13163)
820 10 : const Real sgn = MathUtils::sign(eta.value() - 2.1);
821 :
822 370 : for (std::size_t i = 0; i < _nph2a.size(); ++i)
823 720 : sum += _nph2a[i] * pow(pi, _Iph2a[i]) * pow(abs(eta - 2.1), _Jph2a[i]) * pow(sgn, _Jph2a[i]);
824 :
825 10 : return sum;
826 : }
827 :
828 : ADReal
829 3 : Water97FluidProperties::temperature_from_ph2b(const ADReal & pressure,
830 : const ADReal & enthalpy) const
831 : {
832 : using std::pow, std::abs;
833 :
834 3 : const ADReal pi = pressure / 1.0e6;
835 3 : const ADReal eta = enthalpy / 2000.0e3;
836 3 : ADReal sum = 0.0;
837 :
838 : // Factor out the negatives in pow(pi - 2.0, _Iph2b[i])* pow(eta - 2.6, _Jph2b[i])
839 : // to avoid fpe in dbg (see #13163)
840 3 : const Real sgn0 = MathUtils::sign(pi.value() - 2.0);
841 3 : const Real sgn1 = MathUtils::sign(eta.value() - 2.6);
842 :
843 117 : for (std::size_t i = 0; i < _nph2b.size(); ++i)
844 114 : sum += _nph2b[i] * pow(abs(pi - 2.0), _Iph2b[i]) * pow(sgn0, _Iph2b[i]) *
845 228 : pow(abs(eta - 2.6), _Jph2b[i]) * pow(sgn1, _Jph2b[i]);
846 :
847 3 : return sum;
848 : }
849 :
850 : ADReal
851 4 : Water97FluidProperties::temperature_from_ph2c(const ADReal & pressure,
852 : const ADReal & enthalpy) const
853 : {
854 : using std::pow, std::abs;
855 :
856 4 : const ADReal pi = pressure / 1.0e6;
857 4 : const ADReal eta = enthalpy / 2000.0e3;
858 4 : ADReal sum = 0.0;
859 :
860 : // Factor out the negative in pow(eta - 1.8, _Jph2c[i]) to avoid fpe in dbg (see #13163)
861 4 : const Real sgn = MathUtils::sign(eta.value() - 1.8);
862 :
863 96 : for (std::size_t i = 0; i < _nph2c.size(); ++i)
864 184 : sum += _nph2c[i] * pow(pi + 25.0, _Iph2c[i]) * pow(abs(eta - 1.8), _Jph2c[i]) *
865 184 : pow(sgn, _Jph2c[i]);
866 :
867 4 : return sum;
868 : }
869 :
870 : Real
871 10 : Water97FluidProperties::b3ab(Real pressure) const
872 : {
873 : // Check whether the input pressure is within the region of validity of this equation.
874 10 : if (pressure < b23p(623.15) || pressure > 100.0e6)
875 0 : mooseException(
876 : name(), ": b3ab(): Pressure ", pressure, "is outside range of 16.529 MPa <= p <= 100 MPa");
877 :
878 10 : Real pi = pressure / 1.0e6;
879 10 : Real eta = 0.201464004206875e4 + 0.374696550136983e1 * pi - 0.219921901054187e-1 * pi * pi +
880 10 : 0.87513168600995e-4 * pi * pi * pi;
881 :
882 10 : return eta * 1.0e3;
883 : }
884 :
885 : ADReal
886 4 : Water97FluidProperties::temperature_from_ph3a(const ADReal & pressure,
887 : const ADReal & enthalpy) const
888 : {
889 : using std::pow;
890 :
891 4 : const ADReal pi = pressure / 100.0e6;
892 4 : const ADReal eta = enthalpy / 2300.0e3;
893 4 : ADReal sum = 0.0;
894 :
895 128 : for (std::size_t i = 0; i < _nph3a.size(); ++i)
896 372 : sum += _nph3a[i] * pow(pi + 0.24, _Iph3a[i]) * pow(eta - 0.615, _Jph3a[i]);
897 :
898 8 : return sum * 760.0;
899 : }
900 :
901 : ADReal
902 5 : Water97FluidProperties::temperature_from_ph3b(const ADReal & pressure,
903 : const ADReal & enthalpy) const
904 : {
905 : using std::pow;
906 :
907 5 : const ADReal pi = pressure / 100.0e6;
908 5 : const ADReal eta = enthalpy / 2800.0e3;
909 5 : ADReal sum = 0.0;
910 :
911 170 : for (size_t i = 0; i < _nph3b.size(); ++i)
912 495 : sum += _nph3b[i] * pow(pi + 0.298, _Iph3b[i]) * pow(eta - 0.72, _Jph3b[i]);
913 :
914 10 : return sum * 860.0;
915 : }
916 :
917 : Real
918 9 : Water97FluidProperties::henryConstant(Real temperature, const std::vector<Real> & coeffs) const
919 : {
920 9 : const Real A = coeffs[0];
921 9 : const Real B = coeffs[1];
922 9 : const Real C = coeffs[2];
923 :
924 9 : const Real Tr = temperature / 647.096;
925 9 : const Real tau = 1.0 - Tr;
926 :
927 : const Real lnkh =
928 9 : A / Tr + B * std::pow(tau, 0.355) / Tr + C * std::pow(Tr, -0.41) * std::exp(tau);
929 :
930 : // The vapor pressure used in this formulation
931 : const std::vector<Real> a{
932 9 : -7.85951783, 1.84408259, -11.7866497, 22.6807411, -15.9618719, 1.80122502};
933 9 : const std::vector<Real> b{1.0, 1.5, 3.0, 3.5, 4.0, 7.5};
934 : Real sum = 0.0;
935 :
936 63 : for (std::size_t i = 0; i < a.size(); ++i)
937 54 : sum += a[i] * std::pow(tau, b[i]);
938 :
939 9 : return 22.064e6 * std::exp(sum / Tr) * std::exp(lnkh);
940 9 : }
941 :
942 : void
943 1 : Water97FluidProperties::henryConstant(Real temperature,
944 : const std::vector<Real> & coeffs,
945 : Real & Kh,
946 : Real & dKh_dT) const
947 : {
948 1 : const Real A = coeffs[0];
949 1 : const Real B = coeffs[1];
950 1 : const Real C = coeffs[2];
951 :
952 : const Real pc = 22.064e6;
953 : const Real Tc = 647.096;
954 :
955 1 : const Real Tr = temperature / Tc;
956 1 : const Real tau = 1.0 - Tr;
957 :
958 : const Real lnkh =
959 1 : A / Tr + B * std::pow(tau, 0.355) / Tr + C * std::pow(Tr, -0.41) * std::exp(tau);
960 : const Real dlnkh_dT =
961 1 : (-A / Tr / Tr - B * std::pow(tau, 0.355) / Tr / Tr - 0.355 * B * std::pow(tau, -0.645) / Tr -
962 1 : 0.41 * C * std::pow(Tr, -1.41) * std::exp(tau) - C * std::pow(Tr, -0.41) * std::exp(tau)) /
963 1 : Tc;
964 :
965 : // The vapor pressure used in this formulation
966 : const std::vector<Real> a{
967 1 : -7.85951783, 1.84408259, -11.7866497, 22.6807411, -15.9618719, 1.80122502};
968 1 : const std::vector<Real> b{1.0, 1.5, 3.0, 3.5, 4.0, 7.5};
969 : Real sum = 0.0;
970 : Real dsum = 0.0;
971 :
972 7 : for (std::size_t i = 0; i < a.size(); ++i)
973 : {
974 6 : sum += a[i] * std::pow(tau, b[i]);
975 6 : dsum += a[i] * b[i] * std::pow(tau, b[i] - 1.0);
976 : }
977 :
978 1 : const Real p = pc * std::exp(sum / Tr);
979 1 : const Real dp_dT = -p / Tc / Tr * (sum / Tr + dsum);
980 :
981 : // Henry's constant and its derivative wrt temperature
982 1 : Kh = p * std::exp(lnkh);
983 1 : dKh_dT = (p * dlnkh_dT + dp_dT) * std::exp(lnkh);
984 1 : }
985 :
986 : ADReal
987 0 : Water97FluidProperties::henryConstant(const ADReal & temperature,
988 : const std::vector<Real> & coeffs) const
989 : {
990 0 : const Real T = temperature.value();
991 0 : Real Kh_real = 0.0;
992 0 : Real dKh_dT_real = 0.0;
993 0 : henryConstant(T, coeffs, Kh_real, dKh_dT_real);
994 :
995 0 : ADReal Kh = Kh_real;
996 0 : Kh.derivatives() = temperature.derivatives() * dKh_dT_real;
997 :
998 0 : return Kh;
999 : }
1000 :
1001 : unsigned int
1002 154 : Water97FluidProperties::subregion3(const Real pressure, const Real temperature) const
1003 : {
1004 154 : Real pMPa = pressure / 1.0e6;
1005 : const Real P3cd = 19.00881189173929;
1006 : unsigned int subregion = 0;
1007 :
1008 154 : if (pMPa > 40.0 && pMPa <= 100.0)
1009 : {
1010 16 : if (temperature <= tempXY(pressure, AB))
1011 : subregion = 0;
1012 : else // (temperature > tempXY(pressure, AB))
1013 : subregion = 1;
1014 : }
1015 138 : else if (pMPa > 25.0 && pMPa <= 40.0)
1016 : {
1017 48 : if (temperature <= tempXY(pressure, CD))
1018 : subregion = 2;
1019 12 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, AB))
1020 : subregion = 3;
1021 8 : else if (temperature > tempXY(pressure, AB) && temperature <= tempXY(pressure, EF))
1022 : subregion = 4;
1023 : else // (temperature > tempXY(pressure, EF))
1024 : subregion = 5;
1025 : }
1026 90 : else if (pMPa > 23.5 && pMPa <= 25.0)
1027 : {
1028 16 : if (temperature <= tempXY(pressure, CD))
1029 : subregion = 2;
1030 16 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, GH))
1031 : subregion = 6;
1032 12 : else if (temperature > tempXY(pressure, GH) && temperature <= tempXY(pressure, EF))
1033 : subregion = 7;
1034 8 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, IJ))
1035 : subregion = 8;
1036 4 : else if (temperature > tempXY(pressure, IJ) && temperature <= tempXY(pressure, JK))
1037 : subregion = 9;
1038 : else // (temperature > tempXY(pressure, JK))
1039 : subregion = 10;
1040 : }
1041 74 : else if (pMPa > 23.0 && pMPa <= 23.5)
1042 : {
1043 2 : if (temperature <= tempXY(pressure, CD))
1044 : subregion = 2;
1045 2 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, GH))
1046 : subregion = 11;
1047 2 : else if (temperature > tempXY(pressure, GH) && temperature <= tempXY(pressure, EF))
1048 : subregion = 7;
1049 2 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, IJ))
1050 : subregion = 8;
1051 2 : else if (temperature > tempXY(pressure, IJ) && temperature <= tempXY(pressure, JK))
1052 : subregion = 9;
1053 : else // (temperature > tempXY(pressure, JK))
1054 : subregion = 10;
1055 : }
1056 72 : else if (pMPa > 22.5 && pMPa <= 23.0)
1057 : {
1058 22 : if (temperature <= tempXY(pressure, CD))
1059 : subregion = 2;
1060 22 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, GH))
1061 : subregion = 11;
1062 18 : else if (temperature > tempXY(pressure, GH) && temperature <= tempXY(pressure, MN))
1063 : subregion = 12;
1064 14 : else if (temperature > tempXY(pressure, MN) && temperature <= tempXY(pressure, EF))
1065 : subregion = 13;
1066 10 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, OP))
1067 : subregion = 14;
1068 6 : else if (temperature > tempXY(pressure, OP) && temperature <= tempXY(pressure, IJ))
1069 : subregion = 15;
1070 2 : else if (temperature > tempXY(pressure, IJ) && temperature <= tempXY(pressure, JK))
1071 : subregion = 9;
1072 : else // (temperature > tempXY(pressure, JK))
1073 : subregion = 10;
1074 : }
1075 50 : else if (pMPa > vaporPressure(643.15) * 1.0e-6 &&
1076 : pMPa <= 22.5) // vaporPressure(643.15) = 21.04 MPa
1077 : {
1078 40 : if (temperature <= tempXY(pressure, CD))
1079 : subregion = 2;
1080 40 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, QU))
1081 : subregion = 16;
1082 36 : else if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, RX))
1083 : {
1084 24 : if (pMPa > 22.11 && pMPa <= 22.5)
1085 : {
1086 10 : if (temperature <= tempXY(pressure, UV))
1087 : subregion = 20;
1088 10 : else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1089 : subregion = 21;
1090 6 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1091 : subregion = 22;
1092 : else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1093 : subregion = 23;
1094 : }
1095 14 : else if (pMPa > 22.064 && pMPa <= 22.11)
1096 : {
1097 2 : if (temperature <= tempXY(pressure, UV))
1098 : subregion = 20;
1099 2 : else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1100 : subregion = 24;
1101 2 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1102 : subregion = 25;
1103 : else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1104 : subregion = 23;
1105 : }
1106 12 : else if (temperature <= vaporTemperature(pressure))
1107 : {
1108 8 : if (pMPa > 21.93161551 && pMPa <= 22.064)
1109 6 : if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, UV))
1110 : subregion = 20;
1111 : else
1112 : subregion = 24;
1113 : else // (pMPa > vaporPressure(643.15) * 1.0e-6 && pMPa <= 21.93161551)
1114 : subregion = 20;
1115 : }
1116 4 : else if (temperature > vaporTemperature(pressure))
1117 : {
1118 4 : if (pMPa > 21.90096265 && pMPa <= 22.064)
1119 : {
1120 4 : if (temperature <= tempXY(pressure, WX))
1121 : subregion = 25;
1122 : else
1123 : subregion = 23;
1124 : }
1125 : else
1126 : subregion = 23;
1127 : }
1128 : }
1129 12 : else if (temperature > tempXY(pressure, RX) && temperature <= tempXY(pressure, JK))
1130 : subregion = 17;
1131 : else
1132 : subregion = 10;
1133 : }
1134 10 : else if (pMPa > 20.5 &&
1135 0 : pMPa <= vaporPressure(643.15) * 1.0e-6) // vaporPressure(643.15) = 21.04 MPa
1136 : {
1137 0 : if (temperature <= tempXY(pressure, CD))
1138 : subregion = 2;
1139 0 : else if (temperature > tempXY(pressure, CD) && temperature <= vaporTemperature(pressure))
1140 : subregion = 16;
1141 0 : else if (temperature > vaporTemperature(pressure) && temperature <= tempXY(pressure, JK))
1142 : subregion = 17;
1143 : else // (temperature > tempXY(pressure, JK))
1144 : subregion = 10;
1145 : }
1146 10 : else if (pMPa > P3cd && pMPa <= 20.5) // P3cd = 19.00881189173929
1147 : {
1148 8 : if (temperature <= tempXY(pressure, CD))
1149 : subregion = 2;
1150 6 : else if (temperature > tempXY(pressure, CD) && temperature <= vaporTemperature(pressure))
1151 : subregion = 18;
1152 : else
1153 : subregion = 19;
1154 : }
1155 2 : else if (pMPa > vaporPressure(623.15) * 1.0e-6 && pMPa <= P3cd)
1156 : {
1157 2 : if (temperature < vaporTemperature(pressure))
1158 : subregion = 2;
1159 : else
1160 : subregion = 19;
1161 : }
1162 0 : else if (pMPa > 22.11 && pMPa <= 22.5)
1163 : {
1164 0 : if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, UV))
1165 : subregion = 20;
1166 0 : else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1167 : subregion = 21;
1168 0 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1169 : subregion = 22;
1170 : else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1171 : subregion = 23;
1172 : }
1173 0 : else if (pMPa > 22.064 && pMPa <= 22.11)
1174 : {
1175 0 : if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, UV))
1176 : subregion = 20;
1177 0 : else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1178 : subregion = 24;
1179 0 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1180 : subregion = 25;
1181 : else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1182 : subregion = 23;
1183 : }
1184 : else
1185 0 : mooseError("subregion3(): Shouldn't have got here!");
1186 :
1187 154 : return subregion;
1188 : }
1189 :
1190 : unsigned int
1191 884162 : Water97FluidProperties::inRegion(const Real pressure, const Real temperature) const
1192 : {
1193 : // Valid for 273.15 K <= T <= 1073.15 K, p <= 100 MPa
1194 : // 1073.15 K <= T <= 2273.15 K, p <= 50 Mpa
1195 884162 : if (temperature >= 273.15 && temperature <= 1073.15)
1196 : {
1197 884111 : if (pressure < vaporPressure(273.15) || pressure > 100.0e6)
1198 1 : mooseException("Pressure ", pressure, " is out of range in ", name(), ": inRegion()");
1199 : }
1200 51 : else if (temperature > 1073.15 && temperature <= 2273.15)
1201 : {
1202 51 : if (pressure < 0.0 || pressure > 50.0e6)
1203 1 : mooseException("Pressure ", pressure, " is out of range in ", name(), ": inRegion()");
1204 : }
1205 : else
1206 0 : mooseException("Temperature ", temperature, " is out of range in ", name(), ": inRegion()");
1207 :
1208 : // Determine the phase region that the (P, T) point lies in
1209 : unsigned int region;
1210 :
1211 884160 : if (temperature >= 273.15 && temperature <= 623.15)
1212 : {
1213 883999 : if (pressure > vaporPressure(temperature) && pressure <= 100.0e6)
1214 : region = 1;
1215 : else
1216 : region = 2;
1217 : }
1218 161 : else if (temperature > 623.15 && temperature <= 863.15)
1219 : {
1220 92 : if (pressure <= b23p(temperature))
1221 : region = 2;
1222 : else
1223 : region = 3;
1224 : }
1225 69 : else if (temperature > 863.15 && temperature <= 1073.15)
1226 : region = 2;
1227 : else
1228 : region = 5;
1229 :
1230 884160 : return region;
1231 : }
|
