GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/2p2c/evaporation/constant2p2cfluidsystem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 71 82 86.6%
Functions: 19 21 90.5%
Branches: 69 242 28.5%

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 TwoPTwoCTests
10 * \brief @copybrief Dumux::FluidSystems::ConstantTwoPTwoCFluidsystem
11 */
12 #ifndef DUMUX_CONSTANT_2P2C_FLUID_SYSTEM_HH
13 #define DUMUX_CONSTANT_2P2C_FLUID_SYSTEM_HH
14
15 #include <cassert>
16
17 #include <dumux/material/idealgas.hh>
18
19 #include <dumux/material/fluidsystems/base.hh>
20 #include <dumux/material/binarycoefficients/h2o_constant.hh>
21
22 #include <dumux/common/exceptions.hh>
23
24 namespace Dumux {
25 namespace FluidSystems {
26
27 /*!
28 * \ingroup TwoPTwoCTests
29 *
30 * \brief A two-phase fluid system with two components with constant properties that can be specified via the input file
31 */
32 template <class Scalar, class Component0, class Component1>
33 class ConstantTwoPTwoCFluidsystem
34 : public Base<Scalar, ConstantTwoPTwoCFluidsystem<Scalar, Component0, Component1> >
35 {
36 using ThisType = ConstantTwoPTwoCFluidsystem<Scalar, Component0, Component1>;
37 using Base = Dumux::FluidSystems::Base<Scalar, ThisType>;
38
39 public:
40
41 using ComponentOne = Component0;
42 using ComponentTwo = Component1;
43 //! Number of phases in the fluid system
44 static constexpr int numPhases = 2;
45
46 static constexpr int numComponents = 2; //!< Number of components in the fluid system
47
48 static constexpr int phase0Idx = 0; // index of the wetting phase
49 static constexpr int phase1Idx = 1; // index of the nonwetting phase
50
51 // export component indices to indicate the main component
52 static constexpr int comp0Idx = 0; // index of the wetting phase
53 static constexpr int comp1Idx = 1; // index of the nonwetting phase
54
55 /*!
56 * \brief Returns the human readable name of a fluid phase
57 *
58 * \param phaseIdx The index of the fluid phase to consider
59 */
60 72 static std::string phaseName(int phaseIdx)
61 {
62
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
72 assert(0 <= phaseIdx && phaseIdx < numPhases);
63
2/3
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
72 switch (phaseIdx)
64 {
65 72 case phase0Idx: return "liq";
66 72 case phase1Idx: return "gas";
67 }
68 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
69 }
70
71 /*!
72 * \brief Returns whether a phase is gaseous
73 *
74 * \param phaseIdx The index of the fluid phase to consider
75 */
76 static constexpr bool isGas(int phaseIdx)
77 {
78 assert(0 <= phaseIdx && phaseIdx < numPhases);
79 return phaseIdx == phase1Idx;
80 }
81
82 /*!
83 * \brief Returns true if and only if a fluid phase is assumed to
84 * be an ideal mixture.
85 *
86 * \param phaseIdx The index of the fluid phase to consider
87 */
88 static bool isIdealMixture(int phaseIdx)
89 {
90 assert(0 <= phaseIdx && phaseIdx < numPhases);
91 return true;
92 }
93
94 /*!
95 * \brief Returns true if and only if a fluid phase is assumed to
96 * be compressible.
97 *
98 * Compressible means that the partial derivative of the density
99 * to the fluid pressure is always larger than zero.
100 *
101 * \param phaseIdx The index of the fluid phase to consider
102 */
103 static constexpr bool isCompressible(int phaseIdx)
104 {
105 assert(0 <= phaseIdx && phaseIdx < numPhases);
106 // gases are always compressible
107 if (phaseIdx == phase1Idx)
108 return true;
109 // the liquid component decides for the liquid phase...
110 return ComponentOne::liquidIsCompressible();
111 }
112
113 /*!
114 * \brief Returns true if and only if a fluid phase is assumed to
115 * be an ideal gas.
116 *
117 * \param phaseIdx The index of the fluid phase to consider
118 */
119 static bool isIdealGas(int phaseIdx)
120 {
121 assert(0 <= phaseIdx && phaseIdx < numPhases);
122
123 if (phaseIdx == phase1Idx)
124 // let the components decide
125 return ComponentOne::gasIsIdeal() && ComponentTwo::gasIsIdeal();
126 return false; // not a gas
127 }
128
129 /*!
130 * \brief Returns the human readable name of a component
131 *
132 * \param compIdx The index of the component to consider
133 */
134 32 static std::string componentName(int compIdx)
135 {
136
4/9
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 28 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
36 static std::string name[] = {
137
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
2 ComponentOne::name(),
138
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
4 ComponentTwo::name()
139 };
140
141
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
32 assert(0 <= compIdx && compIdx < numComponents);
142 32 return name[compIdx];
143 }
144
145 /*!
146 * \brief Returns the molar mass of a component in \f$\mathrm{[kg/mol]}\f$.
147 *
148 * \param compIdx The index of the component to consider
149 */
150 17733200 static Scalar molarMass(int compIdx)
151 {
152
4/4
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 17733158 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 38 times.
17733204 static const Scalar M[] = {
153
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
2 ComponentOne::molarMass(),
154
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
4 ComponentTwo::molarMass(),
155 };
156
157
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17733200 times.
17733200 assert(0 <= compIdx && compIdx < numComponents);
158 17733200 return M[compIdx];
159 }
160
161 /*!
162 * \brief Critical temperature of a component \f$\mathrm{[K]}\f$.
163 *
164 * \param compIdx The index of the component to consider
165 */
166 static Scalar criticalTemperature(int compIdx)
167 {
168 static const Scalar Tcrit[] = {
169 ComponentOne::criticalTemperature(), // ComponentOne
170 ComponentTwo::criticalTemperature() // ComponentTwo
171 };
172
173 assert(0 <= compIdx && compIdx < numComponents);
174 return Tcrit[compIdx];
175 }
176
177 /*!
178 * \brief Critical pressure of a component \f$\mathrm{[Pa]}\f$.
179 *
180 * \param compIdx The index of the component to consider
181 */
182 static Scalar criticalPressure(int compIdx)
183 {
184 static const Scalar pcrit[] = {
185 ComponentOne::criticalPressure(),
186 ComponentTwo::criticalPressure()
187 };
188
189 assert(0 <= compIdx && compIdx < numComponents);
190 return pcrit[compIdx];
191 }
192
193
194 /****************************************
195 * thermodynamic relations
196 ****************************************/
197
198 /*!
199 * \brief Initializes the fluid system's static parameters generically
200 *
201 * If a tabulated ComponentOne component is used, we do our best to create
202 * tables that always work.
203 */
204 static void init()
205 {
206 4 init(/*tempMin=*/273.15,
207 /*tempMax=*/623.15,
208 /*numTemp=*/100,
209 /*pMin=*/0.0,
210 /*pMax=*/20e6,
211 /*numP=*/200);
212 }
213
214 /*!
215 * \brief Initializes the fluid system's static parameters using
216 * problem specific temperature and pressure ranges
217 *
218 * \param tempMin The minimum temperature used for tabulation of water \f$\mathrm{[K]}\f$
219 * \param tempMax The maximum temperature used for tabulation of water \f$\mathrm{[K]}\f$
220 * \param nTemp The number of ticks on the temperature axis of the table of water
221 * \param pressMin The minimum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$
222 * \param pressMax The maximum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$
223 * \param nPress The number of ticks on the pressure axis of the table of water
224 */
225 static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
226 Scalar pressMin, Scalar pressMax, unsigned nPress)
227 {
228
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
4 std::cout << "Using very simple constant component fluid system\n";
229 }
230
231 using Base::density;
232 /*!
233 * \brief Calculates the density \f$\mathrm{[kg/m^3]}\f$ of a fluid phase
234 *
235 * \param fluidState An arbitrary fluid state
236 * \param phaseIdx The index of the fluid phase to consider
237 */
238 template <class FluidState>
239 1433760 static Scalar density(const FluidState &fluidState,
240 int phaseIdx)
241 {
242
2/2
✓ Branch 0 taken 716880 times.
✓ Branch 1 taken 716880 times.
1433760 Scalar T = fluidState.temperature(phaseIdx);
243
2/2
✓ Branch 0 taken 716880 times.
✓ Branch 1 taken 716880 times.
1433760 Scalar p = fluidState.pressure(phaseIdx);
244
245 // liquid phase
246
2/2
✓ Branch 0 taken 716880 times.
✓ Branch 1 taken 716880 times.
1433760 if (phaseIdx == phase0Idx)
247 {
248 358440 return ComponentOne::liquidDensity(T, p);
249 }
250
1/2
✓ Branch 0 taken 716880 times.
✗ Branch 1 not taken.
716880 else if (phaseIdx == phase1Idx)// gas phase
251 {
252 716880 return ComponentTwo::gasDensity(T, p);
253 }
254 else DUNE_THROW(Dune::NotImplemented, "Wrong phase index");
255 }
256
257 using Base::molarDensity;
258 /*!
259 * \brief The molar density \f$\rho_{mol,\alpha}\f$
260 * of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
261 *
262 * This is a specific molar density for the combustion test
263 */
264 template <class FluidState>
265 716880 static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
266 {
267
2/2
✓ Branch 0 taken 358440 times.
✓ Branch 1 taken 358440 times.
716880 if (phaseIdx == phase0Idx)
268 {
269 358440 return density(fluidState, phaseIdx)/fluidState.averageMolarMass(phaseIdx);
270 }
271
1/2
✓ Branch 0 taken 358440 times.
✗ Branch 1 not taken.
358440 else if (phaseIdx == phase1Idx)
272 {
273 358440 return density(fluidState, phaseIdx)/fluidState.averageMolarMass(phaseIdx);
274 }
275 else DUNE_THROW(Dune::NotImplemented, "Wrong phase index");
276 }
277
278 using Base::viscosity;
279 /*!
280 * \brief Calculates the dynamic viscosity of a fluid phase \f$\mathrm{[Pa*s]}\f$
281 *
282 * \param fluidState An arbitrary fluid state
283 * \param phaseIdx The index of the fluid phase to consider
284 */
285 template <class FluidState>
286 716880 static Scalar viscosity(const FluidState &fluidState,
287 int phaseIdx)
288 {
289
2/2
✓ Branch 0 taken 358440 times.
✓ Branch 1 taken 358440 times.
716880 Scalar T = fluidState.temperature(phaseIdx);
290
2/2
✓ Branch 0 taken 358440 times.
✓ Branch 1 taken 358440 times.
716880 Scalar p = fluidState.pressure(phaseIdx);
291
292 // liquid phase
293
2/2
✓ Branch 0 taken 358440 times.
✓ Branch 1 taken 358440 times.
716880 if (phaseIdx == phase0Idx)
294 {
295 179220 return ComponentOne::liquidViscosity(T, p);
296 }
297
1/2
✓ Branch 0 taken 358440 times.
✗ Branch 1 not taken.
358440 else if (phaseIdx == phase1Idx) // gas phase
298 {
299 358440 return ComponentTwo::gasViscosity(T, p);
300 }
301 else DUNE_THROW(Dune::NotImplemented, "Wrong phase index");
302 }
303
304 using Base::fugacityCoefficient;
305 /*!
306 * \brief Calculates the fugacity coefficient \f$\mathrm{[-]}\f$ of an individual
307 * component in a fluid phase
308 *
309 * The fugacity coefficient \f$\phi^\kappa_\alpha\f$ of
310 * component \f$\kappa\f$ in phase \f$\alpha\f$ is connected to
311 * the fugacity \f$f^\kappa_\alpha\f$ and the component's mole
312 * fraction \f$x^\kappa_\alpha\f$ by means of the relation
313 *
314 * \f[
315 f^\kappa_\alpha = \phi^\kappa_\alpha\;x^\kappa_\alpha\;p_\alpha
316 \f]
317 * where \f$p_\alpha\f$ is the pressure of the fluid phase.
318 *
319 * The quantity "fugacity" itself is just an other way to express
320 * the chemical potential \f$\zeta^\kappa_\alpha\f$ of the
321 * component. It is defined via
322 *
323 * \f[
324 f^\kappa_\alpha := \exp\left\{\frac{\zeta^\kappa_\alpha}{k_B T_\alpha} \right\}
325 \f]
326 * where \f$k_B = 1.380\cdot10^{-23}\;J/K\f$ is the Boltzmann constant.
327 *
328 * \param fluidState An arbitrary fluid state
329 * \param phaseIdx The index of the fluid phase to consider
330 * \param compIdx The index of the component to consider
331 */
332 template <class FluidState>
333 1433760 static Scalar fugacityCoefficient(const FluidState &fluidState,
334 int phaseIdx,
335 int compIdx)
336 {
337
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1433760 times.
1433760 assert(0 <= phaseIdx && phaseIdx < numPhases);
338
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1433760 times.
1433760 assert(0 <= compIdx && compIdx < numComponents);
339
340
2/2
✓ Branch 0 taken 716880 times.
✓ Branch 1 taken 716880 times.
1433760 Scalar T = fluidState.temperature(phaseIdx);
341
2/2
✓ Branch 0 taken 716880 times.
✓ Branch 1 taken 716880 times.
1433760 Scalar p = fluidState.pressure(phaseIdx);
342 // liquid phase
343
2/2
✓ Branch 0 taken 716880 times.
✓ Branch 1 taken 716880 times.
1433760 if (phaseIdx == phase0Idx)
344 {
345
2/2
✓ Branch 0 taken 358440 times.
✓ Branch 1 taken 358440 times.
716880 if (compIdx == comp0Idx)
346 358440 return ComponentOne::vaporPressure(T)/p;
347 358440 return BinaryCoeff::H2O_Component<Scalar, ComponentTwo>::henryCompInWater(T)/p;
348 }
349
350 // for the gas phase, assume an ideal gas when it comes to
351 // fugacity (-> fugacity == partial pressure)
352 return 1.0;
353 }
354
355 using Base::diffusionCoefficient;
356 /*!
357 * \brief Calculates the molecular diffusion coefficient for a
358 * component in a fluid phase \f$\mathrm{[mol^2 * s / (kg*m^3)]}\f$
359 *
360 * \param fluidState An arbitrary fluid state
361 * \param phaseIdx The index of the fluid phase to consider
362 * \param compIdx The index of the component to consider
363 */
364 template <class FluidState>
365 static Scalar diffusionCoefficient(const FluidState &fluidState,
366 int phaseIdx,
367 int compIdx)
368 {
369 DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients");
370 }
371
372 using Base::binaryDiffusionCoefficient;
373 /*!
374 * \brief Given a phase's composition, temperature and pressure,
375 * returns the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ for components
376 * \f$i\f$ and \f$j\f$ in this phase.
377 *
378 * \param fluidState An arbitrary fluid state
379 * \param phaseIdx The index of the fluid phase to consider
380 * \param compIIdx The index of the first component to consider
381 * \param compJIdx The index of the second component to consider
382 */
383 template <class FluidState>
384 740880 static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
385 int phaseIdx,
386 int compIIdx,
387 int compJIdx)
388
389 {
390 using std::swap;
391
2/2
✓ Branch 0 taken 382440 times.
✓ Branch 1 taken 358440 times.
740880 if (compIIdx > compJIdx)
392 382440 swap(compIIdx, compJIdx);
393
394 // we are in the liquid phase
395
2/2
✓ Branch 0 taken 358440 times.
✓ Branch 1 taken 382440 times.
740880 if (phaseIdx == phase0Idx)
396 {
397
2/4
✓ Branch 0 taken 358440 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 358440 times.
358440 if (compIIdx == comp0Idx && compJIdx == comp1Idx)
398 return 1e-9;
399 else
400 DUNE_THROW(Dune::InvalidStateException,
401 "Binary diffusion coefficient of components "
402 << compIIdx << " and " << compJIdx
403 << " in phase " << phaseIdx << " is undefined!\n");
404 }
405
406 // we are in the gas phase
407
1/2
✓ Branch 0 taken 382440 times.
✗ Branch 1 not taken.
382440 else if (phaseIdx == phase1Idx)
408 {
409
2/4
✓ Branch 0 taken 382440 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 382440 times.
382440 if (compIIdx == comp0Idx && compJIdx == comp1Idx)
410 return 1e-5;
411 else
412 DUNE_THROW(Dune::InvalidStateException,
413 "Binary diffusion coefficient of components "
414 << compIIdx << " and " << compJIdx
415 << " in phase " << phaseIdx << " is undefined!\n");
416 }
417
418 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
419 }
420
421 using Base::enthalpy;
422 /*!
423 * \brief Calculates specific enthalpy \f$\mathrm{[J/kg]}\f$.
424 *
425 * \param fluidState An arbitrary fluid state
426 * \param phaseIdx The index of the fluid phase to consider
427 */
428 template <class FluidState>
429 740880 static Scalar enthalpy(const FluidState &fluidState,
430 int phaseIdx)
431 {
432
2/2
✓ Branch 0 taken 358440 times.
✓ Branch 1 taken 382440 times.
740880 const Scalar T = fluidState.temperature(phaseIdx);
433
2/2
✓ Branch 0 taken 358440 times.
✓ Branch 1 taken 382440 times.
740880 const Scalar p = fluidState.pressure(phaseIdx);
434
435 // liquid phase
436
2/2
✓ Branch 0 taken 358440 times.
✓ Branch 1 taken 382440 times.
740880 if (phaseIdx == phase0Idx)
437 {
438 358440 return Component0::liquidEnthalpy(T, p);
439 }
440
1/2
✓ Branch 0 taken 382440 times.
✗ Branch 1 not taken.
382440 else if (phaseIdx == phase1Idx) // gas phase
441 {
442 // assume ideal mixture: which means
443 // that the total specific enthalpy is the sum of the
444 // "partial specific enthalpies" of the components.
445 382440 Scalar hComponentOne =
446 fluidState.massFraction(phase1Idx, comp0Idx)
447 382440 * ComponentOne::gasEnthalpy(T, p);
448 382440 Scalar hComponentTwo =
449 fluidState.massFraction(phase1Idx, comp1Idx)
450 382440 * ComponentTwo::gasEnthalpy(T, p);
451 382440 return hComponentOne + hComponentTwo;
452 }
453 else DUNE_THROW(Dune::NotImplemented, "Wrong phase index");
454 }
455
456 using Base::thermalConductivity;
457 /*!
458 * \brief Thermal conductivity of a fluid phase \f$\mathrm{[W/(m K)]}\f$.
459 *
460 * Use the conductivity of vapor and liquid water at 100°C
461 *
462 * \param fluidState An arbitrary fluid state
463 * \param phaseIdx The index of the fluid phase to consider
464 */
465 template <class FluidState>
466 358440 static Scalar thermalConductivity(const FluidState &fluidState,
467 const int phaseIdx)
468 {
469
2/2
✓ Branch 1 taken 179220 times.
✓ Branch 2 taken 179220 times.
740880 const Scalar temperature = fluidState.temperature(phaseIdx) ;
470
2/2
✓ Branch 1 taken 179220 times.
✓ Branch 2 taken 179220 times.
920100 const Scalar pressure = fluidState.pressure(phaseIdx);
471
2/2
✓ Branch 0 taken 179220 times.
✓ Branch 1 taken 179220 times.
358440 if (phaseIdx == phase0Idx)
472 179220 return ComponentOne::liquidThermalConductivity(temperature, pressure);
473 else
474 382440 return ComponentTwo::gasThermalConductivity(temperature, pressure);
475 }
476
477 using Base::heatCapacity;
478 /*!
479 * \brief Specific isobaric heat capacity of a fluid phase.
480 * \f$\mathrm{[J/kg / K]}\f$.
481 *
482 * \param fluidState An arbitrary fluid state
483 * \param phaseIdx The index of the fluid phase to consider
484 */
485 template <class FluidState>
486 static Scalar heatCapacity(const FluidState &fluidState,
487 int phaseIdx)
488 {
489 assert(0 <= phaseIdx && phaseIdx < numPhases);
490 // liquid phase
491 const Scalar temperature = fluidState.temperature(phaseIdx);
492 const Scalar pressure = fluidState.pressure(phaseIdx);
493 if (phaseIdx == phase0Idx)
494 {
495 return ComponentOne::liquidHeatCapacity(temperature, pressure);
496 }
497 else if (phaseIdx == phase1Idx) // gas phase
498 {
499 return ComponentOne::gasHeatCapacity(temperature, pressure) * fluidState.moleFraction(phase1Idx, comp0Idx)
500 + ComponentTwo::gasHeatCapacity(temperature, pressure) * fluidState.moleFraction(phase1Idx, comp1Idx);
501 }
502 else DUNE_THROW(Dune::NotImplemented, "Wrong phase index");
503 }
504 };
505
506 } // end namespace FluidSystems
507 } // end namespace Dumux
508
509 #endif
510