GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/material/components/iapws/region2.hh
Date: 2025-06-14 19:21:29
Exec Total Coverage
Lines: 74 75 98.7%
Functions: 8 8 100.0%
Branches: 14 54 25.9%

Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7
8 /*!
9 * \file
10 * \ingroup IAPWS
11 * \brief Implements the equations for region 2 of the IAPWS '97 formulation.
12 * See:
13 *
14 * IAPWS: "Revised Release on the IAPWS Industrial Formulation
15 * 1997 for the Thermodynamic Properties of Water and Steam",
16 * http://www.iapws.org/relguide/IF97-Rev.pdf
17 */
18 #ifndef DUMUX_IAPWS_REGION2_HH
19 #define DUMUX_IAPWS_REGION2_HH
20
21 #include <cmath>
22 #include <iostream>
23 #include <dumux/common/exceptions.hh>
24
25 namespace Dumux {
26 namespace IAPWS {
27
28 /*!
29 * \ingroup IAPWS
30 * \brief Implements the equations for region 2 of the IAPWS '97 formulation.
31 * \tparam Scalar The type used for scalar values
32 * See:
33 *
34 * IAPWS: "Revised Release on the IAPWS Industrial Formulation
35 * 1997 for the Thermodynamic Properties of Water and Steam",
36 * http://www.iapws.org/relguide/IF97-Rev.pdf
37 */
38 template <class Scalar>
39 class Region2
40 {
41 public:
42 /*!
43 * \brief Returns true if IAPWS region 2 applies for a
44 * (temperature, pressure) pair.
45 *
46 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
47 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
48 * \param propertyName the name for which property the check is performed
49 */
50 61849172 static void checkValidityRange(Scalar temperature, Scalar pressure,
51 const std::string& propertyName = "This property")
52 {
53 // actually this is:
54 /* (273.15 <= temperature && temperature <= 623.15 && pressure <= vaporPressure(temperature)) ||
55 (623.15 < temperature && temperature <= 863.15 && pressure <= auxPressure(temperature)) ||
56 (863.15 < temperature && temperature <= 1073.15 && pressure < 100e6); */
57
2/6
✓ Branch 0 taken 61849172 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 61849172 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
61849172 if ((temperature <= 623.15 && pressure <= 100e6) ||
58 (temperature > 623.15 && temperature <= 1073.15 && pressure <= 16.532e6 ))
59 61849172 return;
60
61
0/32
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
61849172 DUNE_THROW(NumericalProblem,
62 propertyName << " of steam is only implemented for temperatures below 623.15K and "
63 "pressures below 100MPa. (T=" << temperature << ", p=" << pressure << ")");
64 }
65
66 /*!
67 * \brief Returns the reduced temperature (dimensionless) for IAPWS region 2.
68 *
69 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
70 */
71 78111507 static constexpr Scalar tau(Scalar temperature)
72 18123006 { return 540.0 / temperature; }
73
74 /*!
75 * \brief Returns the derivative of the reduced temperature to the
76 * temperature for IAPWS region 2.
77 *
78 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
79 */
80 static constexpr Scalar dTau_dt(Scalar temperature)
81 { return - 540.0 / (temperature*temperature); }
82
83 /*!
84 * \brief Returns the reduced pressure (dimensionless) for IAPWS region 2.
85 *
86 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
87 */
88 116859272 static constexpr Scalar pi(Scalar pressure)
89 41844020 { return pressure / 1e6; }
90
91 /*!
92 * \brief Returns the derivative of the reduced pressure to the
93 * pressure for IAPWS region 2 in \f$\mathrm{[1/Pa]}\f$.
94 *
95 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
96 */
97 static constexpr Scalar dPi_dp(Scalar pressure)
98 { return 1.0 / 1e6; }
99
100 /*!
101 * \brief Returns the derivative of the pressure to the
102 * reduced pressure for IAPWS region 2 (dimensionless).
103 *
104 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
105 */
106 static constexpr Scalar dp_dPi(Scalar pressure)
107 { return 1e6; }
108
109 /*!
110 * \brief The Gibbs free energy for IAPWS region 2 (i.e. sub-critical
111 * steam) (dimensionless).
112 *
113 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
114 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
115 *
116 * IAPWS: "Revised Release on the IAPWS Industrial Formulation
117 * 1997 for the Thermodynamic Properties of Water and Steam",
118 * http://www.iapws.org/relguide/IF97-Rev.pdf
119 */
120 static Scalar gamma(Scalar temperature, Scalar pressure)
121 {
122 Scalar tau = tau(temperature); /* reduced temperature */
123 Scalar pi = pi(pressure); /* reduced pressure */
124
125 Scalar result;
126
127 // ideal gas part
128 using std::pow;
129 result = ln(pi);
130 for (int i = 0; i < 9; ++i)
131 result += n_g(i)*pow(tau, J_g(i));
132
133 // residual part
134 for (int i = 0; i < 43; ++i)
135 result += n_r(i)*
136 pow(pi, I_r(i))*
137 pow(tau - 0.5, J_r(i));
138 return result;
139 }
140
141 /*!
142 * \brief The partial derivative of the Gibbs free energy to the
143 * normalized temperature for IAPWS region 2 (i.e. sub-critical
144 * steam) dimensionless).
145 *
146 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
147 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
148 *
149 * IAPWS: "Revised Release on the IAPWS Industrial Formulation
150 * 1997 for the Thermodynamic Properties of Water and Steam",
151 * http://www.iapws.org/relguide/IF97-Rev.pdf
152 */
153 6466527 static Scalar dGamma_dTau(Scalar temperature, Scalar pressure)
154 {
155 6466527 Scalar tau_ = tau(temperature); /* reduced temperature */
156 6466527 Scalar pi_ = pi(pressure); /* reduced pressure */
157
158 // ideal gas part
159 using std::pow;
160 6466527 Scalar result = 0;
161
2/2
✓ Branch 0 taken 58198743 times.
✓ Branch 1 taken 6466527 times.
64665270 for (int i = 0; i < 9; i++) {
162 58198743 result += n_g(i) *
163 58198743 J_g(i) *
164 58198743 pow(tau_, J_g(i) - 1);
165 }
166
167 // residual part
168
2/2
✓ Branch 0 taken 278060661 times.
✓ Branch 1 taken 6466527 times.
284527188 for (int i = 0; i < 43; i++) {
169 278060661 result += n_r(i) *
170 278060661 pow(pi_, I_r(i)) *
171 278060661 J_r(i) *
172 278060661 pow(tau_ - 0.5, J_r(i) - 1);
173 }
174
175 6466527 return result;
176 }
177
178 /*!
179 * \brief The partial derivative of the Gibbs free energy to the
180 * normalized pressure for IAPWS region 2 (i.e. sub-critical
181 * steam) (dimensionless).
182 *
183 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
184 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
185 *
186 * IAPWS: "Revised Release on the IAPWS Industrial Formulation
187 * 1997 for the Thermodynamic Properties of Water and Steam",
188 * http://www.iapws.org/relguide/IF97-Rev.pdf
189 */
190 55947911 static Scalar dGamma_dPi(Scalar temperature, Scalar pressure)
191 {
192 55947911 Scalar tau_ = tau(temperature); /* reduced temperature */
193 55947911 Scalar pi_ = pi(pressure); /* reduced pressure */
194
195 // ideal gas part
196 55947911 Scalar result = 1/pi_;
197
198 // residual part
199 using std::pow;
200
2/2
✓ Branch 0 taken 2405760173 times.
✓ Branch 1 taken 55947911 times.
2461708084 for (int i = 0; i < 43; i++) {
201 2405760173 result += n_r(i) *
202 2405760173 I_r(i) *
203 2405760173 pow(pi_, I_r(i) - 1) *
204 2405760173 pow(tau_ - 0.5, J_r(i));
205 }
206
207 55947911 return result;
208 }
209
210 /*!
211 * \brief The partial derivative of the Gibbs free energy to the
212 * normalized pressure and to the normalized temperature
213 * for IAPWS region 2 (i.e. sub-critical steam) (dimensionless).
214 *
215 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
216 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
217 *
218 * IAPWS: "Revised Release on the IAPWS Industrial Formulation
219 * 1997 for the Thermodynamic Properties of Water and Steam",
220 * http://www.iapws.org/relguide/IF97-Rev.pdf
221 */
222 3619312 static Scalar ddGamma_dTaudPi(Scalar temperature, Scalar pressure)
223 {
224 3619312 Scalar tau_ = tau(temperature); /* reduced temperature */
225 3619312 Scalar pi_ = pi(pressure); /* reduced pressure */
226
227 // ideal gas part
228 3619312 Scalar result = 0;
229
230 // residual part
231 using std::pow;
232
2/2
✓ Branch 0 taken 155630416 times.
✓ Branch 1 taken 3619312 times.
159249728 for (int i = 0; i < 43; i++) {
233 155630416 result += n_r(i) *
234 155630416 I_r(i) *
235 155630416 J_r(i) *
236 155630416 pow(pi_, I_r(i) - 1) *
237 155630416 pow(tau_ - 0.5, J_r(i) - 1);
238 }
239
240 3619312 return result;
241 }
242
243 /*!
244 * \brief The second partial derivative of the Gibbs free energy
245 * to the normalized pressure for IAPWS region 2
246 * (i.e. sub-critical steam) (dimensionless).
247 *
248 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
249 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
250 *
251 * IAPWS: "Revised Release on the IAPWS Industrial Formulation
252 * 1997 for the Thermodynamic Properties of Water and Steam",
253 * http://www.iapws.org/relguide/IF97-Rev.pdf
254 */
255 static Scalar ddGamma_ddPi(Scalar temperature, Scalar pressure)
256 {
257 Scalar tau_ = tau(temperature); /* reduced temperature */
258 Scalar pi_ = pi(pressure); /* reduced pressure */
259
260 // ideal gas part
261 Scalar result = -1/(pi_*pi_);
262
263 // residual part
264 using std::pow;
265 for (int i = 0; i < 43; i++) {
266 result += n_r(i) *
267 I_r(i) *
268 (I_r(i) - 1) *
269 pow(pi_, I_r(i) - 2) *
270 pow(tau_ - 0.5, J_r(i));
271 }
272
273 return result;
274 }
275
276 /*!
277 * \brief The second partial derivative of the Gibbs free energy to the
278 * normalized temperature for IAPWS region 2 (i.e. sub-critical
279 * steam) (dimensionless).
280 *
281 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
282 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
283 *
284 * IAPWS: "Revised Release on the IAPWS Industrial Formulation
285 * 1997 for the Thermodynamic Properties of Water and Steam",
286 * http://www.iapws.org/relguide/IF97-Rev.pdf
287 */
288 2805615 static Scalar ddGamma_ddTau(Scalar temperature, Scalar pressure)
289 {
290 2805615 Scalar tau_ = tau(temperature); /* reduced temperature */
291 2805615 Scalar pi_ = pi(pressure); /* reduced pressure */
292
293 // ideal gas part
294 using std::pow;
295 2805615 Scalar result = 0;
296
2/2
✓ Branch 0 taken 25250535 times.
✓ Branch 1 taken 2805615 times.
28056150 for (int i = 0; i < 9; i++) {
297 25250535 result += n_g(i) *
298 25250535 J_g(i) *
299 25250535 (J_g(i) - 1) *
300 25250535 pow(tau_, J_g(i) - 2);
301 }
302
303 // residual part
304 using std::pow;
305
2/2
✓ Branch 0 taken 120641445 times.
✓ Branch 1 taken 2805615 times.
123447060 for (int i = 0; i < 43; i++) {
306 120641445 result += n_r(i) *
307 120641445 pow(pi_, I_r(i)) *
308 120641445 J_r(i) *
309 120641445 (J_r(i) - 1.) *
310 120641445 pow(tau_ - 0.5, J_r(i) - 2.);
311 }
312
313 2805615 return result;
314 }
315
316 private:
317 83449278 static Scalar n_g(int i)
318 {
319 83449278 constexpr const Scalar n[9] = {
320 -0.96927686500217e1, 0.10086655968018e2, -0.56087911283020e-2,
321 0.71452738081455e-1, -0.40710498223928, 0.14240819171444e1,
322 -0.43839511319450e1, -0.28408632460772, 0.21268463753307e-1
323 };
324 83449278 return n[i];
325 }
326
327 2960092695 static Scalar n_r(int i)
328 {
329 2960092695 constexpr const Scalar n[43] = {
330 -0.17731742473213e-2, -0.17834862292358e-1, -0.45996013696365e-1,
331 -0.57581259083432e-1, -0.50325278727930e-1, -0.33032641670203e-4,
332 -0.18948987516315e-3, -0.39392777243355e-2, -0.43797295650573e-1,
333 -0.26674547914087e-4, 0.20481737692309e-7, 0.43870667284435e-6,
334 -0.32277677238570e-4, -0.15033924542148e-2, -0.40668253562649e-1,
335 -0.78847309559367e-9, 0.12790717852285e-7, 0.48225372718507e-6,
336 0.22922076337661e-5, -0.16714766451061e-10, -0.21171472321355e-2,
337 -0.23895741934104e2, -0.59059564324270e-17, -0.12621808899101e-5,
338 -0.38946842435739e-1, 0.11256211360459e-10, -0.82311340897998e1,
339 0.19809712802088e-7, 0.10406965210174e-18, -0.10234747095929e-12,
340 -0.10018179379511e-8, -0.80882908646985e-10, 0.10693031879409,
341 -0.33662250574171, 0.89185845355421e-24, 0.30629316876232e-12,
342 -0.42002467698208e-5, -0.59056029685639e-25, 0.37826947613457e-5,
343 -0.12768608934681e-14, 0.73087610595061e-28, 0.55414715350778e-16,
344 -0.94369707241210e-6
345 };
346 2960092695 return n[i];
347 }
348
349 2960092695 static Scalar I_r(int i)
350 {
351 2960092695 constexpr const short int I[43] = {
352 1, 1, 1,
353 1, 1, 2,
354 2, 2, 2,
355 2, 3, 3,
356 3, 3, 3,
357 4, 4, 4,
358 5, 6, 6,
359 6, 7, 7,
360 7, 8, 8,
361 9, 10, 10,
362 10, 16, 16,
363 18, 20, 20,
364 20, 21, 22,
365 23, 24, 24,
366 24
367 };
368 2960092695 return I[i];
369 }
370
371 83449278 static Scalar J_g(int i)
372 {
373 83449278 constexpr const short int J[9] = {
374 0, 1, -5,
375 -4, -3, -2,
376 -1, 2, 3
377 };
378 83449278 return J[i];
379 }
380
381 2960092695 static Scalar J_r(int i)
382 {
383 2960092695 constexpr const short int J[43] = {
384 0, 1, 2,
385 3, 6, 1,
386 2, 4, 7,
387 36, 0, 1,
388 3, 6, 35,
389 1, 2, 3,
390 7, 3, 16,
391 35, 0, 11,
392 25, 8, 36,
393 13, 4, 10,
394 14, 29, 50,
395 57, 20, 35,
396 48, 21, 53,
397 39, 26, 40,
398 58
399 };
400 2960092695 return J[i];
401 }
402
403 };
404
405 } // end namespace IAPWS
406 } // end namespace Dumux
407
408 #endif
409