GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/material/fluidsystems/brineair.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 126 134 94.0%
Functions: 36 36 100.0%
Branches: 143 340 42.1%

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::BrineAir
11 */
12 #ifndef DUMUX_BRINE_AIR_FLUID_SYSTEM_HH
13 #define DUMUX_BRINE_AIR_FLUID_SYSTEM_HH
14
15 #include <array>
16 #include <cassert>
17 #include <iomanip>
18
19 #include <dumux/material/idealgas.hh>
20 #include <dumux/material/fluidstates/adapter.hh>
21 #include <dumux/material/fluidsystems/base.hh>
22 #include <dumux/material/fluidsystems/brine.hh>
23 #include <dumux/material/components/air.hh>
24 #include <dumux/material/components/h2o.hh>
25 #include <dumux/material/components/nacl.hh>
26 #include <dumux/material/binarycoefficients/h2o_air.hh>
27 #include <dumux/material/components/tabulatedcomponent.hh>
28
29 #include <dumux/common/exceptions.hh>
30
31 #include <dumux/io/name.hh>
32
33 #include "brine.hh"
34
35 namespace Dumux::FluidSystems {
36
37 /*!
38 * \ingroup FluidSystems
39 * \brief Policy for the brine-air fluid system
40 */
41 template<bool fastButSimplifiedRelations = false>
42 struct BrineAirDefaultPolicy
43 {
44 static constexpr bool useBrineDensityAsLiquidMixtureDensity() { return fastButSimplifiedRelations;}
45 static constexpr bool useIdealGasDensity() { return fastButSimplifiedRelations; }
46 static constexpr bool useKelvinVaporPressure() { return false; }
47 };
48
49 /*!
50 * \ingroup FluidSystems
51 * \brief A compositional two-phase fluid system with a liquid and a gaseous phase
52 * and \f$H_2O\f$, \f$Air\f$ and \f$S\f$ (dissolved minerals) as components.
53 *
54 * \note This fluidsystem is applied by default with the tabulated version of
55 * water of the IAPWS-formulation.
56 */
57 template <class Scalar,
58 class H2Otype = Components::TabulatedComponent<Components::H2O<Scalar>>,
59 class Policy = BrineAirDefaultPolicy<>>
60 class BrineAir
61 : public Base<Scalar, BrineAir<Scalar, H2Otype, Policy>>
62 {
63 using ThisType = BrineAir<Scalar, H2Otype, Policy>;
64 using IdealGas = Dumux::IdealGas<Scalar>;
65
66 public:
67 //! export the involved components
68 using H2O = H2Otype;
69 using Air = Components::Air<Scalar>;
70 using NaCl = Components::NaCl<Scalar>;
71
72 //! export the underlying brine fluid system for the liquid phase
73 using Brine = Dumux::FluidSystems::Brine<Scalar, H2Otype>;
74
75 //! export the binary coefficients between air and water
76 using H2O_Air = BinaryCoeff::H2O_Air;
77
78 //! the type of parameter cache objects
79 using ParameterCache = NullParameterCache;
80
81 /****************************************
82 * Fluid phase related static parameters
83 ****************************************/
84 static constexpr int numPhases = 2; // one liquid and one gas phase
85 static constexpr int numComponents = 3; // H2O, Air, NaCl
86
87 static constexpr int liquidPhaseIdx = 0; // index of the liquid phase
88 static constexpr int gasPhaseIdx = 1; // index of the gas phase
89
90 static constexpr int phase0Idx = liquidPhaseIdx; // index of the first phase
91 static constexpr int phase1Idx = gasPhaseIdx; // index of the second phase
92
93 // export component indices to indicate the main component
94 // of the corresponding phase at atmospheric pressure 1 bar
95 // and room temperature 20°C:
96 static constexpr int H2OIdx = 0;
97 static constexpr int AirIdx = 1;
98 static constexpr int NaClIdx = 2;
99 static constexpr int comp0Idx = H2OIdx;
100 static constexpr int comp1Idx = AirIdx;
101 static constexpr int comp2Idx = NaClIdx;
102
103 private:
104 struct BrineAdapterPolicy
105 {
106 using FluidSystem = Brine;
107
108 static constexpr int phaseIdx(int brinePhaseIdx) { return liquidPhaseIdx; }
109 33258547 static constexpr int compIdx(int brineCompIdx)
110 {
111
17/20
✓ Branch 0 taken 5022751 times.
✓ Branch 1 taken 5022751 times.
✓ Branch 2 taken 7364402 times.
✓ Branch 3 taken 7364401 times.
✓ Branch 4 taken 7364403 times.
✓ Branch 5 taken 88182 times.
✓ Branch 6 taken 343878 times.
✓ Branch 7 taken 343878 times.
✓ Branch 8 taken 343878 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 1 times.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✓ Branch 15 taken 1 times.
✓ Branch 16 taken 6 times.
✓ Branch 17 taken 6 times.
✓ Branch 18 taken 6 times.
✗ Branch 19 not taken.
33258547 switch (brineCompIdx)
112 {
113 case Brine::H2OIdx: return H2OIdx;
114 20439327 case Brine::NaClIdx: return NaClIdx;
115 default: return 0; // this will never be reached, only needed to suppress compiler warning
116 }
117 }
118 };
119
120 template<class FluidState>
121 using BrineAdapter = FluidStateAdapter<FluidState, BrineAdapterPolicy>;
122
123 public:
124
125 /****************************************
126 * phase related static parameters
127 ****************************************/
128
129 /*!
130 * \brief Return the human readable name of a fluid phase
131 * \param phaseIdx index of the phase
132 */
133 8 static std::string phaseName(int phaseIdx)
134 {
135
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 460 times.
464 assert(0 <= phaseIdx && phaseIdx < numPhases);
136
2/2
✓ Branch 0 taken 230 times.
✓ Branch 1 taken 230 times.
464 switch (phaseIdx)
137 {
138 232 case liquidPhaseIdx: return IOName::liquidPhase();
139 232 case gasPhaseIdx: return IOName::gaseousPhase();
140 }
141 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
142 }
143
144 /*!
145 * \brief Returns whether the fluids are miscible
146 */
147 static constexpr bool isMiscible()
148 { return true; }
149
150 /*!
151 * \brief Return whether a phase is gaseous
152 * \param phaseIdx The index of the fluid phase to consider
153 */
154 4 static constexpr bool isGas(int phaseIdx)
155 {
156
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
4 assert(0 <= phaseIdx && phaseIdx < numPhases);
157
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
4 return phaseIdx == gasPhaseIdx;
158 }
159
160 /*!
161 * \brief Returns true if and only if a fluid phase is assumed to
162 * be an ideal mixture.
163 *
164 * We define an ideal mixture as a fluid phase where the fugacity
165 * coefficients of all components times the pressure of the phase
166 * are independent on the fluid composition. This assumption is true
167 * if Henry's law and Raoult's law apply. If you are unsure what
168 * this function should return, it is safe to return false. The
169 * only damage done will be (slightly) increased computation times
170 * in some cases.
171 *
172 * \param phaseIdx The index of the fluid phase to consider
173 */
174 561825 static bool isIdealMixture(int phaseIdx)
175 {
176
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 561825 times.
561825 assert(0 <= phaseIdx && phaseIdx < numPhases);
177 // we assume Henry's and Raoult's laws for the water phase and
178 // and no interaction between gas molecules of different
179 // components, so all phases are ideal mixtures!
180 return true;
181 }
182
183 /*!
184 * \brief Returns true if and only if a fluid phase is assumed to
185 * be compressible.
186 *
187 * Compressible means that the partial derivative of the density
188 * to the fluid pressure is always larger than zero.
189 *
190 * \param phaseIdx The index of the fluid phase to consider
191 */
192 4 static constexpr bool isCompressible(int phaseIdx)
193 {
194 assert(0 <= phaseIdx && phaseIdx < numPhases);
195 // ideal gases are always compressible
196
4/4
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 1 times.
4 if (phaseIdx == gasPhaseIdx)
197 2 return true;
198 // let brine decide for the liquid phase...
199 return Brine::isCompressible(Brine::liquidPhaseIdx);
200 }
201
202 /*!
203 * \brief Returns true if and only if a fluid phase is assumed to
204 * be an ideal gas.
205 * \param phaseIdx The index of the fluid phase to consider
206 */
207 static bool isIdealGas(int phaseIdx)
208 {
209 assert(0 <= phaseIdx && phaseIdx < numPhases);
210 // let the fluids decide
211 if (phaseIdx == gasPhaseIdx)
212 return H2O::gasIsIdeal() && Air::gasIsIdeal();
213 return false; // not a gas
214 }
215
216 /*!
217 * \brief Get the main component of a given phase if possible
218 * \param phaseIdx The index of the fluid phase to consider
219 */
220 80434884 static constexpr int getMainComponent(int phaseIdx)
221 {
222
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 32273976 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 36162240 times.
68436216 assert(0 <= phaseIdx && phaseIdx < numPhases);
223
6/6
✓ Branch 0 taken 16136988 times.
✓ Branch 1 taken 16136988 times.
✓ Branch 2 taken 18081120 times.
✓ Branch 3 taken 18081120 times.
✓ Branch 4 taken 5999334 times.
✓ Branch 5 taken 5999334 times.
80434884 if (phaseIdx == liquidPhaseIdx)
224 return H2OIdx;
225 else
226 40217442 return AirIdx;
227 }
228
229 /****************************************
230 * Component related static parameters
231 ****************************************/
232 /*!
233 * \brief Return the human readable name of a component
234 * \param compIdx The index of the component to consider
235 */
236 52 static std::string componentName(int compIdx)
237 {
238
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 46 times.
52 assert(0 <= compIdx && compIdx < numComponents);
239
3/3
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 18 times.
52 switch (compIdx)
240 {
241 16 case H2OIdx: return H2O::name();
242 16 case AirIdx: return Air::name();
243 20 case NaClIdx: return NaCl::name();
244 }
245 DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
246 }
247
248 /*!
249 * \brief Return the molar mass of a component in \f$\mathrm{[kg/mol]}\f$.
250 * \param compIdx The index of the component to consider
251 */
252
2/2
✓ Branch 0 taken 378 times.
✓ Branch 1 taken 53 times.
4448147 static Scalar molarMass(int compIdx)
253 {
254
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 184491201 times.
184491291 assert(0 <= compIdx && compIdx < numComponents);
255
3/3
✓ Branch 0 taken 47767000 times.
✓ Branch 1 taken 81320184 times.
✓ Branch 2 taken 55404017 times.
184491291 switch (compIdx)
256 {
257 case H2OIdx: return H2O::molarMass();
258 47767020 case AirIdx: return Air::molarMass();
259 81320222 case NaClIdx: return NaCl::molarMass();
260 }
261 DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
262 }
263
264 /*!
265 * \brief Vapor pressure of a component \f$\mathrm{[Pa]}\f$.
266 *
267 * \param fluidState The fluid state
268 * \param compIdx The index of the component to consider
269 */
270 template <class FluidState>
271 static Scalar vaporPressure(const FluidState& fluidState, int compIdx)
272 {
273 // The vapor pressure of the water is affected by the
274 // salinity, thus, we forward to the interface of Brine here
275 if (compIdx == H2OIdx)
276 {
277 if (!Policy::useKelvinVaporPressure())
278 return Brine::vaporPressure(BrineAdapter<FluidState>(fluidState), Brine::H2OIdx);
279 else
280 {
281 //additionally taking the influence of the capillary pressure on the vapor pressure, described by Kelvin's equation, into account.
282 const auto t = fluidState.temperature(liquidPhaseIdx);
283 const auto pc = (fluidState.wettingPhase() == (int) liquidPhaseIdx)
284 ? fluidState.pressure(gasPhaseIdx)-fluidState.pressure(liquidPhaseIdx)
285 : fluidState.pressure(liquidPhaseIdx)-fluidState.pressure(gasPhaseIdx);
286 //influence of the capillary pressure on the vapor pressure, the influence of the salt concentration is already taken into account in the brine vapour pressure
287 return Brine::vaporPressure(BrineAdapter<FluidState>(fluidState), Brine::H2OIdx)
288 *exp( -pc / (Brine::molarDensity(fluidState, liquidPhaseIdx) * (Dumux::Constants<Scalar>::R*t)));
289 }
290 }
291 else if (compIdx == NaClIdx)
292 DUNE_THROW(Dune::NotImplemented, "NaCl::vaporPressure(t)");
293 else
294 DUNE_THROW(Dune::NotImplemented, "Invalid component index " << compIdx);
295 }
296
297 /****************************************
298 * thermodynamic relations
299 ****************************************/
300 /*!
301 * \brief Initialize the fluid system's static parameters generically
302 *
303 * If a tabulated H2O component is used, we do our best to create
304 * tables that always work.
305 */
306 2 static void init()
307 {
308 2 init(/*tempMin=*/273.15,
309 /*tempMax=*/800.0,
310 /*numTemptempSteps=*/200,
311 /*startPressure=*/-10,
312 /*endPressure=*/20e6,
313 /*pressureSteps=*/200);
314 }
315
316 /*!
317 * \brief Initialize the fluid system's static parameters using
318 * problem specific temperature and pressure ranges
319 *
320 * \param tempMin The minimum temperature used for tabulation of water \f$\mathrm{[K]}\f$
321 * \param tempMax The maximum temperature used for tabulation of water \f$\mathrm{[K]}\f$
322 * \param nTemp The number of ticks on the temperature axis of the table of water
323 * \param pressMin The minimum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$
324 * \param pressMax The maximum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$
325 * \param nPress The number of ticks on the pressure axis of the table of water
326 */
327 9 static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
328 Scalar pressMin, Scalar pressMax, unsigned nPress)
329 {
330 9 std::cout << "The brine-air fluid system was configured with the following policy:\n";
331 9 std::cout << " - use brine density as liquid mixture density: " << std::boolalpha << Policy::useBrineDensityAsLiquidMixtureDensity() << "\n";
332 9 std::cout << " - use ideal gas density: " << std::boolalpha << Policy::useIdealGasDensity() << std::endl;
333
334 if (H2O::isTabulated)
335 H2O::init(tempMin, tempMax, nTemp, pressMin, pressMax, nPress);
336 9 }
337
338 using Base<Scalar, ThisType>::density;
339 /*!
340 * \brief Given a phase's composition, temperature, pressure, and
341 * the partial pressures of all components, return its
342 * density \f$\mathrm{[kg/m^3]}\f$.
343 *
344 * \param fluidState the fluid state
345 * \param phaseIdx index of the phase
346 *
347 * Equation given in:
348 * - Batzle & Wang (1992) \cite batzle1992
349 * - cited by: Bachu & Adams (2002)
350 * "Equations of State for basin geofluids" \cite adams2002
351 */
352 template <class FluidState>
353 5022758 static Scalar density(const FluidState &fluidState, int phaseIdx)
354 {
355
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5022754 times.
5022758 assert(0 <= phaseIdx && phaseIdx < numPhases);
356
357
2/2
✓ Branch 0 taken 2511375 times.
✓ Branch 1 taken 2511375 times.
5022758 const auto T = fluidState.temperature(phaseIdx);
358 5022758 const auto p = fluidState.pressure(phaseIdx);
359
360
2/2
✓ Branch 0 taken 2511377 times.
✓ Branch 1 taken 2511377 times.
5022758 if (phaseIdx == liquidPhaseIdx)
361 {
362 // assume pure brine
363 if (Policy::useBrineDensityAsLiquidMixtureDensity())
364 2 return Brine::density(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx);
365
366 // assume one molecule of gas replaces one "brine" molecule
367 else
368 2511377 return Brine::molarDensity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx)
369 2511377 *(H2O::molarMass()*fluidState.moleFraction(liquidPhaseIdx, H2OIdx)
370 2511377 + NaCl::molarMass()*fluidState.moleFraction(liquidPhaseIdx, NaClIdx)
371 2511377 + Air::molarMass()*fluidState.moleFraction(liquidPhaseIdx, AirIdx));
372 }
373 else if (phaseIdx == phase1Idx)
374 {
375 // for the gas phase assume an ideal gas
376 if (Policy::useIdealGasDensity())
377 2 return IdealGas::density(fluidState.averageMolarMass(phase1Idx), T, p);
378
379 // if useComplexRelations = true, compute density. NaCl is assumed
380 // not to be present in gas phase, NaCl has only solid interfaces implemented
381 2511377 return H2O::gasDensity(T, fluidState.partialPressure(phase1Idx, H2OIdx))
382 2511377 + Air::gasDensity(T, fluidState.partialPressure(phase1Idx, AirIdx));
383 }
384 else
385 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
386 }
387
388 using Base<Scalar, ThisType>::molarDensity;
389 //! \copydoc Base<Scalar,ThisType>::molarDensity(const FluidState&,int)
390 template <class FluidState>
391 5022758 static Scalar molarDensity(const FluidState& fluidState, int phaseIdx)
392 {
393
2/2
✓ Branch 0 taken 2511377 times.
✓ Branch 1 taken 2511377 times.
5022758 if (phaseIdx == liquidPhaseIdx)
394 2511379 return Brine::molarDensity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx);
395
1/2
✓ Branch 0 taken 2511377 times.
✗ Branch 1 not taken.
2511379 else if (phaseIdx == phase1Idx)
396 {
397 2511379 const Scalar T = fluidState.temperature(phaseIdx);
398
399 // for the gas phase assume an ideal gas
400 if (Policy::useIdealGasDensity())
401 2 return IdealGas::molarDensity(T, fluidState.pressure(phaseIdx));
402
403 // if useComplexRelations = true, compute density. NaCl is assumed
404 // not to be present in gas phase, NaCl has only solid interfaces implemented
405 2511377 return H2O::gasMolarDensity(T, fluidState.partialPressure(phase1Idx, H2OIdx))
406 2511377 + Air::gasMolarDensity(T, fluidState.partialPressure(phase1Idx, AirIdx));
407 }
408 else
409 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
410 }
411
412 using Base<Scalar, ThisType>::viscosity;
413 /*!
414 * \brief Calculate the dynamic viscosity of a fluid phase \f$\mathrm{[Pa*s]}\f$
415 *
416 * \param fluidState An arbitrary fluid state
417 * \param phaseIdx The index of the fluid phase to consider
418 *
419 * \note For the viscosity of the phases the contribution of the minor
420 * component is neglected. This contribution is probably not big, but somebody
421 * would have to find out its influence.
422 */
423 template <class FluidState>
424 5022758 static Scalar viscosity(const FluidState& fluidState, int phaseIdx)
425 {
426
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5022754 times.
5022758 assert(0 <= phaseIdx && phaseIdx < numPhases);
427
428
2/2
✓ Branch 0 taken 2511377 times.
✓ Branch 1 taken 2511377 times.
5022758 if (phaseIdx == liquidPhaseIdx)
429 2511379 return Brine::viscosity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx);
430 else
431 2511379 return Air::gasViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
432 }
433
434 using Base<Scalar, ThisType>::fugacityCoefficient;
435 /*!
436 * \brief Returns the fugacity coefficient \f$\mathrm{[-]}\f$ of a component in a
437 * phase.
438 * \param fluidState The fluid state
439 * \param phaseIdx Index of the phase
440 * \param compIdx Index of the component
441 *
442 * The fugacity coefficient \f$\mathrm{\phi^\kappa_\alpha}\f$ of
443 * component \f$\mathrm{\kappa}\f$ in phase \f$\mathrm{\alpha}\f$ is connected to
444 * the fugacity \f$\mathrm{f^\kappa_\alpha}\f$ and the component's mole
445 * fraction \f$\mathrm{x^\kappa_\alpha}\f$ by means of the relation
446 *
447 * \f[
448 f^\kappa_\alpha = \phi^\kappa_\alpha\;x^\kappa_\alpha\;p_\alpha
449 \f]
450 * where \f$\mathrm{p_\alpha}\f$ is the pressure of the fluid phase.
451 *
452 * For liquids with very low miscibility this boils down to the
453 * Henry constant for the solutes and the saturated vapor pressure
454 * both divided by phase pressure.
455 */
456 template <class FluidState>
457 15068274 static Scalar fugacityCoefficient(const FluidState& fluidState, int phaseIdx, int compIdx)
458 {
459
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15068262 times.
15068274 assert(0 <= phaseIdx && phaseIdx < numPhases);
460
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15068262 times.
15068274 assert(0 <= compIdx && compIdx < numComponents);
461
462
2/2
✓ Branch 0 taken 7534125 times.
✓ Branch 1 taken 7534125 times.
15068274 Scalar T = fluidState.temperature(phaseIdx);
463 15068274 Scalar p = fluidState.pressure(phaseIdx);
464
465
2/2
✓ Branch 0 taken 7534131 times.
✓ Branch 1 taken 7534131 times.
15068274 if (phaseIdx == gasPhaseIdx)
466 return 1.0;
467
468 else if (phaseIdx == liquidPhaseIdx)
469 {
470 // TODO: Should we use the vapor pressure of the mixture (brine) here?
471 // The presence of NaCl lowers the vapor pressure, thus, we would
472 // expect the fugacity coefficient to be lower as well. However,
473 // with the fugacity coefficient being dependent on the salinity,
474 // the equation system for the phase equilibria becomes non-linear
475 // and our constraint solvers assume linear system of equations.
476
2/2
✓ Branch 0 taken 2511377 times.
✓ Branch 1 taken 5022754 times.
7534137 if (compIdx == H2OIdx)
477 2511379 return H2O::vaporPressure(T)/p;
478
479
2/2
✓ Branch 0 taken 2511377 times.
✓ Branch 1 taken 2511377 times.
5022758 else if (compIdx == AirIdx)
480 2511379 return BinaryCoeff::H2O_Air::henry(T)/p;
481
482 // we assume nacl always stays in the liquid phase
483 else
484 return 0.0;
485 }
486
487 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
488 }
489
490 using Base<Scalar, ThisType>::diffusionCoefficient;
491 //! \copydoc Base<Scalar,ThisType>::diffusionCoefficient(const FluidState&,int,int)
492 template <class FluidState>
493 24 static Scalar diffusionCoefficient(const FluidState &fluidState,
494 int phaseIdx,
495 int compIdx)
496 {
497
10/20
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 12 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 12 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 12 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 12 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 12 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 12 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 12 times.
✗ Branch 29 not taken.
96 DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients");
498 }
499
500 using Base<Scalar, ThisType>::binaryDiffusionCoefficient;
501 //! \copydoc Base<Scalar,ThisType>::binaryDiffusionCoefficient(const FluidState&,int,int,int)
502 template <class FluidState>
503 10064696 static Scalar binaryDiffusionCoefficient(const FluidState& fluidState,
504 int phaseIdx,
505 int compIIdx,
506 int compJIdx)
507 {
508
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10064660 times.
10064696 assert(0 <= phaseIdx && phaseIdx < numPhases);
509
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10064660 times.
10064696 assert(0 <= compIIdx && compIIdx < numComponents);
510
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10064660 times.
10064696 assert(0 <= compJIdx && compJIdx < numComponents);
511
512
2/2
✓ Branch 0 taken 2530499 times.
✓ Branch 1 taken 7534125 times.
10064696 const auto T = fluidState.temperature(phaseIdx);
513 10064696 const auto p = fluidState.pressure(phaseIdx);
514
515
2/2
✓ Branch 0 taken 2530511 times.
✓ Branch 1 taken 7534149 times.
10064696 if (compIIdx > compJIdx)
516 2530523 std::swap(compIIdx, compJIdx);
517
518
2/2
✓ Branch 0 taken 5022768 times.
✓ Branch 1 taken 5041892 times.
10064696 if (phaseIdx == liquidPhaseIdx)
519 {
520
4/4
✓ Branch 0 taken 5022760 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2511379 times.
✓ Branch 3 taken 2511381 times.
5022786 if(compIIdx == H2OIdx && compJIdx == AirIdx)
521 2511383 return H2O_Air::liquidDiffCoeff(T, p);
522
2/2
✓ Branch 0 taken 2511379 times.
✓ Branch 1 taken 2 times.
2511387 else if (compIIdx == H2OIdx && compJIdx == NaClIdx)
523 2511383 return Brine::binaryDiffusionCoefficient(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx, Brine::H2OIdx, Brine::NaClIdx);
524 else
525
15/30
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 10 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 10 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 10 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 10 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 10 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 10 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 10 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 10 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 10 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 10 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 10 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 10 times.
✗ Branch 44 not taken.
100 DUNE_THROW(Dune::NotImplemented, "Binary diffusion coefficient of components "
526 << compIIdx << " and " << compJIdx
527 << " in phase " << phaseIdx);
528 }
529 else if (phaseIdx == gasPhaseIdx)
530 {
531
4/4
✓ Branch 0 taken 2530509 times.
✓ Branch 1 taken 2511383 times.
✓ Branch 2 taken 2530503 times.
✓ Branch 3 taken 6 times.
5041910 if (compIIdx == H2OIdx && compJIdx == AirIdx)
532 2530507 return H2O_Air::gasDiffCoeff(T, p);
533
534 // NaCl is expected to never be present in the gas phase. we need to
535 // return a diffusion coefficient that does not case numerical problems.
536 // We choose a very small value here.
537
4/4
✓ Branch 0 taken 2511381 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2511379 times.
2511403 else if (compIIdx == AirIdx && compJIdx == NaClIdx)
538 return 1e-12;
539
540 else
541
15/30
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 10 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 10 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 10 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 10 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 10 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 10 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 10 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 10 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 10 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 10 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 10 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 10 times.
✗ Branch 44 not taken.
80 DUNE_THROW(Dune::NotImplemented, "Binary diffusion coefficient of components "
542 << compIIdx << " and " << compJIdx
543 << " in phase " << phaseIdx);
544 }
545
546 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
547 }
548
549 using Base<Scalar, ThisType>::enthalpy;
550 /*!
551 * \brief Given a phase's composition, temperature and pressure,
552 * return its specific enthalpy \f$\mathrm{[J/kg]}\f$.
553 * \param fluidState The fluid state
554 * \param phaseIdx The index of the phase
555 *
556 * See:
557 * Class 2000
558 * Theorie und numerische Modellierung nichtisothermer Mehrphasenprozesse in NAPL-kontaminierten porösen Medien
559 * Chapter 2.1.13 Innere Energie, Wäremekapazität, Enthalpie \cite A3:class:2001
560 *
561 * Formula (2.42):
562 * the specific enthalpy of a gas phase result from the sum of (enthalpies*mass fraction) of the components
563 * For the calculation of enthalpy of brine we refer to (Michaelides 1981)
564 *
565 * \note For the phase enthalpy the contribution of gas-molecules in the liquid phase
566 * is neglected. This contribution is probably not big. Somebody would have to
567 * find out the enthalpy of solution for this system. ...
568 */
569 template <class FluidState>
570 171944 static Scalar enthalpy(const FluidState& fluidState, int phaseIdx)
571 {
572
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 171940 times.
171944 assert(0 <= phaseIdx && phaseIdx < numPhases);
573
574
2/2
✓ Branch 0 taken 85968 times.
✓ Branch 1 taken 85968 times.
171944 const Scalar T = fluidState.temperature(phaseIdx);
575 171944 const Scalar p = fluidState.pressure(phaseIdx);
576
577
2/2
✓ Branch 0 taken 85970 times.
✓ Branch 1 taken 85970 times.
171944 if (phaseIdx == liquidPhaseIdx)
578 85972 return Brine::enthalpy(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx);
579 else if (phaseIdx == gasPhaseIdx)
580 {
581 // This assumes NaCl not to be present in the gas phase
582 85972 const Scalar XAir = fluidState.massFraction(gasPhaseIdx, AirIdx);
583 85972 const Scalar XH2O = fluidState.massFraction(gasPhaseIdx, H2OIdx);
584 85972 return XH2O*H2O::gasEnthalpy(T, p) + XAir*Air::gasEnthalpy(T, p);
585 }
586
587 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
588 }
589
590 /*!
591 * \brief Returns the specific enthalpy \f$\mathrm{[J/kg]}\f$ of a component in a specific phase
592 * \param fluidState The fluid state
593 * \param phaseIdx The index of the phase
594 * \param componentIdx The index of the component
595 */
596 template <class FluidState>
597
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 16912 times.
33824 static Scalar componentEnthalpy(const FluidState& fluidState, int phaseIdx, int componentIdx)
598 {
599
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 16912 times.
33824 const Scalar T = fluidState.temperature(gasPhaseIdx);
600 33824 const Scalar p = fluidState.pressure(gasPhaseIdx);
601
602
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 16912 times.
16912 if (phaseIdx == liquidPhaseIdx)
603 DUNE_THROW(Dune::NotImplemented, "The component enthalpies in the liquid phase are not implemented.");
604
605
1/2
✓ Branch 0 taken 16912 times.
✗ Branch 1 not taken.
16912 else if (phaseIdx == gasPhaseIdx)
606 {
607
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 16912 times.
16912 if (componentIdx == H2OIdx)
608 16912 return H2O::gasEnthalpy(T, p);
609
1/2
✓ Branch 0 taken 16912 times.
✗ Branch 1 not taken.
16912 else if (componentIdx == AirIdx)
610 16912 return Air::gasEnthalpy(T, p);
611 else if (componentIdx == NaClIdx)
612 DUNE_THROW(Dune::InvalidStateException, "Implementation assumes NaCl not to be present in gas phase");
613 DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << componentIdx);
614 }
615 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
616 }
617
618 using Base<Scalar, ThisType>::thermalConductivity;
619 /*!
620 * \brief Thermal conductivity of a fluid phase \f$\mathrm{[W/(m K)]}\f$.
621 * \param fluidState An arbitrary fluid state
622 * \param phaseIdx The index of the fluid phase to consider
623 *
624 * \note We assume an ideal mixture for the gas-phase thermal conductivity, as a first approximation.
625 * In porous media, gas-phase thermal conductivities are negligible, as they are
626 * orders of magnitude lower than the thermal conductivity in solid and liquids (gas: 0.0x compared to solid: x.).
627 * However, moisture can have a significant influce on the thermal conductivity of moist air, see e.g. \cite Beirao2012,
628 * which could be relevant for free-flow systems.
629 */
630 template <class FluidState>
631 193280 static Scalar thermalConductivity(const FluidState& fluidState, int phaseIdx)
632 {
633
2/2
✓ Branch 0 taken 88182 times.
✓ Branch 1 taken 105094 times.
193280 if (phaseIdx == liquidPhaseIdx)
634 88184 return Brine::thermalConductivity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx);
635 // This assumes NaCl not to be present in the gas phase
636
1/2
✓ Branch 0 taken 105094 times.
✗ Branch 1 not taken.
105096 else if (phaseIdx == gasPhaseIdx)
637 105096 return Air::gasThermalConductivity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)) * fluidState.massFraction(gasPhaseIdx, AirIdx)
638 105096 + H2O::gasThermalConductivity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)) * fluidState.massFraction(gasPhaseIdx, H2OIdx);
639
640 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
641 }
642
643 using Base<Scalar, ThisType>::heatCapacity;
644 /*!
645 * \brief Specific isobaric heat capacity of a fluid phase.
646 * \f$\mathrm{[J/(kg*K)}\f$.
647 * \param fluidState An arbitrary fluid state
648 * \param phaseIdx The index of the fluid phase to consider
649 *
650 * \note We assume an ideal mixture for the gas-phase specific isobaric heat capacity, as a first approximation.
651 * In porous media, gas-phase heat capacities are, due to the usually low densities of gases, negligible anyway.
652 * A better implementation might be relevant for free-flow systems.
653 */
654 template <class FluidState>
655 4 static Scalar heatCapacity(const FluidState &fluidState, int phaseIdx)
656 {
657 4 const Scalar T = fluidState.temperature(phaseIdx);
658 4 const Scalar p = fluidState.pressure(phaseIdx);
659
660 4 if (phaseIdx == liquidPhaseIdx)
661 2 return Brine::heatCapacity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx);
662
663 // We assume NaCl not to be present in the gas phase here
664 2 else if (phaseIdx == gasPhaseIdx)
665 2 return Air::gasHeatCapacity(T, p)*fluidState.massFraction(gasPhaseIdx, AirIdx)
666 2 + H2O::gasHeatCapacity(T, p)*fluidState.massFraction(gasPhaseIdx, H2OIdx);
667
668 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
669 }
670 };
671
672 } // end namespace Dumux::FluidSystems
673
674 #endif
675