GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/material/components/iapws/common.hh
Date: 2025-06-14 19:21:29
Exec Total Coverage
Lines: 67 67 100.0%
Functions: 2 2 100.0%
Branches: 10 12 83.3%

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 * \file
9 * \ingroup IAPWS
10 * \brief Implements relations common for all regions of the IAPWS '97
11 * 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_COMMON_HH
19 #define DUMUX_IAPWS_COMMON_HH
20
21 #include <cmath>
22 #include <iostream>
23
24 #include <dune/common/math.hh>
25
26 #include <dumux/material/constants.hh>
27
28 namespace Dumux::IAPWS {
29
30 /*!
31 * \ingroup IAPWS
32 * \brief Implements relations which are common for all regions of the IAPWS '97
33 * formulation.
34 *
35 * \tparam Scalar The type used for scalar values
36 *
37 * See:
38 *
39 * IAPWS: "Revised Release on the IAPWS Industrial Formulation
40 * 1997 for the Thermodynamic Properties of Water and Steam",
41 * http://www.iapws.org/relguide/IF97-Rev.pdf
42 */
43 template <class Scalar>
44 class Common
45 {
46 public:
47 //! The molar mass of water \f$\mathrm{[kg/mol]}\f$
48 static constexpr Scalar molarMass = 18.01518e-3;
49
50 //! Specific gas constant of water \f$\mathrm{[J/(kg*K)]}\f$
51 static constexpr Scalar Rs = Constants<Scalar>::R / molarMass;
52
53 //! Critical temperature of water \f$\mathrm{[K]}\f$
54 static constexpr Scalar criticalTemperature = 647.096;
55
56 //! Critical pressure of water \f$\mathrm{[Pa]}\f$
57 static constexpr Scalar criticalPressure = 22.064e6;
58
59 //! Critical molar volume of water \f$\mathrm{[m^3/mol]}\f$
60 static constexpr Scalar criticalMolarVolume = molarMass / 322.0;
61
62 //! The acentric factor of water \f$\mathrm{[-]}\f$
63 static constexpr Scalar acentricFactor = 0.344;
64
65 //! Density of water at the critical point \f$\mathrm{[kg/m^3]}\f$
66 static constexpr Scalar criticalDensity = 322;
67
68 //! Triple temperature of water \f$\mathrm{[K]}\f$
69 static constexpr Scalar tripleTemperature = 273.16;
70
71 //! Triple pressure of water \f$\mathrm{[Pa]}\f$
72 static constexpr Scalar triplePressure = 611.657;
73
74 /*!
75 * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of pure water.
76 *
77 * This relation is valid for all regions of the IAPWS '97
78 * formulation.
79 *
80 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
81 * \param rho density of component in \f$\mathrm{[kg/m^3]}\f$
82 *
83 * See:
84 *
85 * IAPWS: "Release on the IAPWS Formulation 2008 for the Viscosity
86 * of Ordinary Water Substance", http://www.iapws.org/relguide/visc.pdf
87 */
88 7999245 static Scalar viscosity(Scalar temperature, Scalar rho)
89 {
90 7999245 const Scalar rhoBar = rho/322.0;
91 7999245 const Scalar TBar = temperature/criticalTemperature;
92
93 // muBar = muBar_1
94 7999245 constexpr Scalar Hij[6][7] = {
95 { 5.20094e-1, 2.22531e-1,-2.81378e-1, 1.61913e-1,-3.25372e-2, 0, 0 },
96 { 8.50895e-2, 9.99115e-1,-9.06851e-1, 2.57399e-1, 0, 0, 0 },
97 {-1.08374 , 1.88797 ,-7.72479e-1, 0, 0, 0, 0 },
98 {-2.89555e-1, 1.26613 ,-4.89837e-1, 0, 6.98452e-2, 0,-4.35673e-3 },
99 { 0, 0,-2.57040e-1, 0, 0, 8.72102e-3, 0 },
100 { 0, 1.20573e-1, 0, 0, 0, 0,-5.93264e-4 }
101 };
102
103 7999245 Scalar tmp, tmp2, tmp3 = 1;
104 7999245 Scalar muBar = 0;
105
2/2
✓ Branch 0 taken 47995470 times.
✓ Branch 1 taken 7999245 times.
55994715 for (int i = 0; i <= 5; ++i) {
106 tmp = 0;
107 tmp2 = 1;
108
2/2
✓ Branch 0 taken 335968290 times.
✓ Branch 1 taken 47995470 times.
383963760 for (int j = 0; j <= 6; ++j) {
109 335968290 tmp += Hij[i][j]*tmp2;
110 335968290 tmp2 *= (rhoBar - 1);
111 }
112 47995470 muBar += tmp3 * tmp;
113 47995470 tmp3 *= 1.0/TBar - 1;
114 }
115 using std::exp;
116 7999245 muBar *= rhoBar;
117 7999245 muBar = exp(muBar);
118
119 // muBar *= muBar_0
120 using std::sqrt;
121 7999245 muBar *= 100*sqrt(TBar);
122 7999245 constexpr Scalar H[4] = {
123 1.67752, 2.20462, 0.6366564, -0.241605
124 };
125
126 7999245 tmp = 0, tmp2 = 1;
127
2/2
✓ Branch 0 taken 31996980 times.
✓ Branch 1 taken 7999245 times.
39996225 for (int i = 0; i < 4; ++i) {
128 31996980 tmp += H[i]/tmp2;
129 31996980 tmp2 *= TBar;
130 }
131 7999245 muBar /= tmp;
132
133 7999245 return 1e-6*muBar;
134 }
135
136 /*!
137 * \brief Thermal conductivity \f$\mathrm{[[W/(m*K)]}\f$ water (IAPWS) .
138 *
139 * Implementation taken from:
140 * freesteam - IAPWS-IF97 steam tables library
141 * copyright (C) 2004-2009 John Pye
142 *
143 * Appendix B: Recommended Interpolating equation for Industrial Use
144 * see "Revised Release on the IAPWS Formulation 1985 for the Thermal Conductivity of Ordinary Water Substance",
145 * https://doc.modelica.org/Modelica%204.0.0/Resources/Documentation/Media/Water/IF97documentation/thcond.pdf
146 * \cite IAPWS_ThCond
147 *
148 * TODO: implement new formulation (and update dumux.bib entry then) presented in
149 * "Release on the IAPWS Formulation 2011 for the Thermal Conductivity of Ordinary Water Substance",
150 * https://iapws.org/relguide/ThCond.pdf
151 *
152 * \param T absolute temperature in \f$\mathrm{[K]}\f$
153 * \param rho density of water in \f$\mathrm{[kg/m^3]}\f$
154 */
155 6748568 static Scalar thermalConductivityIAPWS(const Scalar T, const Scalar rho)
156 {
157 6748568 constexpr Scalar thcond_tstar = 647.26 ;
158 6748568 constexpr Scalar thcond_rhostar = 317.7 ;
159 /*static constexpr Scalar thcond_kstar = 1.0 ;*/
160
161 6748568 constexpr Scalar thcond_b0 = -0.397070 ;
162 6748568 constexpr Scalar thcond_b1 = 0.400302 ;
163 6748568 constexpr Scalar thcond_b2 = 1.060000 ;
164 6748568 constexpr Scalar thcond_B1 = -0.171587 ;
165 6748568 constexpr Scalar thcond_B2 = 2.392190 ;
166
167 6748568 constexpr Scalar thcond_c1 = 0.642857 ;
168 6748568 constexpr Scalar thcond_c2 = -4.11717 ;
169 6748568 constexpr Scalar thcond_c3 = -6.17937 ;
170 6748568 constexpr Scalar thcond_c4 = 0.00308976 ;
171 6748568 constexpr Scalar thcond_c5 = 0.0822994 ;
172 6748568 constexpr Scalar thcond_c6 = 10.0932 ;
173
174 6748568 constexpr Scalar thcond_d1 = 0.0701309 ;
175 6748568 constexpr Scalar thcond_d2 = 0.0118520 ;
176 6748568 constexpr Scalar thcond_d3 = 0.00169937 ;
177 6748568 constexpr Scalar thcond_d4 = -1.0200 ;
178 6748568 constexpr unsigned int thcond_a_count = 4;
179 6748568 constexpr Scalar thcond_a[thcond_a_count] = {
180 0.0102811
181 ,0.0299621
182 ,0.0156146
183 ,-0.00422464
184 };
185
186 6748568 const Scalar Tbar = T / thcond_tstar;
187 6748568 const Scalar rhobar = rho / thcond_rhostar;
188
189 /* fast implementation... minimised calls to 'pow' routine... */
190 using std::sqrt;
191 6748568 const Scalar Troot = sqrt(Tbar);
192 6748568 Scalar Tpow = Troot;
193 6748568 Scalar lam = 0;
194
195
2/2
✓ Branch 0 taken 26994272 times.
✓ Branch 1 taken 6748568 times.
33742840 for(unsigned int k = 0; k < thcond_a_count; ++k) {
196 26994272 lam += thcond_a[k] * Tpow;
197 26994272 Tpow *= Tbar;
198 }
199
200 using std::exp;
201 6748568 lam += thcond_b0 + thcond_b1
202 6748568 * rhobar + thcond_b2
203 6748568 * exp(thcond_B1 * ((rhobar + thcond_B2)*(rhobar + thcond_B2)));
204
205 using std::abs;
206 using std::pow;
207 using Dune::power;
208
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6748568 times.
6748568 const Scalar DTbar = abs(Tbar - 1) + thcond_c4;
209 6748568 const Scalar DTbarpow = pow(DTbar, 3./5);
210 6748568 const Scalar Q = 2. + thcond_c5 / DTbarpow;
211
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6748568 times.
6748568 const Scalar S = (Tbar >= 1) ? 1. / DTbar : thcond_c6 / DTbarpow;
212
213 6748568 const Scalar rhobar18 = pow(rhobar, 1.8);
214 6748568 const Scalar rhobarQ = pow(rhobar, Q);
215
216 6748568 lam +=
217 6748568 (thcond_d1 / power(Tbar,10) + thcond_d2) * rhobar18 *
218 6748568 exp(thcond_c1 * (1 - rhobar * rhobar18))
219 6748568 + thcond_d3 * S * rhobarQ *
220 6748568 exp((Q/(1+Q))*(1 - rhobar*rhobarQ))
221 6748568 + thcond_d4 *
222 6748568 exp(thcond_c2 * power(Troot,3) + thcond_c3 / power(rhobar,5));
223 6748568 return /*thcond_kstar * */ lam;
224 }
225 };
226
227 } // end namespace Dumux::IAPWS
228
229 #endif
230