GCC Code Coverage Report


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