GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/material/fluidsystems/2pimmiscible.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 68 73 93.2%
Functions: 55 65 84.6%
Branches: 81 128 63.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-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::TwoPImmiscible
11 */
12 #ifndef DUMUX_2P_IMMISCIBLE_FLUID_SYSTEM_HH
13 #define DUMUX_2P_IMMISCIBLE_FLUID_SYSTEM_HH
14
15 #include <limits>
16 #include <cassert>
17
18 #include <dune/common/exceptions.hh>
19
20 #include <dumux/material/fluidsystems/1pliquid.hh>
21 #include <dumux/material/fluidsystems/1pgas.hh>
22 #include <dumux/material/fluidstates/immiscible.hh>
23 #include <dumux/material/components/base.hh>
24 #include <dumux/io/name.hh>
25
26 #include "base.hh"
27
28 namespace Dumux {
29 namespace FluidSystems {
30
31 /*!
32 * \ingroup FluidSystems
33 * \brief A fluid system for two-phase models assuming immiscibility and
34 * thermodynamic equilibrium
35 *
36 * The fluid phases are completely specified by means of their
37 * constituting components.
38 * The fluids can be defined individually via FluidSystem::OnePLiquid<Scalar, Component> and
39 * FluidSystem::OnePGas<Scalar, Component>. These fluids consist of one component.
40 * \tparam Scalar the scalar type
41 * \tparam Fluid0 a one-phase fluid system (use FluidSystem::OnePLiquid<Scalar, Component> / FluidSystem::OnePGas<Scalar, Component>)
42 * \tparam Fluid1 a one-phase fluid system (use FluidSystem::OnePLiquid<Scalar, Component> / FluidSystem::OnePGas<Scalar, Component>)
43 */
44 template <class Scalar, class Fluid0, class Fluid1>
45 class TwoPImmiscible
46 : public Base<Scalar, TwoPImmiscible<Scalar, Fluid0, Fluid1> >
47 {
48 static_assert((Fluid0::numPhases == 1), "Fluid0 has more than one phase");
49 static_assert((Fluid1::numPhases == 1), "Fluid1 has more than one phase");
50 static_assert((Fluid0::numComponents == 1), "Fluid0 has more than one component");
51 static_assert((Fluid1::numComponents == 1), "Fluid1 has more than one component");
52 // two gaseous phases at once do not make sense physically! (but two liquids are fine)
53 static_assert(!Fluid0::isGas() || !Fluid1::isGas(), "One phase has to be a liquid!");
54
55 using ThisType = TwoPImmiscible<Scalar, Fluid0, Fluid1>;
56
57 public:
58 static constexpr int numPhases = 2; //!< Number of phases in the fluid system
59 static constexpr int numComponents = 2; //!< Number of components in the fluid system
60
61 static constexpr int phase0Idx = 0; //!< index of the first phase
62 static constexpr int phase1Idx = 1; //!< index of the second phase
63 static constexpr int comp0Idx = 0; //!< index of the first component
64 static constexpr int comp1Idx = 1; //!< index of the second component
65
66 /****************************************
67 * Fluid phase related static parameters
68 ****************************************/
69 /*!
70 * \brief Return the human readable name of a fluid phase
71 * \param phaseIdx The index of the fluid phase to consider
72 */
73 322 static std::string phaseName(int phaseIdx)
74 {
75
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 314 times.
322 assert(0 <= phaseIdx && phaseIdx < numPhases);
76
77 if (!Fluid0::isGas() && !Fluid1::isGas())
78 {
79
3/4
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 176 times.
✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
204 static const auto name0 = Components::IsAqueous<typename Fluid0::Component>::value ? IOName::aqueousPhase() : IOName::naplPhase();
80
3/4
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 176 times.
✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
204 static const auto name1 = Components::IsAqueous<typename Fluid1::Component>::value ? IOName::aqueousPhase() : IOName::naplPhase();
81
82
2/2
✓ Branch 0 taken 200 times.
✓ Branch 1 taken 2 times.
204 if (name0 != name1)
83
2/2
✓ Branch 0 taken 100 times.
✓ Branch 1 taken 100 times.
300 return (phaseIdx == phase0Idx) ? name0 : name1;
84 else
85
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
4 return (phaseIdx == phase0Idx) ? name0 + "_0" : name1 + "_1";
86 }
87 else
88 {
89
2/2
✓ Branch 0 taken 45 times.
✓ Branch 1 taken 67 times.
118 if (phaseIdx == phase0Idx)
90 48 return Fluid0::isGas() ? IOName::gaseousPhase() : IOName::liquidPhase();
91 else
92 70 return Fluid1::isGas() ? IOName::gaseousPhase() : IOName::liquidPhase();
93 }
94 }
95
96 /*!
97 * \brief Returns whether the fluids are miscible
98 */
99 static constexpr bool isMiscible()
100 { return false; }
101
102 /*!
103 * \brief Return whether a phase is gaseous
104 * \param phaseIdx The index of the fluid phase to consider
105 */
106 static constexpr bool isGas(int phaseIdx)
107 {
108 assert(0 <= phaseIdx && phaseIdx < numPhases);
109
110 if (phaseIdx == phase0Idx)
111 return Fluid0::isGas();
112 return Fluid1::isGas();
113 }
114
115 /*!
116 * \brief Returns true if and only if a fluid phase is assumed to
117 * be an ideal mixture.
118 * \param phaseIdx The index of the fluid phase to consider
119 *
120 * We define an ideal mixture as a fluid phase where the fugacity
121 * coefficients of all components times the pressure of the phase
122 * are independent on the fluid composition. This assumption is true
123 * if immiscibility is assumed. If you are unsure what
124 * this function should return, it is safe to return false. The
125 * only damage done will be (slightly) increased computation times
126 * in some cases.
127 */
128 static bool isIdealMixture(int phaseIdx)
129 {
130 assert(0 <= phaseIdx && phaseIdx < numPhases);
131
132 // we assume immisibility
133 return true;
134 }
135
136 /*!
137 * \brief Returns true if and only if a fluid phase is assumed to
138 * be an ideal gas.
139 *
140 * \param phaseIdx The index of the fluid phase to consider
141 */
142 static constexpr bool isIdealGas(int phaseIdx)
143 {
144 assert(0 <= phaseIdx && phaseIdx < numPhases);
145
146 // let the fluids decide
147 if (phaseIdx == phase0Idx)
148 return Fluid0::isIdealGas();
149 return Fluid1::isIdealGas();
150 }
151
152 /*!
153 * \brief Returns true if and only if a fluid phase is assumed to
154 * be compressible.
155 *
156 * Compressible means. that the partial derivative of the density
157 * to the fluid pressure is always larger than zero.
158 *
159 * \param phaseIdx The index of the fluid phase to consider
160 */
161 static constexpr bool isCompressible(int phaseIdx)
162 {
163 assert(0 <= phaseIdx && phaseIdx < numPhases);
164
165 // let the fluids decide
166 if (phaseIdx == phase0Idx)
167 return Fluid0::isCompressible();
168 return Fluid1::isCompressible();
169 }
170
171 /*!
172 * \brief Returns true if the liquid phase viscostiy is constant
173 *
174 * \param phaseIdx The index of the fluid phase to consider
175 */
176 static constexpr bool viscosityIsConstant(int phaseIdx)
177 {
178 assert(0 <= phaseIdx && phaseIdx < numPhases);
179
180 // let the fluids decide
181 if (phaseIdx == phase0Idx)
182 return Fluid0::viscosityIsConstant();
183 return Fluid1::viscosityIsConstant();
184 }
185
186 /*!
187 * \brief Returns true if and only if a fluid phase is assumed to
188 * be an ideal gas.
189 *
190 * \param phaseIdx The index of the fluid phase to consider
191 */
192 static bool isIdealFluid1(int phaseIdx)
193 {
194 assert(0 <= phaseIdx && phaseIdx < numPhases);
195
196 // let the fluids decide
197 if (phaseIdx == phase0Idx)
198 return Fluid0::isIdealFluid1();
199 return Fluid1::isIdealFluid1();
200 }
201
202 /****************************************
203 * Component related static parameters
204 ****************************************/
205 /*!
206 * \brief Return the human readable name of a component
207 *
208 * \param compIdx index of the component
209 */
210 16 static std::string componentName(int compIdx)
211 {
212
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
16 assert(0 <= compIdx && compIdx < numComponents);
213
214
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
16 if (compIdx == comp0Idx)
215 8 return Fluid0::name();
216 8 return Fluid1::name();
217 }
218
219 /*!
220 * \brief Return the molar mass of a component in \f$\mathrm{[kg/mol]}\f$.
221 * \param compIdx index of the component
222 */
223 static Scalar molarMass(int compIdx)
224 {
225
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
4 assert(0 <= compIdx && compIdx < numComponents);
226
227
6/8
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 4 times.
20 if (compIdx == comp0Idx)
228 return Fluid0::molarMass();
229 10 return Fluid1::molarMass();
230 }
231
232 /*!
233 * \brief Critical temperature of a component \f$\mathrm{[K]}\f$.
234 * \param compIdx index of the component
235 */
236 static Scalar criticalTemperature(int compIdx)
237 {
238 assert(0 <= compIdx && compIdx < numComponents);
239
240 if (compIdx == comp0Idx)
241 return Fluid0::criticalTemperature();
242 return Fluid1::criticalTemperature();
243 }
244
245 /*!
246 * \brief Critical pressure of a component \f$\mathrm{[Pa]}\f$.
247 * \param compIdx index of the component
248 */
249 static Scalar criticalPressure(int compIdx)
250 {
251 assert(0 <= compIdx && compIdx < numComponents);
252
253 if (compIdx == comp0Idx)
254 return Fluid0::criticalPressure();
255 return Fluid1::criticalPressure();
256 }
257
258 /*!
259 * \brief The acentric factor of a component \f$\mathrm{[-]}\f$.
260 * \param compIdx index of the component
261 */
262 static Scalar acentricFactor(int compIdx)
263 {
264 assert(0 <= compIdx && compIdx < numComponents);
265
266 if (compIdx == comp0Idx)
267 return Fluid0::acentricFactor();
268 return Fluid1::acentricFactor();
269 }
270
271 /****************************************
272 * thermodynamic relations
273 ****************************************/
274
275 /*!
276 * \brief Initialize the fluid system's static parameters
277 */
278 static void init()
279 {
280 // initialize with some default values
281 5 init(/*tempMin=*/273.15, /*tempMax=*/623.15, /*numTemp=*/100,
282 /*pMin=*/-10.0, /*pMax=*/20e6, /*numP=*/200);
283 }
284
285 /*!
286 * \brief Initialize the fluid system's static parameters using
287 * problem specific temperature and pressure ranges
288 *
289 * \param tempMin The minimum temperature used for tabulation of water \f$\mathrm{[K]}\f$
290 * \param tempMax The maximum temperature used for tabulation of water\f$\mathrm{[K]}\f$
291 * \param nTemp The number of ticks on the temperature axis of the table of water
292 * \param pressMin The minimum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$
293 * \param pressMax The maximum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$
294 * \param nPress The number of ticks on the pressure axis of the table of water
295 */
296 1 static void init(Scalar tempMin, Scalar tempMax, std::size_t nTemp,
297 Scalar pressMin, Scalar pressMax, std::size_t nPress)
298 {
299 if (Fluid0::Component::isTabulated)
300 1 Fluid0::Component::init(tempMin, tempMax, nTemp, pressMin, pressMax, nPress);
301
302 if (Fluid1::Component::isTabulated)
303
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
2 Fluid1::Component::init(tempMin, tempMax, nTemp, pressMin, pressMax, nPress);
304 1 }
305
306 using Base<Scalar, ThisType>::density;
307 //! \copybrief Base<Scalar,ThisType>::density(const FluidState&,int)
308 template <class FluidState>
309 24476900 static Scalar density(const FluidState &fluidState,
310 int phaseIdx)
311 {
312
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 23281972 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
24476900 assert(0 <= phaseIdx && phaseIdx < numPhases);
313
314
6/6
✓ Branch 0 taken 11926082 times.
✓ Branch 1 taken 11355882 times.
✓ Branch 2 taken 8251166 times.
✓ Branch 3 taken 8251166 times.
✓ Branch 4 taken 57044720 times.
✓ Branch 5 taken 57044720 times.
155068672 Scalar temperature = fluidState.temperature(phaseIdx);
315
6/6
✓ Branch 0 taken 11926082 times.
✓ Branch 1 taken 11355882 times.
✓ Branch 2 taken 8251166 times.
✓ Branch 3 taken 8251166 times.
✓ Branch 4 taken 57044720 times.
✓ Branch 5 taken 57044720 times.
155068672 Scalar pressure = fluidState.pressure(phaseIdx);
316
6/6
✓ Branch 0 taken 11926086 times.
✓ Branch 1 taken 11355886 times.
✓ Branch 2 taken 8251166 times.
✓ Branch 3 taken 8251166 times.
✓ Branch 4 taken 57044720 times.
✓ Branch 5 taken 57044720 times.
155068672 if (phaseIdx == phase0Idx)
317 6999748 return Fluid0::density(temperature, pressure);
318 90913186 return Fluid1::density(temperature, pressure);
319 }
320
321 using Base<Scalar, ThisType>::molarDensity;
322 //! \copybrief Base<Scalar,ThisType>::molarDensity(const FluidState&,int)
323 template <class FluidState>
324 16 static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
325 {
326 16 Scalar temperature = fluidState.temperature(phaseIdx);
327 16 Scalar pressure = fluidState.pressure(phaseIdx);
328
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
16 if (phaseIdx == phase0Idx)
329 12 return Fluid0::molarDensity(temperature, pressure);
330 10 return Fluid1::molarDensity(temperature, pressure);
331 }
332
333 using Base<Scalar, ThisType>::viscosity;
334 //! \copybrief Base<Scalar,ThisType>::viscosity(const FluidState&,int)
335 template <class FluidState>
336 16107836 static Scalar viscosity(const FluidState &fluidState,
337 int phaseIdx)
338 {
339
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 14912908 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
16107836 assert(0 <= phaseIdx && phaseIdx < numPhases);
340
341
4/6
✓ Branch 0 taken 3869378 times.
✓ Branch 1 taken 11043522 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 65295886 times.
✓ Branch 5 taken 65295886 times.
146699608 Scalar temperature = fluidState.temperature(phaseIdx);
342
4/6
✓ Branch 0 taken 3869378 times.
✓ Branch 1 taken 11043522 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 65295886 times.
✓ Branch 5 taken 65295886 times.
146699608 Scalar pressure = fluidState.pressure(phaseIdx);
343
4/6
✓ Branch 0 taken 3869382 times.
✓ Branch 1 taken 11043526 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 65295886 times.
✓ Branch 5 taken 65295886 times.
146699608 if (phaseIdx == phase0Idx)
344 6999748 return Fluid0::viscosity(temperature, pressure);
345 74799778 return Fluid1::viscosity(temperature, pressure);
346 }
347
348 using Base<Scalar, ThisType>::fugacityCoefficient;
349 //! \copybrief Base<Scalar,ThisType>::fugacityCoefficient(const FluidState&,int,int)
350 template <class FluidState>
351 static Scalar fugacityCoefficient(const FluidState &fluidState,
352 int phaseIdx,
353 int compIdx)
354 {
355 assert(0 <= phaseIdx && phaseIdx < numPhases);
356 assert(0 <= compIdx && compIdx < numComponents);
357
358 if (phaseIdx == compIdx)
359 // We could calculate the real fugacity coefficient of
360 // the component in the fluid. Probably that's not worth
361 // the effort, since the fugacity coefficient of the other
362 // component is infinite anyway...
363 return 1.0;
364 return std::numeric_limits<Scalar>::infinity();
365 }
366
367 using Base<Scalar, ThisType>::diffusionCoefficient;
368 //! \copybrief Base<Scalar,ThisType>::diffusionCoefficient(const FluidState&,int,int)
369 template <class FluidState>
370 28 static Scalar diffusionCoefficient(const FluidState &fluidState,
371 int phaseIdx,
372 int compIdx)
373 {
374
7/16
✓ Branch 2 taken 14 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 14 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 14 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 14 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 14 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 14 times.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 14 times.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
308 DUNE_THROW(Dune::InvalidStateException,
375 "Diffusion coefficients of components are meaningless if"
376 " immiscibility is assumed");
377 }
378
379 using Base<Scalar, ThisType>::binaryDiffusionCoefficient;
380 //! \copybrief Base<Scalar,ThisType>::binaryDiffusionCoefficient(const FluidState&,int,int,int)
381 template <class FluidState>
382 52 static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
383 int phaseIdx,
384 int compIIdx,
385 int compJIdx)
386
387 {
388
7/16
✓ Branch 2 taken 26 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 26 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 26 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 26 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 26 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 26 times.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 26 times.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
572 DUNE_THROW(Dune::InvalidStateException,
389 "Binary diffusion coefficients of components are meaningless if"
390 " immiscibility is assumed");
391 }
392
393 using Base<Scalar, ThisType>::enthalpy;
394 //! \copybrief Base<Scalar,ThisType>::enthalpy(const FluidState&,int)
395 template <class FluidState>
396 16 static Scalar enthalpy(const FluidState &fluidState,
397 int phaseIdx)
398 {
399
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
16 assert(0 <= phaseIdx && phaseIdx < numPhases);
400
401 16 Scalar temperature = fluidState.temperature(phaseIdx);
402 16 Scalar pressure = fluidState.pressure(phaseIdx);
403
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
16 if (phaseIdx == phase0Idx)
404 12 return Fluid0::enthalpy(temperature, pressure);
405 14 return Fluid1::enthalpy(temperature, pressure);
406 }
407
408 using Base<Scalar, ThisType>::thermalConductivity;
409 //! \copybrief Base<Scalar,ThisType>::thermalConductivity(const FluidState&,int)
410 template <class FluidState>
411 16 static Scalar thermalConductivity(const FluidState &fluidState,
412 int phaseIdx)
413 {
414
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
16 assert(0 <= phaseIdx && phaseIdx < numPhases);
415
416 16 Scalar temperature = fluidState.temperature(phaseIdx);
417 16 Scalar pressure = fluidState.pressure(phaseIdx);
418
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
16 if (phaseIdx == phase0Idx)
419 16 return Fluid0::thermalConductivity(temperature, pressure);
420 16 return Fluid1::thermalConductivity(temperature, pressure);
421 }
422
423 using Base<Scalar, ThisType>::heatCapacity;
424 //! \copybrief Base<Scalar,ThisType>::heatCapacity(const FluidState&,int)
425 template <class FluidState>
426 8 static Scalar heatCapacity(const FluidState &fluidState,
427 int phaseIdx)
428 {
429 8 assert(0 <= phaseIdx && phaseIdx < numPhases);
430
431 8 Scalar temperature = fluidState.temperature(phaseIdx);
432 8 Scalar pressure = fluidState.pressure(phaseIdx);
433 8 if (phaseIdx == phase0Idx)
434 6 return Fluid0::heatCapacity(temperature, pressure);
435 7 return Fluid1::heatCapacity(temperature, pressure);
436 }
437 };
438
439 } // end namespace FluidSystems
440 } // end namespace Dumux
441
442 #endif
443