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 166 : Water97FluidProperties::validParams()
18 : {
19 166 : InputParameters params = SinglePhaseFluidProperties::validParams();
20 166 : params.addClassDescription("Fluid properties for water and steam (H2O) using IAPWS-IF97");
21 166 : return params;
22 0 : }
23 :
24 84 : Water97FluidProperties::Water97FluidProperties(const InputParameters & parameters)
25 : : SinglePhaseFluidProperties(parameters),
26 84 : _Mh2o(18.015e-3),
27 84 : _Rw(461.526),
28 84 : _p_critical(22.064e6),
29 84 : _T_critical(647.096),
30 84 : _rho_critical(322.0),
31 84 : _p_triple(611.657),
32 8904 : _T_triple(273.16)
33 : {
34 504 : }
35 :
36 84 : Water97FluidProperties::~Water97FluidProperties() {}
37 :
38 : std::string
39 12 : Water97FluidProperties::fluidName() const
40 : {
41 12 : return "water";
42 : }
43 :
44 : Real
45 29 : Water97FluidProperties::molarMass() const
46 : {
47 29 : 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 100654 : Water97FluidProperties::rho_from_p_T(Real pressure, Real temperature) const
82 : {
83 100654 : 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 100186 : Water97FluidProperties::e_from_p_T(Real pressure, Real temperature) const
178 : {
179 100186 : 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 159 : Water97FluidProperties::c_from_p_T(const Real p, const Real T) const
233 : {
234 159 : 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 450 : Water97FluidProperties::cp_from_p_T(Real pressure, Real temperature) const
245 : {
246 450 : 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 151 : Water97FluidProperties::cv_from_p_T(Real pressure, Real temperature) const
264 : {
265 151 : 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 100172 : Water97FluidProperties::mu_from_p_T(Real pressure, Real temperature) const
283 : {
284 100172 : 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 157 : Water97FluidProperties::k_from_p_T(Real pressure, Real temperature) const
396 : {
397 157 : 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 191 : Water97FluidProperties::s_from_p_T(Real pressure, Real temperature) const
432 : {
433 191 : 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 100649 : Water97FluidProperties::h_from_p_T(Real pressure, Real temperature) const
493 : {
494 100649 : 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 806117 : 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 806117 : 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 806117 : theta = temperature + _n4[8] / (temperature - _n4[9]);
534 806117 : theta2 = theta * theta;
535 806117 : a = theta2 + _n4[0] * theta + _n4[1];
536 806117 : b = _n4[2] * theta2 + _n4[3] * theta + _n4[4];
537 806117 : c = _n4[5] * theta2 + _n4[6] * theta + _n4[7];
538 806117 : p = Utility::pow<4>(2.0 * c / (-b + std::sqrt(b * b - 4.0 * a * c)));
539 :
540 806117 : return p * 1.e6;
541 : }
542 :
543 : ADReal
544 72 : 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 72 : 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 144 : const ADReal beta = pow(pressure / 1.e6, 0.25);
555 : const ADReal beta2 = beta * beta;
556 72 : const ADReal e = beta2 + _n4[2] * beta + _n4[5];
557 72 : const ADReal f = _n4[0] * beta2 + _n4[3] * beta + _n4[6];
558 72 : const ADReal g = _n4[1] * beta2 + _n4[4] * beta + _n4[7];
559 216 : const ADReal d = 2.0 * g / (-f - sqrt(f * f - 4.0 * e * g));
560 :
561 288 : 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 72 : Water97FluidProperties::vaporTemperature(Real pressure) const
566 : {
567 72 : const ADReal p = pressure;
568 :
569 72 : 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 110 : 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 110 : 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 110 : return (_n23[0] + _n23[1] * temperature + _n23[2] * temperature * temperature) * 1.e6;
596 : }
597 :
598 : Real
599 24 : 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 24 : 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 24 : return _n23[3] + std::sqrt((pressure / 1.e6 - _n23[4]) / _n23[2]);
608 : }
609 :
610 : unsigned int
611 57 : 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 57 : Real p273 = vaporPressure(273.15);
617 57 : Real p623 = vaporPressure(623.15);
618 :
619 57 : if (enthalpy < h_from_p_T(pressure, 273.15))
620 0 : mooseException("Enthalpy ", enthalpy, " is out of range in ", name(), ": inRegionPH()");
621 :
622 56 : if (pressure >= p273 && pressure <= p623)
623 : {
624 : // Use region 1 definition of h_from_p_T to get the lower bound
625 28 : if (enthalpy <=
626 56 : _Rw * _T_star[0] *
627 28 : dgamma1_dtau(pressure / _p_star[0], _T_star[0] / vaporTemperature(pressure)))
628 : region = 1;
629 : // Use region 2 definition of h_from_p_T to get the upper bound
630 17 : else if (enthalpy <=
631 34 : _Rw * _T_star[1] *
632 17 : dgamma2_dtau(pressure / _p_star[1], _T_star[1] / vaporTemperature(pressure)))
633 : region = 4;
634 15 : else if (enthalpy <= h_from_p_T(pressure, 1073.15))
635 : region = 2;
636 2 : else if (enthalpy <= h_from_p_T(pressure, 2273.15))
637 : region = 5;
638 : else
639 1 : mooseException("Enthalpy ", enthalpy, " is out of range in ", name(), ": inRegionPH()");
640 : }
641 28 : else if (pressure > p623 && pressure <= 50.0e6)
642 : {
643 15 : if (enthalpy <= h_from_p_T(pressure, 623.15))
644 : region = 1;
645 15 : else if (enthalpy <= h_from_p_T(pressure, b23T(pressure)))
646 : region = 3;
647 7 : else if (enthalpy <= h_from_p_T(pressure, 1073.15))
648 : region = 2;
649 3 : else if (enthalpy <= h_from_p_T(pressure, 2273.15))
650 : region = 5;
651 : else
652 1 : mooseException("Enthalpy ", enthalpy, " is out of range in ", name(), ": inRegionPH()");
653 : }
654 13 : else if (pressure > 50.0e6 && pressure <= 100.0e6)
655 : {
656 13 : if (enthalpy <= h_from_p_T(pressure, 623.15))
657 : region = 1;
658 8 : else if (enthalpy <= h_from_p_T(pressure, b23T(pressure)))
659 : region = 3;
660 4 : else if (enthalpy <= h_from_p_T(pressure, 1073.15))
661 : region = 2;
662 : else
663 1 : mooseException("Enthalpy ", enthalpy, " is out of range in ", name(), ": inRegionPH()");
664 : }
665 : else
666 0 : mooseException("Pressure ", pressure, " is out of range in ", name(), ": inRegionPH()");
667 :
668 53 : return region;
669 : }
670 :
671 : unsigned int
672 17 : Water97FluidProperties::subregion2ph(Real pressure, Real enthalpy) const
673 : {
674 : unsigned int subregion;
675 :
676 17 : if (pressure <= 4.0e6)
677 : subregion = 1;
678 7 : else if (pressure > 4.0e6 && pressure < 6.5467e6)
679 : subregion = 2;
680 : else
681 : {
682 5 : if (enthalpy >= b2bc(pressure))
683 : subregion = 2;
684 : else
685 : subregion = 3;
686 : }
687 :
688 17 : return subregion;
689 : }
690 :
691 : unsigned int
692 9 : Water97FluidProperties::subregion3ph(Real pressure, Real enthalpy) const
693 : {
694 : unsigned int subregion;
695 :
696 9 : if (enthalpy <= b3ab(pressure))
697 : subregion = 1;
698 : else
699 : subregion = 2;
700 :
701 9 : return subregion;
702 : }
703 :
704 : Real
705 6 : Water97FluidProperties::b2bc(Real pressure) const
706 : {
707 : // Check whether the input pressure is within the region of validity of this equation.
708 6 : if (pressure < 6.5467e6 || pressure > 100.0e6)
709 0 : mooseException(
710 : name(), ": b2bc(): Pressure ", pressure, " is outside range of 6.5467 MPa <= p <= 100 MPa");
711 :
712 6 : Real pi = pressure / 1.0e6;
713 :
714 6 : return (0.26526571908428e4 + std::sqrt((pi - 0.45257578905948e1) / 0.12809002730136e-3)) * 1.0e3;
715 : }
716 :
717 : Real
718 36 : Water97FluidProperties::T_from_p_h(Real pressure, Real enthalpy) const
719 : {
720 36 : const ADReal p = pressure;
721 36 : const ADReal h = enthalpy;
722 :
723 36 : return T_from_p_h_ad(p, h).value();
724 : }
725 :
726 : void
727 0 : Water97FluidProperties::T_from_p_h(
728 : Real pressure, Real enthalpy, Real & temperature, Real & dT_dp, Real & dT_dh) const
729 : {
730 0 : ADReal p = pressure;
731 : Moose::derivInsert(p.derivatives(), 0, 1.0);
732 0 : ADReal h = enthalpy;
733 : Moose::derivInsert(h.derivatives(), 1, 1.0);
734 :
735 0 : const ADReal T = T_from_p_h_ad(p, h);
736 :
737 0 : temperature = T.value();
738 0 : dT_dp = T.derivatives()[0];
739 0 : dT_dh = T.derivatives()[1];
740 0 : }
741 :
742 : ADReal
743 2 : Water97FluidProperties::T_from_p_h(const ADReal & pressure, const ADReal & enthalpy) const
744 : {
745 2 : return T_from_p_h_ad(pressure, enthalpy);
746 : }
747 :
748 : ADReal
749 39 : Water97FluidProperties::T_from_p_h_ad(const ADReal & pressure, const ADReal & enthalpy) const
750 : {
751 39 : ADReal temperature = 0.0;
752 :
753 : // Determine which region the point is in
754 39 : const unsigned int region = inRegionPH(pressure.value(), enthalpy.value());
755 :
756 39 : switch (region)
757 : {
758 13 : case 1:
759 13 : temperature = temperature_from_ph1(pressure, enthalpy);
760 13 : break;
761 :
762 : case 2:
763 : {
764 : // First, determine which subregion the point is in:
765 17 : const unsigned int subregion = subregion2ph(pressure.value(), enthalpy.value());
766 :
767 17 : if (subregion == 1)
768 10 : temperature = temperature_from_ph2a(pressure, enthalpy);
769 7 : else if (subregion == 2)
770 3 : temperature = temperature_from_ph2b(pressure, enthalpy);
771 : else
772 4 : temperature = temperature_from_ph2c(pressure, enthalpy);
773 : break;
774 : }
775 :
776 : case 3:
777 : {
778 : // First, determine which subregion the point is in:
779 9 : const unsigned int subregion = subregion3ph(pressure.value(), enthalpy.value());
780 :
781 9 : if (subregion == 1)
782 4 : temperature = temperature_from_ph3a(pressure, enthalpy);
783 : else
784 5 : temperature = temperature_from_ph3b(pressure, enthalpy);
785 : break;
786 : }
787 :
788 0 : case 4:
789 0 : return vaporTemperature_ad(pressure);
790 : break;
791 :
792 0 : case 5:
793 : {
794 : Real T_min = 1073.15;
795 : Real T_max = 2273.15;
796 0 : Real T_find = 0.;
797 : Real T_error = 1.0;
798 :
799 : // Bisection solve
800 0 : while (T_error > libMesh::TOLERANCE)
801 : {
802 0 : T_find = 0.5 * (T_min + T_max);
803 0 : Real h_find = h_from_p_T(pressure.value(), T_find);
804 :
805 0 : if (h_find > enthalpy.value())
806 : T_max = T_find;
807 : else
808 : T_min = T_find;
809 :
810 0 : T_error = std::abs((T_max - T_min) / T_find);
811 : }
812 0 : if (!_allow_imperfect_jacobians)
813 0 : mooseDoOnce(
814 : mooseWarning("The derivatives for T_from_p_h_ad are currently neglected, set to 0."));
815 0 : temperature = T_find;
816 :
817 : break;
818 : }
819 0 : default:
820 0 : mooseError("inRegionPH() has given an incorrect region");
821 : }
822 :
823 39 : return temperature;
824 : }
825 :
826 : ADReal
827 13 : Water97FluidProperties::temperature_from_ph1(const ADReal & pressure, const ADReal & enthalpy) const
828 : {
829 : using std::pow;
830 :
831 13 : const ADReal pi = pressure / 1.0e6;
832 13 : const ADReal eta = enthalpy / 2500.0e3;
833 13 : ADReal sum = 0.0;
834 :
835 273 : for (std::size_t i = 0; i < _nph1.size(); ++i)
836 520 : sum += _nph1[i] * pow(pi, _Iph1[i]) * pow(eta + 1.0, _Jph1[i]);
837 :
838 13 : return sum;
839 : }
840 :
841 : ADReal
842 10 : Water97FluidProperties::temperature_from_ph2a(const ADReal & pressure,
843 : const ADReal & enthalpy) const
844 : {
845 : using std::pow, std::abs;
846 :
847 10 : const ADReal pi = pressure / 1.0e6;
848 10 : const ADReal eta = enthalpy / 2000.0e3;
849 10 : ADReal sum = 0.0;
850 :
851 : // Factor out the negative in pow(eta - 2.1, _Jph2a[i]) to avoid fpe in dbg (see #13163)
852 10 : const Real sgn = MathUtils::sign(eta.value() - 2.1);
853 :
854 370 : for (std::size_t i = 0; i < _nph2a.size(); ++i)
855 720 : sum += _nph2a[i] * pow(pi, _Iph2a[i]) * pow(abs(eta - 2.1), _Jph2a[i]) * pow(sgn, _Jph2a[i]);
856 :
857 10 : return sum;
858 : }
859 :
860 : ADReal
861 3 : Water97FluidProperties::temperature_from_ph2b(const ADReal & pressure,
862 : const ADReal & enthalpy) const
863 : {
864 : using std::pow, std::abs;
865 :
866 3 : const ADReal pi = pressure / 1.0e6;
867 3 : const ADReal eta = enthalpy / 2000.0e3;
868 3 : ADReal sum = 0.0;
869 :
870 : // Factor out the negatives in pow(pi - 2.0, _Iph2b[i])* pow(eta - 2.6, _Jph2b[i])
871 : // to avoid fpe in dbg (see #13163)
872 3 : const Real sgn0 = MathUtils::sign(pi.value() - 2.0);
873 3 : const Real sgn1 = MathUtils::sign(eta.value() - 2.6);
874 :
875 117 : for (std::size_t i = 0; i < _nph2b.size(); ++i)
876 114 : sum += _nph2b[i] * pow(abs(pi - 2.0), _Iph2b[i]) * pow(sgn0, _Iph2b[i]) *
877 228 : pow(abs(eta - 2.6), _Jph2b[i]) * pow(sgn1, _Jph2b[i]);
878 :
879 3 : return sum;
880 : }
881 :
882 : ADReal
883 4 : Water97FluidProperties::temperature_from_ph2c(const ADReal & pressure,
884 : const ADReal & enthalpy) const
885 : {
886 : using std::pow, std::abs;
887 :
888 4 : const ADReal pi = pressure / 1.0e6;
889 4 : const ADReal eta = enthalpy / 2000.0e3;
890 4 : ADReal sum = 0.0;
891 :
892 : // Factor out the negative in pow(eta - 1.8, _Jph2c[i]) to avoid fpe in dbg (see #13163)
893 4 : const Real sgn = MathUtils::sign(eta.value() - 1.8);
894 :
895 96 : for (std::size_t i = 0; i < _nph2c.size(); ++i)
896 184 : sum += _nph2c[i] * pow(pi + 25.0, _Iph2c[i]) * pow(abs(eta - 1.8), _Jph2c[i]) *
897 184 : pow(sgn, _Jph2c[i]);
898 :
899 4 : return sum;
900 : }
901 :
902 : Real
903 10 : Water97FluidProperties::b3ab(Real pressure) const
904 : {
905 : // Check whether the input pressure is within the region of validity of this equation.
906 10 : if (pressure < b23p(623.15) || pressure > 100.0e6)
907 0 : mooseException(
908 : name(), ": b3ab(): Pressure ", pressure, "is outside range of 16.529 MPa <= p <= 100 MPa");
909 :
910 10 : Real pi = pressure / 1.0e6;
911 10 : Real eta = 0.201464004206875e4 + 0.374696550136983e1 * pi - 0.219921901054187e-1 * pi * pi +
912 10 : 0.87513168600995e-4 * pi * pi * pi;
913 :
914 10 : return eta * 1.0e3;
915 : }
916 :
917 : ADReal
918 4 : Water97FluidProperties::temperature_from_ph3a(const ADReal & pressure,
919 : const ADReal & enthalpy) const
920 : {
921 : using std::pow;
922 :
923 4 : const ADReal pi = pressure / 100.0e6;
924 4 : const ADReal eta = enthalpy / 2300.0e3;
925 4 : ADReal sum = 0.0;
926 :
927 128 : for (std::size_t i = 0; i < _nph3a.size(); ++i)
928 372 : sum += _nph3a[i] * pow(pi + 0.24, _Iph3a[i]) * pow(eta - 0.615, _Jph3a[i]);
929 :
930 8 : return sum * 760.0;
931 : }
932 :
933 : ADReal
934 5 : Water97FluidProperties::temperature_from_ph3b(const ADReal & pressure,
935 : const ADReal & enthalpy) const
936 : {
937 : using std::pow;
938 :
939 5 : const ADReal pi = pressure / 100.0e6;
940 5 : const ADReal eta = enthalpy / 2800.0e3;
941 5 : ADReal sum = 0.0;
942 :
943 170 : for (size_t i = 0; i < _nph3b.size(); ++i)
944 495 : sum += _nph3b[i] * pow(pi + 0.298, _Iph3b[i]) * pow(eta - 0.72, _Jph3b[i]);
945 :
946 10 : return sum * 860.0;
947 : }
948 :
949 : Real
950 9 : Water97FluidProperties::henryConstant(Real temperature, const std::vector<Real> & coeffs) const
951 : {
952 9 : const Real A = coeffs[0];
953 9 : const Real B = coeffs[1];
954 9 : const Real C = coeffs[2];
955 :
956 9 : const Real Tr = temperature / 647.096;
957 9 : const Real tau = 1.0 - Tr;
958 :
959 : const Real lnkh =
960 9 : A / Tr + B * std::pow(tau, 0.355) / Tr + C * std::pow(Tr, -0.41) * std::exp(tau);
961 :
962 : // The vapor pressure used in this formulation
963 : const std::vector<Real> a{
964 9 : -7.85951783, 1.84408259, -11.7866497, 22.6807411, -15.9618719, 1.80122502};
965 9 : const std::vector<Real> b{1.0, 1.5, 3.0, 3.5, 4.0, 7.5};
966 : Real sum = 0.0;
967 :
968 63 : for (std::size_t i = 0; i < a.size(); ++i)
969 54 : sum += a[i] * std::pow(tau, b[i]);
970 :
971 9 : return 22.064e6 * std::exp(sum / Tr) * std::exp(lnkh);
972 9 : }
973 :
974 : void
975 1 : Water97FluidProperties::henryConstant(Real temperature,
976 : const std::vector<Real> & coeffs,
977 : Real & Kh,
978 : Real & dKh_dT) const
979 : {
980 1 : const Real A = coeffs[0];
981 1 : const Real B = coeffs[1];
982 1 : const Real C = coeffs[2];
983 :
984 : const Real pc = 22.064e6;
985 : const Real Tc = 647.096;
986 :
987 1 : const Real Tr = temperature / Tc;
988 1 : const Real tau = 1.0 - Tr;
989 :
990 : const Real lnkh =
991 1 : A / Tr + B * std::pow(tau, 0.355) / Tr + C * std::pow(Tr, -0.41) * std::exp(tau);
992 : const Real dlnkh_dT =
993 1 : (-A / Tr / Tr - B * std::pow(tau, 0.355) / Tr / Tr - 0.355 * B * std::pow(tau, -0.645) / Tr -
994 1 : 0.41 * C * std::pow(Tr, -1.41) * std::exp(tau) - C * std::pow(Tr, -0.41) * std::exp(tau)) /
995 1 : Tc;
996 :
997 : // The vapor pressure used in this formulation
998 : const std::vector<Real> a{
999 1 : -7.85951783, 1.84408259, -11.7866497, 22.6807411, -15.9618719, 1.80122502};
1000 1 : const std::vector<Real> b{1.0, 1.5, 3.0, 3.5, 4.0, 7.5};
1001 : Real sum = 0.0;
1002 : Real dsum = 0.0;
1003 :
1004 7 : for (std::size_t i = 0; i < a.size(); ++i)
1005 : {
1006 6 : sum += a[i] * std::pow(tau, b[i]);
1007 6 : dsum += a[i] * b[i] * std::pow(tau, b[i] - 1.0);
1008 : }
1009 :
1010 1 : const Real p = pc * std::exp(sum / Tr);
1011 1 : const Real dp_dT = -p / Tc / Tr * (sum / Tr + dsum);
1012 :
1013 : // Henry's constant and its derivative wrt temperature
1014 1 : Kh = p * std::exp(lnkh);
1015 1 : dKh_dT = (p * dlnkh_dT + dp_dT) * std::exp(lnkh);
1016 1 : }
1017 :
1018 : ADReal
1019 0 : Water97FluidProperties::henryConstant(const ADReal & temperature,
1020 : const std::vector<Real> & coeffs) const
1021 : {
1022 0 : const Real T = temperature.value();
1023 0 : Real Kh_real = 0.0;
1024 0 : Real dKh_dT_real = 0.0;
1025 0 : henryConstant(T, coeffs, Kh_real, dKh_dT_real);
1026 :
1027 0 : ADReal Kh = Kh_real;
1028 0 : Kh.derivatives() = temperature.derivatives() * dKh_dT_real;
1029 :
1030 0 : return Kh;
1031 : }
1032 :
1033 : unsigned int
1034 157 : Water97FluidProperties::subregion3(const Real pressure, const Real temperature) const
1035 : {
1036 157 : Real pMPa = pressure / 1.0e6;
1037 : const Real P3cd = 19.00881189173929;
1038 : unsigned int subregion = 0;
1039 :
1040 157 : if (pMPa > 40.0 && pMPa <= 100.0)
1041 : {
1042 17 : if (temperature <= tempXY(pressure, AB))
1043 : subregion = 0;
1044 : else // (temperature > tempXY(pressure, AB))
1045 : subregion = 1;
1046 : }
1047 140 : else if (pMPa > 25.0 && pMPa <= 40.0)
1048 : {
1049 49 : if (temperature <= tempXY(pressure, CD))
1050 : subregion = 2;
1051 12 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, AB))
1052 : subregion = 3;
1053 8 : else if (temperature > tempXY(pressure, AB) && temperature <= tempXY(pressure, EF))
1054 : subregion = 4;
1055 : else // (temperature > tempXY(pressure, EF))
1056 : subregion = 5;
1057 : }
1058 91 : else if (pMPa > 23.5 && pMPa <= 25.0)
1059 : {
1060 16 : if (temperature <= tempXY(pressure, CD))
1061 : subregion = 2;
1062 16 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, GH))
1063 : subregion = 6;
1064 12 : else if (temperature > tempXY(pressure, GH) && temperature <= tempXY(pressure, EF))
1065 : subregion = 7;
1066 8 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, IJ))
1067 : subregion = 8;
1068 4 : else if (temperature > tempXY(pressure, IJ) && temperature <= tempXY(pressure, JK))
1069 : subregion = 9;
1070 : else // (temperature > tempXY(pressure, JK))
1071 : subregion = 10;
1072 : }
1073 75 : else if (pMPa > 23.0 && pMPa <= 23.5)
1074 : {
1075 2 : if (temperature <= tempXY(pressure, CD))
1076 : subregion = 2;
1077 2 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, GH))
1078 : subregion = 11;
1079 2 : else if (temperature > tempXY(pressure, GH) && temperature <= tempXY(pressure, EF))
1080 : subregion = 7;
1081 2 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, IJ))
1082 : subregion = 8;
1083 2 : else if (temperature > tempXY(pressure, IJ) && temperature <= tempXY(pressure, JK))
1084 : subregion = 9;
1085 : else // (temperature > tempXY(pressure, JK))
1086 : subregion = 10;
1087 : }
1088 73 : else if (pMPa > 22.5 && pMPa <= 23.0)
1089 : {
1090 22 : if (temperature <= tempXY(pressure, CD))
1091 : subregion = 2;
1092 22 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, GH))
1093 : subregion = 11;
1094 18 : else if (temperature > tempXY(pressure, GH) && temperature <= tempXY(pressure, MN))
1095 : subregion = 12;
1096 14 : else if (temperature > tempXY(pressure, MN) && temperature <= tempXY(pressure, EF))
1097 : subregion = 13;
1098 10 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, OP))
1099 : subregion = 14;
1100 6 : else if (temperature > tempXY(pressure, OP) && temperature <= tempXY(pressure, IJ))
1101 : subregion = 15;
1102 2 : else if (temperature > tempXY(pressure, IJ) && temperature <= tempXY(pressure, JK))
1103 : subregion = 9;
1104 : else // (temperature > tempXY(pressure, JK))
1105 : subregion = 10;
1106 : }
1107 51 : else if (pMPa > vaporPressure(643.15) * 1.0e-6 &&
1108 : pMPa <= 22.5) // vaporPressure(643.15) = 21.04 MPa
1109 : {
1110 41 : if (temperature <= tempXY(pressure, CD))
1111 : subregion = 2;
1112 41 : else if (temperature > tempXY(pressure, CD) && temperature <= tempXY(pressure, QU))
1113 : subregion = 16;
1114 37 : else if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, RX))
1115 : {
1116 24 : if (pMPa > 22.11 && pMPa <= 22.5)
1117 : {
1118 10 : if (temperature <= tempXY(pressure, UV))
1119 : subregion = 20;
1120 10 : else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1121 : subregion = 21;
1122 6 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1123 : subregion = 22;
1124 : else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1125 : subregion = 23;
1126 : }
1127 14 : else if (pMPa > 22.064 && pMPa <= 22.11)
1128 : {
1129 2 : if (temperature <= tempXY(pressure, UV))
1130 : subregion = 20;
1131 2 : else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1132 : subregion = 24;
1133 2 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1134 : subregion = 25;
1135 : else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1136 : subregion = 23;
1137 : }
1138 12 : else if (temperature <= vaporTemperature(pressure))
1139 : {
1140 8 : if (pMPa > 21.93161551 && pMPa <= 22.064)
1141 6 : if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, UV))
1142 : subregion = 20;
1143 : else
1144 : subregion = 24;
1145 : else // (pMPa > vaporPressure(643.15) * 1.0e-6 && pMPa <= 21.93161551)
1146 : subregion = 20;
1147 : }
1148 4 : else if (temperature > vaporTemperature(pressure))
1149 : {
1150 4 : if (pMPa > 21.90096265 && pMPa <= 22.064)
1151 : {
1152 4 : if (temperature <= tempXY(pressure, WX))
1153 : subregion = 25;
1154 : else
1155 : subregion = 23;
1156 : }
1157 : else
1158 : subregion = 23;
1159 : }
1160 : }
1161 13 : else if (temperature > tempXY(pressure, RX) && temperature <= tempXY(pressure, JK))
1162 : subregion = 17;
1163 : else
1164 : subregion = 10;
1165 : }
1166 10 : else if (pMPa > 20.5 &&
1167 0 : pMPa <= vaporPressure(643.15) * 1.0e-6) // vaporPressure(643.15) = 21.04 MPa
1168 : {
1169 0 : if (temperature <= tempXY(pressure, CD))
1170 : subregion = 2;
1171 0 : else if (temperature > tempXY(pressure, CD) && temperature <= vaporTemperature(pressure))
1172 : subregion = 16;
1173 0 : else if (temperature > vaporTemperature(pressure) && temperature <= tempXY(pressure, JK))
1174 : subregion = 17;
1175 : else // (temperature > tempXY(pressure, JK))
1176 : subregion = 10;
1177 : }
1178 10 : else if (pMPa > P3cd && pMPa <= 20.5) // P3cd = 19.00881189173929
1179 : {
1180 8 : if (temperature <= tempXY(pressure, CD))
1181 : subregion = 2;
1182 6 : else if (temperature > tempXY(pressure, CD) && temperature <= vaporTemperature(pressure))
1183 : subregion = 18;
1184 : else
1185 : subregion = 19;
1186 : }
1187 2 : else if (pMPa > vaporPressure(623.15) * 1.0e-6 && pMPa <= P3cd)
1188 : {
1189 2 : if (temperature < vaporTemperature(pressure))
1190 : subregion = 2;
1191 : else
1192 : subregion = 19;
1193 : }
1194 0 : else if (pMPa > 22.11 && pMPa <= 22.5)
1195 : {
1196 0 : if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, UV))
1197 : subregion = 20;
1198 0 : else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1199 : subregion = 21;
1200 0 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1201 : subregion = 22;
1202 : else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1203 : subregion = 23;
1204 : }
1205 0 : else if (pMPa > 22.064 && pMPa <= 22.11)
1206 : {
1207 0 : if (temperature > tempXY(pressure, QU) && temperature <= tempXY(pressure, UV))
1208 : subregion = 20;
1209 0 : else if (temperature > tempXY(pressure, UV) && temperature <= tempXY(pressure, EF))
1210 : subregion = 24;
1211 0 : else if (temperature > tempXY(pressure, EF) && temperature <= tempXY(pressure, WX))
1212 : subregion = 25;
1213 : else // (temperature > tempXY(pressure, WX) && temperature <= tempXY(pressure, RX))
1214 : subregion = 23;
1215 : }
1216 : else
1217 0 : mooseError("subregion3(): Shouldn't have got here!");
1218 :
1219 157 : return subregion;
1220 : }
1221 :
1222 : unsigned int
1223 403095 : Water97FluidProperties::inRegion(const Real pressure, const Real temperature) const
1224 : {
1225 : // Valid for 273.15 K <= T <= 1073.15 K, p <= 100 MPa
1226 : // 1073.15 K <= T <= 2273.15 K, p <= 50 Mpa
1227 403095 : if (temperature >= 273.15 && temperature <= 1073.15)
1228 : {
1229 403036 : if (pressure < vaporPressure(273.15) || pressure > 100.0e6)
1230 2 : mooseException("Pressure ", pressure, " is out of range in ", name(), ": inRegion()");
1231 : }
1232 59 : else if (temperature > 1073.15 && temperature <= 2273.15)
1233 : {
1234 59 : if (pressure < 0.0 || pressure > 50.0e6)
1235 1 : mooseException("Pressure ", pressure, " is out of range in ", name(), ": inRegion()");
1236 : }
1237 : else
1238 0 : mooseException("Temperature ", temperature, " is out of range in ", name(), ": inRegion()");
1239 :
1240 : // Determine the phase region that the (P, T) point lies in
1241 : unsigned int region;
1242 :
1243 403092 : if (temperature >= 273.15 && temperature <= 623.15)
1244 : {
1245 402905 : if (pressure > vaporPressure(temperature) && pressure <= 100.0e6)
1246 : region = 1;
1247 : else
1248 : region = 2;
1249 : }
1250 187 : else if (temperature > 623.15 && temperature <= 863.15)
1251 : {
1252 99 : if (pressure <= b23p(temperature))
1253 : region = 2;
1254 : else
1255 : region = 3;
1256 : }
1257 88 : else if (temperature > 863.15 && temperature <= 1073.15)
1258 : region = 2;
1259 : else
1260 : region = 5;
1261 :
1262 403092 : return region;
1263 : }
|
