GCC Code Coverage Report


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