GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/material/fluidsystems/2p1c.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 66 69 95.7%
Functions: 18 19 94.7%
Branches: 57 106 53.8%
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 FluidSystems
10 * \copybrief Dumux::FluidSystems::TwoPOneC
11 */
12 #ifndef DUMUX_2P_1C_FLUID_SYSTEM_HH
13 #define DUMUX_2P_1C_FLUID_SYSTEM_HH
14
15 #include <limits>
16 #include <cassert>
17 #include <iostream>
18
19 #include <dune/common/exceptions.hh>
20
21 #include <dumux/io/name.hh>
22
23 #include "base.hh"
24
25 namespace Dumux {
26 namespace FluidSystems {
27
28 /*!
29 * \ingroup FluidSystems
30 * \brief A two-phase fluid system with only one component.
31 */
32 template <class Scalar, class ComponentType>
33 class TwoPOneC
34 : public Base<Scalar, TwoPOneC<Scalar, ComponentType> >
35 {
36 using ThisType = TwoPOneC<Scalar, ComponentType>;
37 using Component = ComponentType;
38
39 public:
40 static constexpr int numPhases = 2; //!< Number of phases in the fluid system
41 static constexpr int numComponents = 1; //!< Number of components in the fluid system
42
43 static constexpr int liquidPhaseIdx = 0; //!< index of the liquid phase
44 static constexpr int gasPhaseIdx = 1; //!< index of the gas phase
45 static constexpr int phase0Idx = liquidPhaseIdx; //!< index of the first phase
46 static constexpr int phase1Idx = gasPhaseIdx; //!< index of the second phase
47
48 static constexpr int comp0Idx = 0; //!< index of the only component
49
50 /****************************************
51 * Fluid phase related static parameters
52 ****************************************/
53 /*!
54 * \brief Return the human readable name of a fluid phase
55 *
56 * \param phaseIdx The index of the fluid phase to consider
57 */
58 34 static std::string phaseName(int phaseIdx)
59 {
60
3/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 29 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
 34 static std::string name[] = {
61 std::string(IOName::liquidPhase()),
62 std::string(IOName::gaseousPhase()),
63 };
64
65
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 34 times.
 34 assert(0 <= phaseIdx && phaseIdx < numPhases);
66 34 return name[phaseIdx];
67 }
68
69 /*!
70 * \brief Returns whether the fluids are miscible
71 */
72 static constexpr bool isMiscible()
73 { return false; }
74
75 /*!
76 * \brief Return whether a phase is gaseous
77 *
78 * \param phaseIdx The index of the fluid phase to consider
79 */
80 static constexpr bool isGas(int phaseIdx)
81 {
82 assert(0 <= phaseIdx && phaseIdx < numPhases);
83 return phaseIdx == gasPhaseIdx;
84 }
85
86 /*!
87 * \brief Returns true if and only if a fluid phase is assumed to
88 * be an ideal mixture.
89 *
90 * We define an ideal mixture as a fluid phase where the fugacity
91 * coefficients of all components times the pressure of the phase
92 * are independent on the fluid composition. This assumption is true
93 * if Henry's law and Raoult's law apply. If you are unsure what
94 * this function should return, it is safe to return false. The
95 * only damage done will be (slightly) increased computation times
96 * in some cases.
97 *
98 * \param phaseIdx The index of the fluid phase to consider
99 */
100 static constexpr bool isIdealMixture(int phaseIdx)
101 {
102 assert(0 <= phaseIdx && phaseIdx < numPhases);
103 // we assume Henry's and Raoult's laws for the water phase and
104 // and no interaction between gas molecules of different
105 // components, so all phases are ideal mixtures!
106 return true;
107 }
108
109 /*!
110 * \brief Returns true if and only if a fluid phase is assumed to
111 * be compressible.
112 *
113 * Compressible means that the partial derivative of the density
114 * to the fluid pressure is always larger than zero.
115 *
116 * \param phaseIdx The index of the fluid phase to consider
117 */
118 static constexpr bool isCompressible(int phaseIdx)
119 {
120 assert(0 <= phaseIdx && phaseIdx < numPhases);
121 // gases are always compressible
122 if (phaseIdx == gasPhaseIdx)
123 return true;
124 // the component decides for the liquid phase...
125 return Component::liquidIsCompressible();
126 }
127
128 /*!
129 * \brief Returns true if and only if a fluid phase is assumed to
130 * be an ideal gas.
131 *
132 * \param phaseIdx The index of the fluid phase to consider
133 */
134 static constexpr bool isIdealGas(int phaseIdx)
135 {
136 assert(0 <= phaseIdx && phaseIdx < numPhases);
137
138 if (phaseIdx == gasPhaseIdx)
139 // let the components decide
140 return Component::gasIsIdeal();
141 return false; // not a gas
142 }
143
144 /****************************************
145 * Component related static parameters
146 ****************************************/
147
148 /*!
149 * \brief Return the human readable name of a component
150 *
151 * \param compIdx The index of the component to consider
152 */
153 1 static std::string componentName(int compIdx)
154 {
155
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 assert(0 <= compIdx && compIdx < numComponents);
156 1 return Component::name();
157 }
158
159 /*!
160 * \brief Return the molar mass of a component in [kg/mol].
161 *
162 * \param compIdx The index of the component to consider
163 */
164 static Scalar molarMass(int compIdx)
165 {
166 assert(0 <= compIdx && compIdx < numComponents);
167 return Component::molarMass();
168 }
169
170 /*!
171 * \brief Critical temperature of a component [K].
172 *
173 * \param compIdx The index of the component to consider
174 */
175 static Scalar criticalTemperature(int compIdx)
176 {
177 assert(0 <= compIdx && compIdx < numComponents);
178 return Component::criticalTemperature();
179 }
180
181 /*!
182 * \brief Critical pressure of a component [Pa].
183 *
184 * \param compIdx The index of the component to consider
185 */
186 static Scalar criticalPressure(int compIdx)
187 {
188 assert(0 <= compIdx && compIdx < numComponents);
189 return Component::criticalPressure();
190 }
191
192 /*!
193 * \brief Molar volume of a component at the critical point [m^3/mol].
194 *
195 * \param compIdx The index of the component to consider
196 */
197 static Scalar criticalMolarVolume(int compIdx)
198 {
199 DUNE_THROW(Dune::NotImplemented,
200 "TwoPOneC::criticalMolarVolume()");
201 return 0;
202 }
203
204 /*!
205 * \brief The acentric factor of a component [].
206 *
207 * \param compIdx The index of the component to consider
208 */
209 static Scalar acentricFactor(int compIdx)
210 {
211 assert(0 <= compIdx && compIdx < numComponents);
212 return Component::acentricFactor();
213 }
214
215 /****************************************
216 * thermodynamic relations
217 ****************************************/
218
219 /*!
220 * \brief Initialize the fluid system's static parameters generically
221 *
222 * If a tabulated H2O component is used, we do our best to create
223 * tables that always work.
224 */
225 static void init()
226 {
227 5 init(/*tempMin=*/273.15,
228 /*tempMax=*/623.15,
229 /*numTemp=*/100,
230 /*pMin=*/0.0,
231 /*pMax=*/20e6,
232 /*numP=*/200);
233 }
234
235 /*!
236 * \brief Initialize the fluid system's static parameters using
237 * problem specific temperature and pressure ranges
238 *
239 * \param tempMin The minimum temperature used for tabulation of water [K]
240 * \param tempMax The maximum temperature used for tabulation of water [K]
241 * \param nTemp The number of ticks on the temperature axis of the table of water
242 * \param pressMin The minimum pressure used for tabulation of water [Pa]
243 * \param pressMax The maximum pressure used for tabulation of water [Pa]
244 * \param nPress The number of ticks on the pressure axis of the table of water
245 */
246 ✗ static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
247 Scalar pressMin, Scalar pressMax, unsigned nPress)
248 {
249 if (Component::isTabulated)
250 {
251
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 Component::init(tempMin, tempMax, nTemp,
252 pressMin, pressMax, nPress);
253 }
254 ✗ }
255
256 using Base<Scalar, ThisType>::density;
257 //! \copydoc Base<Scalar,ThisType>::density(const FluidState&,int)
258 template <class FluidState>
259 2421212 static Scalar density(const FluidState &fluidState,
260 int phaseIdx)
261 {
262
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2421212 times.
 2421212 assert(0 <= phaseIdx && phaseIdx < numPhases);
263
264
2/2
✓ Branch 0 taken 1210607 times.
✓ Branch 1 taken 1210603 times.
 2421212 Scalar temperature = fluidState.temperature(phaseIdx);
265
2/2
✓ Branch 0 taken 1210607 times.
✓ Branch 1 taken 1210603 times.
 2421212 Scalar pressure = fluidState.pressure(phaseIdx);
266
267 // liquid phase
268
2/2
✓ Branch 0 taken 1210608 times.
✓ Branch 1 taken 1210604 times.
 2421212 if (phaseIdx == liquidPhaseIdx) {
269 1210608 return Component::liquidDensity(temperature, pressure);
270 }
271 else if (phaseIdx == gasPhaseIdx)// gas phase
272 {
273 1210604 return Component::gasDensity(temperature, pressure);
274 }
275 else
276 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index.");
277 }
278
279 using Base<Scalar, ThisType>::molarDensity;
280 //! \copydoc Base<Scalar,ThisType>::molarDensity(const FluidState&,int)
281 template <class FluidState>
282 2421208 static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
283 {
284
2/2
✓ Branch 0 taken 1210603 times.
✓ Branch 1 taken 1210603 times.
 2421208 Scalar temperature = fluidState.temperature(phaseIdx);
285
2/2
✓ Branch 0 taken 1210603 times.
✓ Branch 1 taken 1210603 times.
 2421208 Scalar pressure = fluidState.temperature(phaseIdx);
286
2/2
✓ Branch 0 taken 1210604 times.
✓ Branch 1 taken 1210604 times.
 2421208 if (phaseIdx == liquidPhaseIdx)
287 1210604 return Component::liquidMolarDensity(temperature, pressure);
288
1/2
✓ Branch 0 taken 1210604 times.
✗ Branch 1 not taken.
 1210604 else if (phaseIdx == gasPhaseIdx)
289 1210604 return Component::gasMolarDensity(temperature, pressure);
290 else
291 ✗ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index.");
292 }
293
294 using Base<Scalar, ThisType>::viscosity;
295 //! \copydoc Base<Scalar,ThisType>::viscosity(const FluidState&,int)
296 template <class FluidState>
297 2421208 static Scalar viscosity(const FluidState &fluidState,
298 int phaseIdx)
299 {
300
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2421208 times.
 2421208 assert(0 <= phaseIdx && phaseIdx < numPhases);
301
302
2/2
✓ Branch 0 taken 1210603 times.
✓ Branch 1 taken 1210603 times.
 2421208 Scalar temperature = fluidState.temperature(phaseIdx);
303
2/2
✓ Branch 0 taken 1210603 times.
✓ Branch 1 taken 1210603 times.
 2421208 Scalar pressure = fluidState.pressure(phaseIdx);
304
305 // liquid phase
306
2/2
✓ Branch 0 taken 1210604 times.
✓ Branch 1 taken 1210604 times.
 2421208 if (phaseIdx == liquidPhaseIdx)
307 1210604 return Component::liquidViscosity(temperature, pressure);
308 else if (phaseIdx == gasPhaseIdx) // gas phase
309 1210604 return Component::gasViscosity(temperature, pressure);
310 else
311 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index.");
312 }
313
314 /*!
315 * \brief calculate the temperature of vapor at a given pressure on the vapor pressure curve.
316 *
317 * \param fluidState An arbitrary fluid state
318 * \param phaseIdx The index of the fluid phase to consider
319 */
320 template <class FluidState>
321 2643 static Scalar vaporTemperature(const FluidState &fluidState,
322 const unsigned int phaseIdx)
323 {
324
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2643 times.
 2643 assert(0 <= phaseIdx && phaseIdx < numPhases);
325 99951 Scalar pressure = fluidState.pressure(gasPhaseIdx);
326
327
2/4
✓ Branch 2 taken 47 times.
✓ Branch 3 taken 97261 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 99951 return Component::vaporTemperature(pressure);
328 }
329
330 using Base<Scalar, ThisType>::fugacityCoefficient;
331 /*!
332 * \brief Calculate the fugacity coefficient [-] of an individual
333 * component in a fluid phase
334 *
335 * The fugacity coefficient \f$\phi^\kappa_\alpha\f$ of
336 * component \f$\kappa\f$ in phase \f$\alpha\f$ is connected to
337 * the fugacity \f$f^\kappa_\alpha\f$ and the component's mole
338 * fraction \f$x^\kappa_\alpha\f$ by means of the relation
339 *
340 * \f[
341 f^\kappa_\alpha = \phi^\kappa_\alpha\;x^\kappa_\alpha\;p_\alpha
342 \f]
343 * where \f$p_\alpha\f$ is the pressure of the fluid phase.
344 *
345 * The quantity "fugacity" itself is just an other way to express
346 * the chemical potential \f$\zeta^\kappa_\alpha\f$ of the
347 * component. It is defined via
348 *
349 * \f[
350 f^\kappa_\alpha := \exp\left\{\frac{\zeta^\kappa_\alpha}{k_B T_\alpha} \right\}
351 \f]
352 * where \f$k_B = 1.380\cdot10^{-23}\;J/K\f$ is the Boltzmann constant.
353 *
354 * \param fluidState An arbitrary fluid state
355 * \param phaseIdx The index of the fluid phase to consider
356 * \param compIdx The index of the component to consider
357 */
358 template <class FluidState>
359 2 static Scalar fugacityCoefficient(const FluidState &fluidState,
360 int phaseIdx,
361 int compIdx)
362 {
363
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 assert(0 <= phaseIdx && phaseIdx < numPhases);
364
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 assert(0 <= compIdx && compIdx < numComponents);
365
366 2 Scalar temperature = fluidState.temperature(phaseIdx);
367 2 Scalar pressure = fluidState.pressure(phaseIdx);
368
369 // liquid phase
370
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if (phaseIdx == liquidPhaseIdx)
371 1 return Component::vaporPressure(temperature)/pressure;
372
373 // for the gas phase, assume an ideal gas when it comes to
374 // fugacity (-> fugacity == partial pressure)
375 else if (phaseIdx == gasPhaseIdx)
376 return 1.0;
377
378 else
379 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index.");
380 }
381
382
383 using Base<Scalar, ThisType>::diffusionCoefficient;
384 //! \copydoc Base<Scalar,ThisType>::diffusionCoefficient(const FluidState&,int,int)
385 template <class FluidState>
386 2 static Scalar diffusionCoefficient(const FluidState &fluidState,
387 int phaseIdx,
388 int compIdx)
389 {
390 // TODO!
391
7/16
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 2 times.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 2 times.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
 22 DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients");
392 }
393
394 using Base<Scalar, ThisType>::binaryDiffusionCoefficient;
395 //! \copydoc Base<Scalar,ThisType>::binaryDiffusionCoefficient(const FluidState&,int,int,int)
396 template <class FluidState>
397 2 static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
398 int phaseIdx,
399 int compIIdx,
400 int compJIdx)
401
402
7/16
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 2 times.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 2 times.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
 22 { DUNE_THROW(Dune::NotImplemented, "Binary Diffusion coefficients"); }
403
404 using Base<Scalar, ThisType>::enthalpy;
405 //! \copydoc Base<Scalar,ThisType>::enthalpy(const FluidState&,int)
406 template <class FluidState>
407 2421208 static Scalar enthalpy(const FluidState &fluidState,
408 int phaseIdx)
409 {
410
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2421208 times.
 2421208 assert(0 <= phaseIdx && phaseIdx < numPhases);
411
412 // liquid phase
413
2/2
✓ Branch 0 taken 1210604 times.
✓ Branch 1 taken 1210604 times.
 2421208 if (phaseIdx == liquidPhaseIdx)
414 3631810 return Component::liquidEnthalpy(fluidState.temperature(phaseIdx),
415 1210604 fluidState.pressure(phaseIdx));
416
417 else if (phaseIdx == gasPhaseIdx) // gas phase
418 3631810 return Component::gasEnthalpy(fluidState.temperature(phaseIdx),
419 1210604 fluidState.pressure(phaseIdx));
420
421 else
422 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index.");
423 }
424
425 using Base<Scalar, ThisType>::thermalConductivity;
426 //! \copydoc Base<Scalar,ThisType>::thermalConductivity(const FluidState&,int)
427 template <class FluidState>
428 2421204 static Scalar thermalConductivity(const FluidState &fluidState,
429 const int phaseIdx)
430 {
431
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2421204 times.
 2421204 assert(0 <= phaseIdx && phaseIdx < numPhases);
432 // liquid phase
433
2/2
✓ Branch 0 taken 1210604 times.
✓ Branch 1 taken 1210600 times.
 2421204 if (phaseIdx == liquidPhaseIdx)
434 3631810 return Component::liquidThermalConductivity(fluidState.temperature(phaseIdx),
435 1210600 fluidState.pressure(phaseIdx)); //0.68 ;
436
437 else if (phaseIdx == gasPhaseIdx) // gas phase
438 3631798 return Component::gasThermalConductivity(fluidState.temperature(phaseIdx),
439 1210600 fluidState.pressure(phaseIdx)); //0.0248;
440
441 else
442 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index.");
443 }
444
445 using Base<Scalar, ThisType>::heatCapacity;
446 //! \copydoc Base<Scalar,ThisType>::heatCapacity(const FluidState&,int)
447 template <class FluidState>
448 2 static Scalar heatCapacity(const FluidState &fluidState,
449 int phaseIdx)
450 {
451
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 assert(0 <= phaseIdx && phaseIdx < numPhases);
452 // liquid phase
453
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if (phaseIdx == liquidPhaseIdx)
454 1 return Component::liquidHeatCapacity(fluidState.temperature(phaseIdx),
455 1 fluidState.pressure(phaseIdx));//4.217e3 ;
456
457 else if (phaseIdx == gasPhaseIdx) // gas phase
458 1 return Component::gasHeatCapacity(fluidState.temperature(phaseIdx),
459 1 fluidState.pressure(phaseIdx));//2.029e3;
460
461 else
462 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index.");
463 }
464 };
465
466 } // end namespace FluidSystems
467 } // end namespace Dumux
468
469 #endif
470