GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/mpnc/thermalnonequilibrium/combustionfluidsystem.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 0 44 0.0%
Functions: 0 5 0.0%
Branches: 0 119 0.0%
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 MPNCTests
10 * \brief @copybrief Dumux::FluidSystems::CombustionFluidsystem
11 */
12 #ifndef DUMUX_PURE_WATER_FLUID_SYSTEM_HH
13 #define DUMUX_PURE_WATER_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/components/n2.hh>
21 #include <dumux/material/components/simpleh2o.hh>
22 #include <dumux/material/binarycoefficients/h2o_n2.hh>
23
24 #include <dumux/common/exceptions.hh>
25
26 namespace Dumux {
27 namespace FluidSystems {
28
29 /*!
30 * \ingroup MPNCTests
31 *
32 * \brief A two-phase fluid system with water as sole component.
33 *
34 * Values are taken from Shi and Wang, 2011 \cite Shi2011.
35 */
36 template <class Scalar>
37 class CombustionFluidsystem
38 : public Base<Scalar, CombustionFluidsystem<Scalar> >
39 {
40 using ThisType = CombustionFluidsystem<Scalar>;
41 using Base = Dumux::FluidSystems::Base<Scalar, ThisType>;
42
43 // convenience using declarations
44 using IdealGas = Dumux::IdealGas<Scalar>;
45 using SimpleH2O = Dumux::Components::SimpleH2O<Scalar>;
46 using SimpleN2 = Dumux::Components::N2<Scalar>;
47
48 public:
49 /****************************************
50 * Fluid phase related static parameters
51 ****************************************/
52
53 //! Number of phases in the fluid system
54 static constexpr int numPhases = 2;
55
56 static constexpr int wPhaseIdx = 0; // index of the wetting phase
57 static constexpr int nPhaseIdx = 1; // index of the nonwetting phase
58 static constexpr int phase0Idx = 0; // index of the wetting phase
59 static constexpr int phase1Idx = 1; // index of the nonwetting phase
60
61 // export component indices to indicate the main component
62 // of the corresponding phase at atmospheric pressure 1 bar
63 // and room temperature 20°C:
64 static constexpr int wCompIdx = wPhaseIdx;
65 static constexpr int nCompIdx = nPhaseIdx;
66 static constexpr int comp0Idx = 0; // index of the wetting phase
67 static constexpr int comp1Idx = 1; // index of the nonwetting phase
68
69 /*!
70 * \brief Returns the human readable name of a fluid phase
71 *
72 * \param phaseIdx The index of the fluid phase to consider
73 */
74 ✗ static std::string phaseName(int phaseIdx)
75 {
76 ✗ assert(0 <= phaseIdx && phaseIdx < numPhases);
77 ✗ switch (phaseIdx)
78 {
79 ✗ case wPhaseIdx: return "liq";
80 ✗ case nPhaseIdx: return "gas";
81 }
82 ✗ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
83 }
84
85 /*!
86 * \brief Returns whether a phase is gaseous
87 *
88 * \param phaseIdx The index of the fluid phase to consider
89 */
90 static constexpr bool isGas(int phaseIdx)
91 {
92 assert(0 <= phaseIdx && phaseIdx < numPhases);
93 return phaseIdx == nPhaseIdx;
94 }
95
96 /*!
97 * \brief Returns true if and only if a fluid phase is assumed to
98 * be an ideal mixture.
99 *
100 * We define an ideal mixture as a fluid phase where the fugacity
101 * coefficients of all components times the pressure of the phase
102 * are independent on the fluid composition. This assumption is true
103 * if Henry's law and Raoult's law apply. If you are unsure what
104 * this function should return, it is safe to return false. The
105 * only damage done will be (slightly) increased computation times
106 * in some cases.
107 *
108 * \param phaseIdx The index of the fluid phase to consider
109 */
110 static bool isIdealMixture(int phaseIdx)
111 {
112 ✗ assert(0 <= phaseIdx && phaseIdx < numPhases);
113 // we assume Henry's and Raoult's laws for the water phase and
114 // and no interaction between gas molecules of different
115 // components, so all phases are ideal mixtures!
116 return true;
117 }
118
119 /*!
120 * \brief Returns true if and only if a fluid phase is assumed to
121 * be compressible.
122 *
123 * Compressible means that the partial derivative of the density
124 * to the fluid pressure is always larger than zero.
125 *
126 * \param phaseIdx The index of the fluid phase to consider
127 */
128 static constexpr bool isCompressible(int phaseIdx)
129 {
130 assert(0 <= phaseIdx && phaseIdx < numPhases);
131 // gases are always compressible
132 if (phaseIdx == nPhaseIdx)
133 return true;
134 // the water component decides for the liquid phase...
135 return H2O::liquidIsCompressible();
136 }
137
138 /*!
139 * \brief Returns true if and only if a fluid phase is assumed to
140 * be an ideal gas.
141 *
142 * \param phaseIdx The index of the fluid phase to consider
143 */
144 static bool isIdealGas(int phaseIdx)
145 {
146 assert(0 <= phaseIdx && phaseIdx < numPhases);
147
148 if (phaseIdx == nPhaseIdx)
149 // let the components decide
150 return H2O::gasIsIdeal() && N2::gasIsIdeal();
151 return false; // not a gas
152 }
153
154 /****************************************
155 * Component related static parameters
156 ****************************************/
157
158 //! Number of components in the fluid system
159 static constexpr int numComponents = 2;
160
161 static constexpr int H2OIdx = wCompIdx;
162 static constexpr int N2Idx = nCompIdx;
163
164 //! The components for pure water
165 using H2O = SimpleH2O;
166
167 //! The components for pure nitrogen
168 using N2 = SimpleN2;
169
170 /*!
171 * \brief Returns the human readable name of a component
172 *
173 * \param compIdx The index of the component to consider
174 */
175 ✗ static std::string componentName(int compIdx)
176 {
177 ✗ static std::string name[] = {
178 H2O::name(),
179 N2::name()
180 };
181
182 ✗ assert(0 <= compIdx && compIdx < numComponents);
183 ✗ return name[compIdx];
184 }
185
186 /*!
187 * \brief Returns the molar mass of a component in \f$\mathrm{[kg/mol]}\f$.
188 *
189 * \param compIdx The index of the component to consider
190 */
191 static Scalar molarMass(int compIdx)
192 {
193 static const Scalar M[] = {
194 H2O::molarMass(),
195 N2::molarMass(),
196 };
197
198 assert(0 <= compIdx && compIdx < numComponents);
199 ✗ return M[compIdx];
200 }
201
202 /*!
203 * \brief Critical temperature of a component \f$\mathrm{[K]}\f$.
204 *
205 * \param compIdx The index of the component to consider
206 */
207 static Scalar criticalTemperature(int compIdx)
208 {
209 static const Scalar Tcrit[] = {
210 H2O::criticalTemperature(), // H2O
211 N2::criticalTemperature() // N2
212 };
213
214 assert(0 <= compIdx && compIdx < numComponents);
215 return Tcrit[compIdx];
216 }
217
218 /*!
219 * \brief Critical pressure of a component \f$\mathrm{[Pa]}\f$.
220 *
221 * \param compIdx The index of the component to consider
222 */
223 static Scalar criticalPressure(int compIdx)
224 {
225 static const Scalar pcrit[] = {
226 H2O::criticalPressure(),
227 N2::criticalPressure()
228 };
229
230 assert(0 <= compIdx && compIdx < numComponents);
231 return pcrit[compIdx];
232 }
233
234 /*!
235 * \brief Molar volume of a component at the critical point \f$\mathrm{[m^3/mol]}\f$.
236 *
237 * \param compIdx The index of the component to consider
238 */
239 static Scalar criticalMolarVolume(int compIdx)
240 {
241 DUNE_THROW(Dune::NotImplemented, "criticalMolarVolume()");
242 }
243
244 /*!
245 * \brief The acentric factor of a component \f$\mathrm{[-]}\f$.
246 *
247 * \param compIdx The index of the component to consider
248 */
249 static Scalar acentricFactor(int compIdx)
250 {
251 static const Scalar accFac[] = {
252 H2O::acentricFactor(), // H2O (from Reid, et al.)
253 N2::acentricFactor()
254 };
255
256 assert(0 <= compIdx && compIdx < numComponents);
257 return accFac[compIdx];
258 }
259
260 /****************************************
261 * thermodynamic relations
262 ****************************************/
263
264 /*!
265 * \brief Initializes the fluid system's static parameters generically
266 *
267 * If a tabulated H2O component is used, we do our best to create
268 * tables that always work.
269 */
270 static void init()
271 {
272 init(/*tempMin=*/273.15,
273 /*tempMax=*/623.15,
274 /*numTemp=*/100,
275 /*pMin=*/0.0,
276 /*pMax=*/20e6,
277 /*numP=*/200);
278 }
279
280 /*!
281 * \brief Initializes the fluid system's static parameters using
282 * problem specific temperature and pressure ranges
283 *
284 * \param tempMin The minimum temperature used for tabulation of water \f$\mathrm{[K]}\f$
285 * \param tempMax The maximum temperature used for tabulation of water \f$\mathrm{[K]}\f$
286 * \param nTemp The number of ticks on the temperature axis of the table of water
287 * \param pressMin The minimum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$
288 * \param pressMax The maximum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$
289 * \param nPress The number of ticks on the pressure axis of the table of water
290 */
291 static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
292 Scalar pressMin, Scalar pressMax, unsigned nPress)
293 {
294 std::cout << "Using very simple pure water fluid system\n";
295 }
296
297 using Base::density;
298 /*!
299 * \brief Calculates the density \f$\mathrm{[kg/m^3]}\f$ of a fluid phase
300 *
301 * \param fluidState An arbitrary fluid state
302 * \param phaseIdx The index of the fluid phase to consider
303 */
304 template <class FluidState>
305 static Scalar density(const FluidState &fluidState,
306 int phaseIdx)
307 {
308 ✗ assert(0 <= phaseIdx && phaseIdx < numPhases);
309
310 // liquid phase
311 ✗ if (phaseIdx == wPhaseIdx)
312 {
313 return 1044.0;
314 }
315 else if (phaseIdx == nPhaseIdx)// gas phase
316 {
317 return 1.679;
318 }
319 else DUNE_THROW(Dune::NotImplemented, "Wrong phase index");
320 }
321
322 using Base::molarDensity;
323 /*!
324 * \brief The molar density \f$\rho_{mol,\alpha}\f$
325 * of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
326 *
327 * This is a specific molar density for the combustion test
328 */
329 template <class FluidState>
330 ✗ static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
331 {
332 ✗ if (phaseIdx == wPhaseIdx)
333 {
334 ✗ return density(fluidState, phaseIdx)/fluidState.averageMolarMass(phaseIdx);
335 }
336 ✗ else if (phaseIdx == nPhaseIdx)
337 {
338 ✗ return density(fluidState, phaseIdx)/fluidState.averageMolarMass(phaseIdx);
339 }
340 ✗ else DUNE_THROW(Dune::NotImplemented, "Wrong phase index");
341 }
342
343 using Base::viscosity;
344 /*!
345 * \brief Calculates the dynamic viscosity of a fluid phase \f$\mathrm{[Pa*s]}\f$
346 *
347 * \param fluidState An arbitrary fluid state
348 * \param phaseIdx The index of the fluid phase to consider
349 */
350 template <class FluidState>
351 static Scalar viscosity(const FluidState &fluidState,
352 int phaseIdx)
353 {
354 ✗ assert(0 <= phaseIdx && phaseIdx < numPhases);
355
356 // liquid phase
357 ✗ if (phaseIdx == wPhaseIdx)
358 {
359 return 2.694e-7 * density(fluidState, phaseIdx);
360 }
361 else if (phaseIdx == nPhaseIdx) // gas phase
362 {
363 return 7.16e-6 * density(fluidState, phaseIdx);
364 }
365 else DUNE_THROW(Dune::NotImplemented, "Wrong phase index");
366 }
367
368 /*!
369 * \brief Calculates the temperature \f$\mathrm{[K]}\f$ of vapor at a given pressure on the vapor pressure curve.
370 *
371 * \param fluidState An arbitrary fluid state
372 * \param phaseIdx The index of the fluid phase to consider
373 */
374 template <class FluidState>
375 static Scalar vaporTemperature(const FluidState &fluidState,
376 const unsigned int phaseIdx)
377 {
378 assert(0 <= phaseIdx && phaseIdx < numPhases);
379 ✗ Scalar pressure = fluidState.pressure(nPhaseIdx);
380
381 ✗ return IAPWS::Region4<Scalar>::vaporTemperature( pressure );
382 }
383
384 using Base::fugacityCoefficient;
385 /*!
386 * \brief Calculates the fugacity coefficient \f$\mathrm{[-]}\f$ of an individual
387 * component in a fluid phase
388 *
389 * The fugacity coefficient \f$\phi^\kappa_\alpha\f$ of
390 * component \f$\kappa\f$ in phase \f$\alpha\f$ is connected to
391 * the fugacity \f$f^\kappa_\alpha\f$ and the component's mole
392 * fraction \f$x^\kappa_\alpha\f$ by means of the relation
393 *
394 * \f[
395 f^\kappa_\alpha = \phi^\kappa_\alpha\;x^\kappa_\alpha\;p_\alpha
396 \f]
397 * where \f$p_\alpha\f$ is the pressure of the fluid phase.
398 *
399 * The quantity "fugacity" itself is just an other way to express
400 * the chemical potential \f$\zeta^\kappa_\alpha\f$ of the
401 * component. It is defined via
402 *
403 * \f[
404 f^\kappa_\alpha := \exp\left\{\frac{\zeta^\kappa_\alpha}{k_B T_\alpha} \right\}
405 \f]
406 * where \f$k_B = 1.380\cdot10^{-23}\;J/K\f$ is the Boltzmann constant.
407 *
408 * \param fluidState An arbitrary fluid state
409 * \param phaseIdx The index of the fluid phase to consider
410 * \param compIdx The index of the component to consider
411 */
412 template <class FluidState>
413 ✗ static Scalar fugacityCoefficient(const FluidState &fluidState,
414 int phaseIdx,
415 int compIdx)
416 {
417 ✗ assert(0 <= phaseIdx && phaseIdx < numPhases);
418 ✗ assert(0 <= compIdx && compIdx < numComponents);
419
420 ✗ Scalar T = fluidState.temperature(phaseIdx);
421 ✗ Scalar p = fluidState.pressure(phaseIdx);
422
423 // liquid phase
424 ✗ if (phaseIdx == wPhaseIdx)
425 {
426 ✗ if (compIdx == H2OIdx)
427 ✗ return H2O::vaporPressure(T)/p;
428 ✗ return BinaryCoeff::H2O_N2::henry(T)/p;
429 }
430
431 // for the gas phase, assume an ideal gas when it comes to
432 // fugacity (-> fugacity == partial pressure)
433 return 1.0;
434 }
435
436 using Base::diffusionCoefficient;
437 /*!
438 * \brief Calculates the molecular diffusion coefficient for a
439 * component in a fluid phase \f$\mathrm{[mol^2 * s / (kg*m^3)]}\f$
440 *
441 * \param fluidState An arbitrary fluid state
442 * \param phaseIdx The index of the fluid phase to consider
443 * \param compIdx The index of the component to consider
444 */
445 template <class FluidState>
446 static Scalar diffusionCoefficient(const FluidState &fluidState,
447 int phaseIdx,
448 int compIdx)
449 {
450 DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients");
451 }
452
453 using Base::binaryDiffusionCoefficient;
454 /*!
455 * \brief Given a phase's composition, temperature and pressure,
456 * returns the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ for components
457 * \f$i\f$ and \f$j\f$ in this phase.
458 *
459 * \param fluidState An arbitrary fluid state
460 * \param phaseIdx The index of the fluid phase to consider
461 * \param compIIdx The index of the first component to consider
462 * \param compJIdx The index of the second component to consider
463 */
464 template <class FluidState>
465 static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
466 int phaseIdx,
467 int compIIdx,
468 int compJIdx)
469
470 {
471 DUNE_THROW(Dune::NotImplemented, "Binary Diffusion coefficients");
472 }
473
474 using Base::enthalpy;
475 /*!
476 * \brief Calculates specific enthalpy \f$\mathrm{[J/kg]}\f$.
477 *
478 * \param fluidState An arbitrary fluid state
479 * \param phaseIdx The index of the fluid phase to consider
480 */
481 template <class FluidState>
482 ✗ static Scalar enthalpy(const FluidState &fluidState,
483 int phaseIdx)
484 {
485 ✗ assert(0 <= phaseIdx && phaseIdx < numPhases);
486 ✗ Scalar temperature = fluidState.temperature(phaseIdx);
487
488 ✗ const Scalar cp = heatCapacity(fluidState, phaseIdx);
489
490 // liquid phase
491 ✗ if (phaseIdx == wPhaseIdx)
492 {
493 ✗ return cp * (temperature - 373.15);
494 }
495 else if (phaseIdx == nPhaseIdx) // gas phase
496 {
497 ✗ return cp * (temperature - 373.15) + 2.257e6;
498 }
499 else DUNE_THROW(Dune::NotImplemented, "Wrong phase index");
500 }
501
502 using Base::thermalConductivity;
503 /*!
504 * \brief Thermal conductivity of a fluid phase \f$\mathrm{[W/(m K)]}\f$.
505 *
506 * Use the conductivity of vapor and liquid water at 100°C
507 *
508 * \param fluidState An arbitrary fluid state
509 * \param phaseIdx The index of the fluid phase to consider
510 */
511 template <class FluidState>
512 static Scalar thermalConductivity(const FluidState &fluidState,
513 const int phaseIdx)
514 {
515 ✗ assert(0 <= phaseIdx && phaseIdx < numPhases);
516 // liquid phase
517 ✗ if (phaseIdx == wPhaseIdx)
518 {
519 return 0.68;
520 }
521 else if (phaseIdx == nPhaseIdx) // gas phase
522 {
523 return 0.0248;
524 }
525 else DUNE_THROW(Dune::NotImplemented, "Wrong phase index");
526 }
527
528 using Base::heatCapacity;
529 /*!
530 * \brief Specific isobaric heat capacity of a fluid phase.
531 * \f$\mathrm{[J/kg / K]}\f$.
532 *
533 * \param fluidState An arbitrary fluid state
534 * \param phaseIdx The index of the fluid phase to consider
535 */
536 template <class FluidState>
537 static Scalar heatCapacity(const FluidState &fluidState,
538 int phaseIdx)
539 {
540 ✗ assert(0 <= phaseIdx && phaseIdx < numPhases);
541 // liquid phase
542 ✗ if (phaseIdx == wPhaseIdx)
543 {
544 return 4.217e3;
545 }
546 else if (phaseIdx == nPhaseIdx) // gas phase
547 {
548 return 2.029e3;
549 }
550 else DUNE_THROW(Dune::NotImplemented, "Wrong phase index");
551 }
552 };
553
554 } // end namespace FluidSystems
555 } // end namespace Dumux
556
557 #endif
558