GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/material/components/co2.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 37 45 82.2%
Functions: 5 6 83.3%
Branches: 39 68 57.4%

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-FileCopyrightInfo: 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 Components
10 * \brief A class for the CO2 fluid properties
11 */
12 #ifndef DUMUX_CO2_HH
13 #define DUMUX_CO2_HH
14
15 #include <cmath>
16 #include <iostream>
17
18 #include <dune/common/stdstreams.hh>
19 #include <dune/common/math.hh>
20
21 #include <dumux/common/exceptions.hh>
22
23 #include <dumux/material/constants.hh>
24 #include <dumux/material/components/base.hh>
25 #include <dumux/material/components/liquid.hh>
26 #include <dumux/material/components/gas.hh>
27
28 namespace Dumux {
29 namespace Components {
30
31 /*!
32 * \ingroup Components
33 * \brief A class for the CO2 fluid properties
34 *
35 * Under reservoir conditions, CO2 is typically in supercritical state. These
36 * properties can be provided in tabulated form, which is necessary for this
37 * component implementation. The template is passed through the fluidsystem
38 * brineco2fluidsystem.hh.
39 * Depending on the used tabulation, the fluidsystem can also be used for gaseous CO2
40 */
41
42 template <class Scalar, class CO2Tables>
43 class CO2
44 : public Components::Base<Scalar, CO2<Scalar, CO2Tables> >
45 , public Components::Liquid<Scalar, CO2<Scalar, CO2Tables> >
46 , public Components::Gas<Scalar, CO2<Scalar, CO2Tables> >
47 {
48 static const Scalar R;
49
50 static bool warningThrown;
51
52 public:
53 /*!
54 * \brief A human readable name for the CO2.
55 */
56 static std::string name()
57
24/48
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 1 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 1 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 1 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 1 times.
✗ Branch 34 not taken.
✓ Branch 36 taken 1 times.
✗ Branch 37 not taken.
✓ Branch 41 taken 1 times.
✗ Branch 42 not taken.
✓ Branch 44 taken 1 times.
✗ Branch 45 not taken.
✓ Branch 47 taken 1 times.
✗ Branch 48 not taken.
✓ Branch 50 taken 1 times.
✗ Branch 51 not taken.
✓ Branch 53 taken 1 times.
✗ Branch 54 not taken.
✓ Branch 56 taken 1 times.
✗ Branch 57 not taken.
✓ Branch 61 taken 1 times.
✗ Branch 62 not taken.
✓ Branch 64 taken 1 times.
✗ Branch 65 not taken.
✓ Branch 67 taken 1 times.
✗ Branch 68 not taken.
✓ Branch 70 taken 1 times.
✗ Branch 71 not taken.
✓ Branch 73 taken 1 times.
✗ Branch 74 not taken.
✓ Branch 76 taken 1 times.
✗ Branch 77 not taken.
44 { return "CO2"; }
58
59 /*!
60 * \brief The mass in \f$\mathrm{[kg/mol]}\f$ of one mole of CO2.
61 */
62 static constexpr Scalar molarMass()
63 { return 44e-3; /* [kg/mol] */ }
64
65 /*!
66 * \brief Returns the critical temperature \f$\mathrm{[K]}\f$ of CO2
67 */
68 static Scalar criticalTemperature()
69 { return 273.15 + 30.95; /* [K] */ }
70
71 /*!
72 * \brief Returns the critical pressure \f$\mathrm{[Pa]}\f$ of CO2
73 */
74 static Scalar criticalPressure()
75 { return 73.8e5; /* [Pa] */ }
76
77 /*!
78 * \brief Returns the temperature \f$\mathrm{[K]}\f$ at CO2's triple point.
79 */
80 static Scalar tripleTemperature()
81 { return 273.15 - 56.35; /* [K] */ }
82
83 /*!
84 * \brief Returns the pressure \f$\mathrm{[Pa]}\f$ at CO2's triple point.
85 */
86 static Scalar triplePressure()
87 { return 5.11e5; /* [N/m^2] */ }
88
89 /*!
90 * \brief Returns the minimal tabulated pressure \f$\mathrm{[Pa]}\f$ of the used table
91 */
92 static Scalar minTabulatedPressure()
93 { return CO2Tables::tabulatedEnthalpy.minPress(); /* [Pa] */ }
94
95 /*!
96 * \brief Returns the maximal tabulated pressure \f$\mathrm{[Pa]}\f$ of the used table
97 */
98 static Scalar maxTabulatedPressure()
99 { return CO2Tables::tabulatedEnthalpy.maxPress(); /* [Pa] */ }
100
101 /*!
102 * \brief Returns the minimal tabulated temperature \f$\mathrm{[K]}\f$ of the used table
103 */
104 static Scalar minTabulatedTemperature()
105 { return CO2Tables::tabulatedEnthalpy.minTemp(); /* [K] */ }
106
107 /*!
108 * \brief Returns the maximal tabulated temperature \f$\mathrm{[K]}\f$ of the used table
109 */
110 static Scalar maxTabulatedTemperature()
111 { return CO2Tables::tabulatedEnthalpy.maxTemp(); /* [K] */ }
112
113 /*!
114 * \brief Returns true if the gas phase is assumed to be ideal
115 */
116 static constexpr bool gasIsIdeal()
117 { return false; }
118
119 /*!
120 * \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of pure CO2
121 * at a given temperature.
122 * \param T the temperature \f$\mathrm{[K]}\f$
123 * See:
124 *
125 * R. Span and W. Wagner (1996, pp. 1509-1596) \cite span1996
126 */
127 static Scalar vaporPressure(Scalar T)
128 {
129 static const Scalar a[4] =
130 { -7.0602087, 1.9391218, -1.6463597, -3.2995634 };
131 static const Scalar t[4] =
132 { 1.0, 1.5, 2.0, 4.0 };
133
134 // this is on page 1524 of the reference
135 Scalar exponent = 0;
136 Scalar Tred = T/criticalTemperature();
137
138 using std::pow;
139 for (int i = 0; i < 4; ++i)
140 exponent += a[i]*pow(1 - Tred, t[i]);
141 exponent *= 1.0/Tred;
142
143 using std::exp;
144 return exp(exponent)*criticalPressure();
145 }
146
147 /*!
148 * \brief Specific enthalpy of gaseous CO2 \f$\mathrm{[J/kg]}\f$.
149 * \param temperature the temperature \f$\mathrm{[K]}\f$
150 * \param pressure the pressure \f$\mathrm{[Pa]}\f$
151 */
152 92849396 static Scalar gasEnthalpy(Scalar temperature,
153 Scalar pressure)
154 {
155
4/6
✓ Branch 0 taken 92849364 times.
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 92849364 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 32 times.
92849396 if ((temperature < criticalTemperature() || pressure < criticalPressure()) && !warningThrown)
156 {
157 Dune::dwarn << "Subcritical values: Be aware to use "
158 <<"Tables with sufficient resolution!"<< std::endl;
159 warningThrown=true;
160 }
161 return
162 92849396 CO2Tables::tabulatedEnthalpy.at(temperature, pressure);
163 }
164
165 /*!
166 * \brief Specific enthalpy of liquid CO2 \f$\mathrm{[J/kg]}\f$.
167 * \param temperature the temperature \f$\mathrm{[K]}\f$
168 * \param pressure the pressure \f$\mathrm{[Pa]}\f$
169 */
170 46343538 static Scalar liquidEnthalpy(Scalar temperature,
171 Scalar pressure)
172 {
173
5/6
✓ Branch 0 taken 46343514 times.
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 46343514 times.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 23 times.
46343538 if ((temperature < criticalTemperature() || pressure < criticalPressure()) && !warningThrown)
174 {
175 1 Dune::dwarn << "Subcritical values: Be aware to use "
176 1 <<"Tables with sufficient resolution!"<< std::endl;
177 1 warningThrown=true;
178 }
179
180 46343538 return gasEnthalpy(temperature, pressure);
181 }
182
183 /*!
184 * \brief Specific internal energy of CO2 \f$\mathrm{[J/kg]}\f$.
185 * \param temperature the temperature \f$\mathrm{[K]}\f$
186 * \param pressure the pressure \f$\mathrm{[Pa]}\f$
187 */
188 static Scalar gasInternalEnergy(Scalar temperature,
189 Scalar pressure)
190 {
191 Scalar h = gasEnthalpy(temperature, pressure);
192 Scalar rho = gasDensity(temperature, pressure);
193
194 return h - (pressure / rho);
195 }
196
197 /*!
198 * \brief Specific internal energy of liquid CO2 \f$\mathrm{[J/kg]}\f$.
199 * \param temperature the temperature \f$\mathrm{[K]}\f$
200 * \param pressure the pressure \f$\mathrm{[Pa]}\f$
201 */
202 static Scalar liquidInternalEnergy(Scalar temperature,
203 Scalar pressure)
204 {
205 Scalar h = liquidEnthalpy(temperature, pressure);
206 Scalar rho = liquidDensity(temperature, pressure);
207
208 return h - (pressure / rho);
209 }
210
211 /*!
212 * \brief The density of CO2 at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$.
213 * \param temperature the temperature \f$\mathrm{[K]}\f$
214 * \param pressure the pressure \f$\mathrm{[Pa]}\f$
215 */
216 531440676 static Scalar gasDensity(Scalar temperature, Scalar pressure)
217 {
218
5/6
✓ Branch 0 taken 526122495 times.
✓ Branch 1 taken 5318181 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 526122495 times.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 5318178 times.
531440676 if ((temperature < criticalTemperature() || pressure < criticalPressure()) && !warningThrown)
219 {
220 3 Dune::dwarn << "Subcritical values: Be aware to use "
221 3 <<"Tables with sufficient resolution!"<< std::endl;
222 3 warningThrown=true;
223 }
224 531440676 return CO2Tables::tabulatedDensity.at(temperature, pressure);
225 }
226
227 /*!
228 * \brief The molar density of CO2 gas in \f$\mathrm{[mol/m^3]}\f$ at a given pressure and temperature.
229 *
230 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
231 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
232 *
233 */
234 static Scalar gasMolarDensity(Scalar temperature, Scalar pressure)
235 71188469 { return gasDensity(temperature, pressure)/molarMass(); }
236
237 /*!
238 * \brief The density of pure CO2 at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$.
239 * \param temperature the temperature \f$\mathrm{[K]}\f$
240 * \param pressure the pressure \f$\mathrm{[Pa]}\f$
241 */
242 static Scalar liquidDensity(Scalar temperature, Scalar pressure)
243 {
244 if ((temperature < criticalTemperature() || pressure < criticalPressure()) && !warningThrown)
245 {
246 Dune::dwarn << "Subcritical values: Be aware to use "
247 <<"Tables with sufficient resolution!"<< std::endl;
248 warningThrown=true;
249 }
250 return CO2Tables::tabulatedDensity.at(temperature, pressure);
251 }
252
253 /*!
254 * \brief The molar density of CO2 in \f$\mathrm{[mol/m^3]}\f$ at a given pressure and temperature.
255 *
256 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
257 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
258 *
259 */
260 static Scalar liquidMolarDensity(Scalar temperature, Scalar pressure)
261 { return liquidDensity(temperature, pressure)/molarMass(); }
262
263 /*!
264 * \brief The pressure of steam in \f$\mathrm{[Pa]}\f$ at a given density and temperature.
265 *
266 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
267 * \param density density of component in \f$\mathrm{[kg/m^3]}\f$
268 */
269 static Scalar gasPressure(Scalar temperature, Scalar density)
270 {
271 DUNE_THROW(NumericalProblem, "CO2::gasPressure()");
272 }
273
274 /*!
275 * \brief The pressure of liquid water in \f$\mathrm{[Pa]}\f$ at a given density and
276 * temperature.
277 *
278 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
279 * \param density density of component in \f$\mathrm{[kg/m^3]}\f$
280 */
281 static Scalar liquidPressure(Scalar temperature, Scalar density)
282 {
283 DUNE_THROW(NumericalProblem, "CO2::liquidPressure()");
284 }
285
286 /*!
287 * \brief Specific isobaric heat capacity of the component \f$\mathrm{[J/(kg*K)]}\f$ as a liquid.
288 * USE WITH CAUTION! Exploits enthalpy function with artificial increment
289 * of the temperature!
290 * Equation with which the specific heat capacity is calculated : \f$ c_p = \frac{dh}{dT}\f$
291 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
292 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
293 */
294 8 static Scalar liquidHeatCapacity(Scalar temperature, Scalar pressure)
295 {
296 //temperature difference :
297 8 Scalar dT = 1.; // 1K temperature increment
298 8 Scalar temperature2 = temperature+dT;
299
300 // enthalpy difference
301 8 Scalar hold = liquidEnthalpy(temperature, pressure);
302 8 Scalar hnew = liquidEnthalpy(temperature2, pressure);
303 8 Scalar dh = hold-hnew;
304
305 //specific heat capacity
306 8 return dh/dT ;
307 }
308
309
310 /*!
311 * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of CO2.
312 * Equations given in: - Vesovic et al., 1990
313 * - Fenhour et al., 1998
314 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
315 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
316 */
317 145035987 static Scalar gasViscosity(Scalar temperature, Scalar pressure)
318 {
319 static const double a0 = 0.235156;
320 static const double a1 = -0.491266;
321 static const double a2 = 5.211155E-2;
322 static const double a3 = 5.347906E-2;
323 static const double a4 = -1.537102E-2;
324
325 static const double d11 = 0.4071119E-2;
326 static const double d21 = 0.7198037E-4;
327 static const double d64 = 0.2411697E-16;
328 static const double d81 = 0.2971072E-22;
329 static const double d82 = -0.1627888E-22;
330
331 static const double ESP = 251.196;
332
333 double mu0, SigmaStar, TStar;
334 double dmu, rho;
335 double visco_CO2;
336
337
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 145035987 times.
145035987 if(temperature < 275.) // regularisation
338 {
339 temperature = 275;
340 Dune::dgrave << "Temperature below 275K in viscosity function:"
341 << "Regularizing temperature to 275K. " << std::endl;
342 }
343
344
345 145035987 TStar = temperature/ESP;
346
347 /* mu0: viscosity in zero-density limit */
348 using std::exp;
349 using std::log;
350 using std::sqrt;
351 290071974 SigmaStar = exp(a0 + a1*log(TStar)
352 145035987 + a2*log(TStar)*log(TStar)
353 145035987 + a3*log(TStar)*log(TStar)*log(TStar)
354 145035987 + a4*log(TStar)*log(TStar)*log(TStar)*log(TStar) );
355 145035987 mu0 = 1.00697*sqrt(temperature) / SigmaStar;
356
357 /* dmu : excess viscosity at elevated density */
358 145035987 rho = gasDensity(temperature, pressure); /* CO2 mass density [kg/m^3] */
359
360 using Dune::power;
361 145035987 dmu = d11*rho + d21*rho*rho + d64*power(rho,6)/(TStar*TStar*TStar)
362 290071974 + d81*power(rho,8) + d82*power(rho,8)/TStar;
363
364 145035987 visco_CO2 = (mu0 + dmu)/1.0E6; /* conversion to [Pa s] */
365
366 145035987 return visco_CO2;
367 }
368
369 /*!
370 * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of pure CO2.
371 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
372 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
373 */
374 static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
375 {
376 // no difference for supercritical CO2
377 return gasViscosity(temperature, pressure);
378 }
379
380 /*!
381 * \brief Thermal conductivity \f$\mathrm{[[W/(m*K)]}\f$ of CO2.
382 *
383 * Thermal conductivity of CO2 at T=20°C, see:
384 * http://www.engineeringtoolbox.com/carbon-dioxide-d_1000.html
385 *
386 * \param temperature absolute temperature in \f$\mathrm{[K]}\f$
387 * \param pressure of the phase in \f$\mathrm{[Pa]}\f$
388 */
389 static Scalar gasThermalConductivity(Scalar temperature, Scalar pressure)
390 {
391 return 0.087;
392 }
393 };
394
395 template <class Scalar, class CO2Tables>
396 const Scalar CO2<Scalar, CO2Tables>::R = Constants<Scalar>::R;
397
398 template <class Scalar, class CO2Tables>
399 bool CO2<Scalar, CO2Tables>::warningThrown = false;
400
401 } // end namespace Components
402
403 } // end namespace Dumux
404
405 #endif
406