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 |