GCC Code Coverage Report

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