GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/material/fluidsystems/2pimmiscible.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 68 69 98.6%
Functions: 58 58 100.0%
Branches: 72 106 67.9%

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