GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/porousmediumflow/nonequilibrium/volumevariables.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 116 137 84.7%
Functions: 13 36 36.1%
Branches: 17 45 37.8%

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 NonEquilibriumModel
10 * \brief This class contains the volume variables required for the
11 * modules which require the specific interfacial area between
12 * fluid phases.
13 *
14 * This files contains all specializations which use 'real'
15 * interfacial areas.
16 */
17
18 #ifndef DUMUX_NONEQUILIBRIUM_VOLUME_VARIABLES_HH
19 #define DUMUX_NONEQUILIBRIUM_VOLUME_VARIABLES_HH
20
21 #include <cassert>
22 #include <array>
23
24 #include <dumux/common/dimensionlessnumbers.hh>
25 #include <dumux/common/parameters.hh>
26
27 namespace Dumux {
28
29 /*!
30 * \ingroup NonEquilibriumModel
31 * \brief This class contains the volume variables required for the
32 * modules which require the specific interfacial area between
33 * fluid phases.
34 */
35 template<class Traits, class EquilibriumVolumeVariables, bool enableChemicalNonEquilibrium,
36 bool enableThermalNonEquilibrium, int numEnergyEqFluid>
37 class NonEquilibriumVolumeVariablesImplementation;
38
39 template<class Traits, class EquilibriumVolumeVariables>
40 using NonEquilibriumVolumeVariables =
41 NonEquilibriumVolumeVariablesImplementation<Traits,
42 EquilibriumVolumeVariables,
43 Traits::ModelTraits::enableChemicalNonEquilibrium(),
44 Traits::ModelTraits::enableThermalNonEquilibrium(),
45 Traits::ModelTraits::numEnergyEqFluid()>;
46
47 /////////////////////////////////////////////////////////////////////////////////////////////////////////////
48 // specialization for the case of NO kinetic mass but kinetic energy transfer of two fluids and a solid
49 /////////////////////////////////////////////////////////////////////////////////////////////////////////////
50 template<class Traits, class EquilibriumVolumeVariables>
51 class NonEquilibriumVolumeVariablesImplementation<Traits,
52 EquilibriumVolumeVariables,
53 false/*chemicalNonEquilibrium?*/,
54 true/*thermalNonEquilibrium?*/,
55 2>
56 : public EquilibriumVolumeVariables
57 {
58 using ParentType = EquilibriumVolumeVariables;
59 using ParameterCache = typename Traits::FluidSystem::ParameterCache;
60 using Scalar = typename Traits::PrimaryVariables::value_type;
61
62 using ModelTraits = typename Traits::ModelTraits;
63
64 using FS = typename Traits::FluidSystem;
65 static constexpr auto numEnergyEqFluid = ModelTraits::numEnergyEqFluid();
66 static constexpr auto numEnergyEqSolid = ModelTraits::numEnergyEqSolid();
67
68 static constexpr auto phase0Idx = FS::phase0Idx;
69 static constexpr auto phase1Idx = FS::phase1Idx;
70 static constexpr auto sPhaseIdx = FS::numPhases;
71
72 using DimLessNum = DimensionlessNumbers<Scalar>;
73
74 using NumFluidPhasesArray = std::array<Scalar, ModelTraits::numFluidPhases()>;
75 using InterfacialAreasArray = std::array<std::array<Scalar, ModelTraits::numFluidPhases()+numEnergyEqSolid>,
76 ModelTraits::numFluidPhases()+numEnergyEqSolid>;
77
78 public:
79 using FluidState = typename Traits::FluidState;
80 using Indices = typename ModelTraits::Indices;
81 /*!
82 * \brief Update all quantities for a given control volume
83 *
84 * \param elemSol A vector containing all primary variables connected to the element
85 * \param problem The object specifying the problem which ought to be simulated
86 * \param element An element which contains part of the control volume
87 * \param scv The sub-control volume
88 */
89 template<class ElemSol, class Problem, class Element, class Scv>
90 void update(const ElemSol &elemSol,
91 const Problem &problem,
92 const Element &element,
93 const Scv& scv)
94 {
95 // Update parent type (also completes the fluid state)
96 ParentType::update(elemSol, problem, element, scv);
97
98 ParameterCache paramCache;
99 paramCache.updateAll(this->fluidState());
100 updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
101 updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
102 }
103
104 /*!
105 * \brief Updates dimensionless numbers
106 *
107 * \param elemSol A vector containing all primary variables connected to the element
108 * \param fluidState Container for all the secondary variables concerning the fluids
109 * \param paramCache The parameter cache corresponding to the fluid state
110 * \param problem The problem to be solved
111 * \param element An element which contains part of the control volume
112 * \param scv The sub-control volume
113 */
114 template<class ElemSol, class Problem, class Element, class Scv>
115 void updateDimLessNumbers(const ElemSol& elemSol,
116 const FluidState& fluidState,
117 const ParameterCache& paramCache,
118 const Problem& problem,
119 const Element& element,
120 const Scv& scv)
121 {
122 const auto& spatialParams = problem.spatialParams();
123 factorEnergyTransfer_ = spatialParams.factorEnergyTransfer(element, scv, elemSol);
124 characteristicLength_ = spatialParams.characteristicLength(element, scv, elemSol);
125
126 // set the dimensionless numbers and obtain respective quantities
127 const unsigned int vIdxGlobal = scv.dofIndex();
128 for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
129 {
130 const auto darcyMagVelocity = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, vIdxGlobal);
131 const auto dynamicViscosity = fluidState.viscosity(phaseIdx);
132 const auto density = fluidState.density(phaseIdx);
133 const auto kinematicViscosity = dynamicViscosity/density;
134
135 using FluidSystem = typename Traits::FluidSystem;
136 const auto heatCapacity = FluidSystem::heatCapacity(fluidState, paramCache, phaseIdx);
137 const auto thermalConductivity = FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);
138 const auto porosity = this->porosity();
139
140 reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_,kinematicViscosity);
141 prandtlNumber_[phaseIdx] = DimLessNum::prandtlNumber(dynamicViscosity, heatCapacity, thermalConductivity);
142 nusseltNumber_[phaseIdx] = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
143 prandtlNumber_[phaseIdx],
144 porosity,
145 ModelTraits::nusseltFormulation());
146
147 }
148 }
149
150 /*!
151 * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
152 *
153 * \param elemSol A vector containing all primary variables connected to the element
154 * \param fluidState Container for all the secondary variables concerning the fluids
155 * \param paramCache The parameter cache corresponding to the fluid state
156 * \param problem The problem to be solved
157 * \param element An element which contains part of the control volume
158 * \param scv The sub-control volume
159 */
160 template<class ElemSol, class Problem, class Element, class Scv>
161 void updateInterfacialArea(const ElemSol& elemSol,
162 const FluidState& fluidState,
163 const ParameterCache& paramCache,
164 const Problem& problem,
165 const Element& element,
166 const Scv& scv)
167 {
168 const Scalar pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
169 const Scalar Sw = fluidState.saturation(phase0Idx);
170
171 const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
172 const auto areaWN = fluidMatrixInteraction.wettingNonwettingInterface().area(Sw, pc);
173 const auto areaNS = fluidMatrixInteraction.nonwettingSolidInterface().area(Sw, pc);
174 interfacialArea_[phase0Idx][phase1Idx] = areaWN;
175 interfacialArea_[phase1Idx][phase0Idx] = interfacialArea_[phase0Idx][phase1Idx];
176 interfacialArea_[phase0Idx][phase0Idx] = 0.0;
177
178 // Switch for using a a_{wn} relations that has some "maximum capillary pressure" as parameter
179 // That value is obtained by regularization of the pc(Sw) function.
180 static const bool computeAwsFromAnsAndPcMax = getParam<bool>("SpatialParams.ComputeAwsFromAnsAndPcMax", true);
181 if (computeAwsFromAnsAndPcMax)
182 {
183 // I know the solid surface from the pore network. But it is more consistent to use the fit value.
184 const Scalar pcMax = fluidMatrixInteraction.wettingNonwettingInterface().basicParams().pcMax();
185 const auto solidSurface = fluidMatrixInteraction.nonwettingSolidInterface().area(/*Sw=*/0., pcMax);
186 interfacialArea_[phase0Idx][sPhaseIdx] = solidSurface - areaNS;
187 }
188 else
189 interfacialArea_[phase0Idx][sPhaseIdx] = fluidMatrixInteraction.wettingSolidInterface().area(Sw, pc);
190
191 interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
192 interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.0;
193
194 interfacialArea_[phase1Idx][sPhaseIdx] = areaNS;
195 interfacialArea_[sPhaseIdx][phase1Idx] = interfacialArea_[phase1Idx][sPhaseIdx];
196 interfacialArea_[phase1Idx][phase1Idx] = 0.0;
197 }
198
199 /*!
200 * \brief The specific interfacial area between two fluid phases [m^2 / m^3]
201 * \note This is _only_ required by the kinetic mass/energy modules
202 */
203 const Scalar interfacialArea(const unsigned int phaseIIdx, const unsigned int phaseJIdx) const
204 {
205 // there is no interfacial area between a phase and itself
206 assert(phaseIIdx not_eq phaseJIdx);
207 return interfacialArea_[phaseIIdx][phaseJIdx];
208 }
209
210 //! access function Reynolds Number
211 const Scalar reynoldsNumber(const unsigned int phaseIdx) const
212 { return reynoldsNumber_[phaseIdx]; }
213
214 //! access function Prandtl Number
215 const Scalar prandtlNumber(const unsigned int phaseIdx) const
216 { return prandtlNumber_[phaseIdx]; }
217
218 //! access function Nusselt Number
219 const Scalar nusseltNumber(const unsigned int phaseIdx) const
220 { return nusseltNumber_[phaseIdx]; }
221
222 //! access function characteristic length
223 const Scalar characteristicLength() const
224 { return characteristicLength_; }
225
226 //! access function pre factor energy transfer
227 const Scalar factorEnergyTransfer() const
228 { return factorEnergyTransfer_; }
229
230 private:
231 //! dimensionless numbers
232 NumFluidPhasesArray reynoldsNumber_;
233 NumFluidPhasesArray prandtlNumber_;
234 NumFluidPhasesArray nusseltNumber_;
235
236 Scalar characteristicLength_;
237 Scalar factorEnergyTransfer_;
238 InterfacialAreasArray interfacialArea_;
239 };
240
241 /////////////////////////////////////////////////////////////////////////////////////////////////////////////
242 // specialization for the case of NO kinetic mass but kinetic energy transfer of a fluid mixture and a solid
243 /////////////////////////////////////////////////////////////////////////////////////////////////////////////
244 template<class Traits, class EquilibriumVolumeVariables>
245 1070190 class NonEquilibriumVolumeVariablesImplementation<Traits,
246 EquilibriumVolumeVariables,
247 false/*chemicalNonEquilibrium?*/,
248 true/*thermalNonEquilibrium?*/,
249 1>
250 : public EquilibriumVolumeVariables
251 {
252 using ParentType = EquilibriumVolumeVariables;
253 using ParameterCache = typename Traits::FluidSystem::ParameterCache;
254 using Scalar = typename Traits::PrimaryVariables::value_type;
255
256 using ModelTraits = typename Traits::ModelTraits;
257 using FS = typename Traits::FluidSystem;
258 static constexpr auto numEnergyEqFluid = ModelTraits::numEnergyEqFluid();
259 static constexpr auto numEnergyEqSolid = ModelTraits::numEnergyEqSolid();
260
261 using DimLessNum = DimensionlessNumbers<Scalar>;
262
263 using NumFluidPhasesArray = std::array<Scalar, ModelTraits::numFluidPhases()>;
264
265 public:
266 using Indices = typename ModelTraits::Indices;
267 using FluidState = typename Traits::FluidState;
268
269 /*!
270 * \brief Update all quantities for a given control volume
271 *
272 * \param elemSol A vector containing all primary variables connected to the element
273 * \param problem The object specifying the problem which ought to
274 * be simulated
275 * \param element An element which contains part of the control volume
276 * \param scv The sub-control volume
277 */
278 template<class ElemSol, class Problem, class Element, class Scv>
279 1391823 void update(const ElemSol &elemSol,
280 const Problem &problem,
281 const Element &element,
282 const Scv& scv)
283 {
284 // Update parent type (also completes the fluid state)
285 1391823 ParentType::update(elemSol, problem, element, scv);
286
287 ParameterCache paramCache;
288 2783646 paramCache.updateAll(this->fluidState());
289 //only update of DimLessNumbers is necessary here, as interfacial area
290 // is easy due to only one fluid with a solid and is directly computed in localresidual
291 1391823 updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
292 2783646 updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
293 1391823 }
294
295 /*!
296 * \brief Updates dimensionless numbers
297 *
298 * \param elemSol A vector containing all primary variables connected to the element
299 * \param fluidState Container for all the secondary variables concerning the fluids
300 * \param paramCache The parameter cache corresponding to the fluid state
301 * \param problem The problem to be solved
302 * \param element An element which contains part of the control volume
303 * \param scv The sub-control volume
304 */
305 template<class ElemSol, class Problem, class Element, class Scv>
306 1391823 void updateDimLessNumbers(const ElemSol& elemSol,
307 const FluidState& fluidState,
308 const ParameterCache& paramCache,
309 const Problem& problem,
310 const Element& element,
311 const Scv& scv)
312 {
313 1391823 const auto& spatialParams = problem.spatialParams();
314 1391823 factorEnergyTransfer_ = spatialParams.factorEnergyTransfer(element, scv, elemSol);
315 1391823 characteristicLength_ = spatialParams.characteristicLength(element, scv, elemSol);
316
317 // set the dimensionless numbers and obtain respective quantities
318 1391823 const unsigned int vIdxGlobal = scv.dofIndex();
319
2/2
✓ Branch 0 taken 1391823 times.
✓ Branch 1 taken 1391823 times.
2783646 for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
320 {
321
0/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
2783646 const auto darcyMagVelocity = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, vIdxGlobal);
322
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
1391823 const auto dynamicViscosity = fluidState.viscosity(phaseIdx);
323
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
1391823 const auto density = fluidState.density(phaseIdx);
324 1391823 const auto kinematicViscosity = dynamicViscosity/density;
325
326 using FluidSystem = typename Traits::FluidSystem;
327
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
1391823 const auto heatCapacity = FluidSystem::heatCapacity(fluidState, paramCache, phaseIdx);
328
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
1391823 const auto thermalConductivity = FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);
329 1391823 const auto porosity = this->porosity();
330
331 2783646 reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_,kinematicViscosity);
332 2783646 prandtlNumber_[phaseIdx] = DimLessNum::prandtlNumber(dynamicViscosity, heatCapacity, thermalConductivity);
333 1391823 nusseltNumber_[phaseIdx] = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
334 2783646 prandtlNumber_[phaseIdx],
335 porosity,
336 ModelTraits::nusseltFormulation());
337 }
338 1391823 }
339
340
341 /*!
342 * \brief Updates the volume specific interfacial area [m^2 / m^3] between the solid and the fluid phase.
343 *
344 * \param elemSol A vector containing all primary variables connected to the element
345 * \param fluidState Container for all the secondary variables concerning the fluids
346 * \param paramCache The parameter cache corresponding to the fluid state
347 * \param problem The problem to be solved
348 * \param element An element which contains part of the control volume
349 * \param scv The sub-control volume
350 */
351 template<class ElemSol, class Problem, class Element, class Scv>
352 void updateInterfacialArea(const ElemSol& elemSol,
353 const FluidState& fluidState,
354 const ParameterCache& paramCache,
355 const Problem& problem,
356 const Element& element,
357 const Scv& scv)
358 {
359 using FluidSolidInterfacialAreaFormulation = typename Problem::SpatialParams::FluidSolidInterfacialAreaFormulation;
360 4175469 interfacialArea_ = FluidSolidInterfacialAreaFormulation::fluidSolidInterfacialArea(this->porosity(), characteristicLength());
361 }
362
363 //! access function Reynolds Number
364 const Scalar reynoldsNumber(const unsigned int phaseIdx) const
365 91600 { return reynoldsNumber_[phaseIdx]; }
366
367 //! access function Prandtl Number
368 const Scalar prandtlNumber(const unsigned int phaseIdx) const
369 91600 { return prandtlNumber_[phaseIdx]; }
370
371 //! access function Nusselt Number
372 const Scalar nusseltNumber(const unsigned int phaseIdx) const
373
0/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
5010400 { return nusseltNumber_[phaseIdx]; }
374
375 //! access function characteristic length
376 const Scalar characteristicLength() const
377 { return characteristicLength_; }
378
379 //! access function pre factor energy transfer
380 const Scalar factorEnergyTransfer() const
381 { return factorEnergyTransfer_; }
382
383 const Scalar fluidSolidInterfacialArea() const
384 {return interfacialArea_;}
385
386 private:
387 //! dimensionless numbers
388 NumFluidPhasesArray reynoldsNumber_;
389 NumFluidPhasesArray prandtlNumber_;
390 NumFluidPhasesArray nusseltNumber_;
391
392 Scalar characteristicLength_;
393 Scalar factorEnergyTransfer_;
394 Scalar interfacialArea_ ;
395 };
396
397 ////////////////////////////////////////////////////////////////////////////////////////////////////
398 // specialization for the case of (only) kinetic mass transfer. Be careful, we do not test this!
399 ////////////////////////////////////////////////////////////////////////////////////////////////////
400 template<class Traits, class EquilibriumVolumeVariables>
401
2/3
✓ Branch 5 taken 1600 times.
✓ Branch 6 taken 6400 times.
✗ Branch 7 not taken.
173272 class NonEquilibriumVolumeVariablesImplementation<Traits,
402 EquilibriumVolumeVariables,
403 true/*chemicalNonEquilibrium?*/,
404 false/*thermalNonEquilibrium?*/,
405 0>
406 : public EquilibriumVolumeVariables
407 {
408 using ParentType = EquilibriumVolumeVariables;
409 using FluidState = typename Traits::FluidState;
410 using FS = typename Traits::FluidSystem;
411 using ParameterCache = typename Traits::FluidSystem::ParameterCache;
412 using Scalar = typename Traits::PrimaryVariables::value_type;
413
414 using ModelTraits = typename Traits::ModelTraits;
415
416 static constexpr auto phase0Idx = FS::phase0Idx;
417 static constexpr auto phase1Idx = FS::phase1Idx;
418 static constexpr auto wCompIdx = FS::comp0Idx;
419 static constexpr auto nCompIdx = FS::comp1Idx;
420
421 using DimLessNum = DimensionlessNumbers<Scalar>;
422
423 using NumFluidPhasesArray = std::array<Scalar, ModelTraits::numFluidPhases()>;
424
425 public:
426 using Indices = typename ModelTraits::Indices;
427 /*!
428 * \brief Update all quantities for a given control volume
429 *
430 * \param elemSol A vector containing all primary variables connected to the element
431 * \param problem The object specifying the problem which ought to
432 * be simulated
433 * \param element An element which contains part of the control volume
434 * \param scv The sub-control volume
435 */
436 template<class ElemSol, class Problem, class Element, class Scv>
437 394268 void update(const ElemSol &elemSol,
438 const Problem &problem,
439 const Element &element,
440 const Scv& scv)
441 {
442 // Update parent type (also completes the fluid state)
443 394268 ParentType::update(elemSol, problem, element, scv);
444
445 ParameterCache paramCache;
446 788536 paramCache.updateAll(this->fluidState());
447 394268 updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
448 394268 updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
449 394268 }
450
451 /*!
452 * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
453 *
454 * \param elemSol A vector containing all primary variables connected to the element
455 * \param fluidState Container for all the secondary variables concerning the fluids
456 * \param paramCache The parameter cache corresponding to the fluid state
457 * \param problem The problem to be solved
458 * \param element An element which contains part of the control volume
459 * \param scv The sub-control volume
460 */
461 template<class ElemSol, class Problem, class Element, class Scv>
462 394268 void updateDimLessNumbers(const ElemSol& elemSol,
463 const FluidState& fluidState,
464 const ParameterCache& paramCache,
465 const Problem& problem,
466 const Element& element,
467 const Scv& scv)
468 {
469 394268 const auto& spatialParams = problem.spatialParams();
470 394268 factorMassTransfer_ = spatialParams.factorMassTransfer(element, scv, elemSol);
471 394268 characteristicLength_ = spatialParams.characteristicLength(element, scv, elemSol);
472
473 // set the dimensionless numbers and obtain respective quantities.
474 394268 const unsigned int vIdxGlobal = scv.dofIndex();
475
2/2
✓ Branch 0 taken 788536 times.
✓ Branch 1 taken 394268 times.
1182804 for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
476 {
477 1577072 const auto darcyMagVelocity = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, vIdxGlobal);
478 788536 const auto dynamicViscosity = fluidState.viscosity(phaseIdx);
479 788536 const auto density = fluidState.density(phaseIdx);
480 788536 const auto kinematicViscosity = dynamicViscosity/density;
481
482 // diffusion coefficient of nonwetting component in wetting phase
483 using FluidSystem = typename Traits::FluidSystem;
484 788536 const auto diffCoeff = FluidSystem::binaryDiffusionCoefficient(fluidState,
485 paramCache,
486 phaseIdx,
487 wCompIdx,
488 nCompIdx);
489
490 1577072 reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_, kinematicViscosity);
491 1577072 schmidtNumber_[phaseIdx] = DimLessNum::schmidtNumber(dynamicViscosity, density, diffCoeff);
492 788536 sherwoodNumber_[phaseIdx] = DimLessNum::sherwoodNumber(reynoldsNumber_[phaseIdx],
493 1577072 schmidtNumber_[phaseIdx],
494 ModelTraits::sherwoodFormulation());
495 }
496 394268 }
497
498 /*!
499 * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
500 *
501 * \param elemSol A vector containing all primary variables connected to the element
502 * \param fluidState Container for all the secondary variables concerning the fluids
503 * \param paramCache The parameter cache corresponding to the fluid state
504 * \param problem The problem to be solved
505 * \param element An element which contains part of the control volume
506 * \param scv The sub-control volume
507 */
508 template<class ElemSol, class Problem, class Element, class Scv>
509 394268 void updateInterfacialArea(const ElemSol& elemSol,
510 const FluidState& fluidState,
511 const ParameterCache& paramCache,
512 const Problem& problem,
513 const Element& element,
514 const Scv& scv)
515 {
516 788536 const auto Sw = fluidState.saturation(phase0Idx) ;
517 1182804 const auto pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
518
519 // when we only consider chemical non-equilibrium there is only mass transfer between
520 // the fluid phases, so in 2p only interfacial area between wetting and non-wetting
521 1182804 const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
522 788536 interfacialArea_ = fluidMatrixInteraction.wettingNonwettingInterface().area(Sw, pc);
523 394268 }
524
525 /*!
526 * \brief The specific interfacial area between two fluid phases [m^2 / m^3]
527 */
528 const Scalar interfacialArea(const unsigned int phaseIIdx, const unsigned int phaseJIdx) const
529 {
530 // so far there is only a model for kinetic mass transfer between fluid phases
531 assert( (phaseIIdx == phase1Idx && phaseJIdx == phase0Idx)
532 || (phaseIIdx == phase0Idx && phaseJIdx == phase1Idx) );
533 628800 return interfacialArea_;
534 }
535
536 //! access function Reynolds Number
537 const Scalar reynoldsNumber(const unsigned int phaseIdx) const
538 24000 { return reynoldsNumber_[phaseIdx]; }
539
540 //! access function Schmidt Number
541 const Scalar schmidtNumber(const unsigned int phaseIdx) const
542 { return schmidtNumber_[phaseIdx]; }
543
544 //! access function Sherwood Number
545 const Scalar sherwoodNumber(const unsigned int phaseIdx) const
546 2515200 { return sherwoodNumber_[phaseIdx]; }
547
548 //! access function characteristic length
549 const Scalar characteristicLength() const
550 { return characteristicLength_; }
551
552 //! access function pre factor mass transfer
553 const Scalar factorMassTransfer() const
554 { return factorMassTransfer_; }
555
556 private:
557 Scalar characteristicLength_;
558 Scalar factorMassTransfer_;
559 Scalar interfacialArea_ ;
560 NumFluidPhasesArray sherwoodNumber_;
561 NumFluidPhasesArray schmidtNumber_;
562 NumFluidPhasesArray reynoldsNumber_;
563 };
564
565 // Specialization where everything is in non-equilibrium.
566 template<class Traits, class EquilibriumVolumeVariables>
567 191720 class NonEquilibriumVolumeVariablesImplementation<Traits,
568 EquilibriumVolumeVariables,
569 true/*chemicalNonEquilibrium?*/,
570 true/*thermalNonEquilibrium?*/,
571 2>
572 : public EquilibriumVolumeVariables
573 {
574 using ParentType = EquilibriumVolumeVariables;
575 using FluidState = typename Traits::FluidState;
576 using FS = typename Traits::FluidSystem;
577 using ParameterCache = typename Traits::FluidSystem::ParameterCache;
578 using Scalar = typename Traits::PrimaryVariables::value_type;
579
580 using ModelTraits = typename Traits::ModelTraits;
581 static constexpr auto numEnergyEqFluid = ModelTraits::numEnergyEqFluid();
582 static constexpr auto numEnergyEqSolid = ModelTraits::numEnergyEqSolid();
583
584 static constexpr auto phase0Idx = FS::phase0Idx;
585 static constexpr auto phase1Idx = FS::phase1Idx;
586 static constexpr auto sPhaseIdx = FS::numPhases;
587 static constexpr auto wCompIdx = FS::comp0Idx;
588 static constexpr auto nCompIdx = FS::comp1Idx;
589
590 using DimLessNum = DimensionlessNumbers<Scalar>;
591 static_assert((numEnergyEqFluid > 1), "This model only deals with energy transfer between two fluids and one solid phase");
592
593 using NumFluidPhasesArray = std::array<Scalar, ModelTraits::numFluidPhases()>;
594 using InterfacialAreasArray = std::array<std::array<Scalar, ModelTraits::numFluidPhases()+numEnergyEqSolid>,
595 ModelTraits::numFluidPhases()+numEnergyEqSolid>;
596
597 public:
598 using Indices = typename ModelTraits::Indices;
599
600 /*!
601 * \brief Update all quantities for a given control volume
602 *
603 * \param elemSol A vector containing all primary variables connected to the element
604 * \param problem The object specifying the problem which ought to
605 * be simulated
606 * \param element An element which contains part of the control volume
607 * \param scv The sub-control volume
608 */
609 template<class ElemSol, class Problem, class Element, class Scv>
610 786240 void update(const ElemSol &elemSol,
611 const Problem &problem,
612 const Element &element,
613 const Scv& scv)
614 {
615 // Update parent type (also completes the fluid state)
616 786240 ParentType::update(elemSol, problem, element, scv);
617
618 ParameterCache paramCache;
619 1572480 paramCache.updateAll(this->fluidState());
620 786240 updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
621 786240 updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
622 786240 }
623
624 /*!
625 * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
626 *
627 * \param elemSol A vector containing all primary variables connected to the element
628 * \param fluidState Container for all the secondary variables concerning the fluids
629 * \param paramCache The parameter cache corresponding to the fluid state
630 * \param problem The problem to be solved
631 * \param element An element which contains part of the control volume
632 * \param scv The sub-control volume
633 */
634 template<class ElemSol, class Problem, class Element, class Scv>
635 786240 void updateDimLessNumbers(const ElemSol& elemSol,
636 const FluidState& fluidState,
637 const ParameterCache& paramCache,
638 const Problem& problem,
639 const Element& element,
640 const Scv& scv)
641 {
642 786240 const auto& spatialParams = problem.spatialParams();
643 786240 factorMassTransfer_ = spatialParams.factorMassTransfer(element, scv, elemSol);
644 786240 factorEnergyTransfer_ = spatialParams.factorEnergyTransfer(element, scv, elemSol);
645 786240 characteristicLength_ = spatialParams.characteristicLength(element, scv, elemSol);
646
647 786240 const auto vIdxGlobal = scv.dofIndex();
648 using FluidSystem = typename Traits::FluidSystem;
649
2/2
✓ Branch 0 taken 1572480 times.
✓ Branch 1 taken 786240 times.
2358720 for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
650 {
651 3144960 const auto darcyMagVelocity = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, vIdxGlobal);
652 1572480 const auto dynamicViscosity = fluidState.viscosity(phaseIdx);
653 1572480 const auto density = fluidState.density(phaseIdx);
654 1572480 const auto kinematicViscosity = dynamicViscosity/density;
655 1572480 const auto heatCapacity = FluidSystem::heatCapacity(fluidState, paramCache, phaseIdx);
656 1572480 const auto thermalConductivity = FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);
657
658 // diffusion coefficient of nonwetting component in wetting phase
659 1572480 const auto porosity = this->porosity();
660 1572480 const auto diffCoeff = FluidSystem::binaryDiffusionCoefficient(fluidState,
661 paramCache,
662 phaseIdx,
663 wCompIdx,
664 nCompIdx);
665
666 3144960 reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_, kinematicViscosity);
667 3144960 prandtlNumber_[phaseIdx] = DimLessNum::prandtlNumber(dynamicViscosity, heatCapacity, thermalConductivity);
668 3144960 schmidtNumber_[phaseIdx] = DimLessNum::schmidtNumber(dynamicViscosity, density, diffCoeff);
669 1572480 nusseltNumber_[phaseIdx] = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
670 3144960 prandtlNumber_[phaseIdx],
671 porosity,
672 ModelTraits::nusseltFormulation());
673 // If Diffusion is not enabled, Sherwood is divided by zero
674 1572480 sherwoodNumber_[phaseIdx] = DimLessNum::sherwoodNumber(reynoldsNumber_[phaseIdx],
675 3144960 schmidtNumber_[phaseIdx],
676 ModelTraits::sherwoodFormulation());
677 }
678 786240 }
679
680 /*!
681 * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
682 *
683 * \param elemSol A vector containing all primary variables connected to the element
684 * \param fluidState Container for all the secondary variables concerning the fluids
685 * \param paramCache The parameter cache corresponding to the fluid state
686 * \param problem The problem to be solved
687 * \param element An element which contains part of the control volume
688 * \param scv The sub-control volume
689 */
690 template<class ElemSol, class Problem, class Element, class Scv>
691 786240 void updateInterfacialArea(const ElemSol& elemSol,
692 const FluidState& fluidState,
693 const ParameterCache& paramCache,
694 const Problem& problem,
695 const Element& element,
696 const Scv& scv)
697 {
698 1572480 const Scalar pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
699 786240 const Scalar Sw = fluidState.saturation(phase0Idx);
700
701 1572480 const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
702
703 786240 const auto awn = fluidMatrixInteraction.wettingNonwettingInterface().area(Sw, pc);
704 1572480 interfacialArea_[phase0Idx][phase1Idx] = awn;
705 3144960 interfacialArea_[phase1Idx][phase0Idx] = interfacialArea_[phase0Idx][phase1Idx];
706 1572480 interfacialArea_[phase0Idx][phase0Idx] = 0.;
707
708 786240 const auto areaNS = fluidMatrixInteraction.nonwettingSolidInterface().area(Sw, pc);
709
710 // Switch for using a a_{wn} relations that has some "maximum capillary pressure" as parameter.
711 // That value is obtained by regularization of the pc(Sw) function.
712
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 786239 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
786240 static const bool computeAwsFromAnsAndPcMax = getParam<bool>("SpatialParams.ComputeAwsFromAnsAndPcMax", true);
713
1/2
✓ Branch 0 taken 786240 times.
✗ Branch 1 not taken.
786240 if (computeAwsFromAnsAndPcMax)
714 {
715 // I know the solid surface from the pore network. But it is more consistent to use the fit value.
716 786240 const Scalar pcMax = fluidMatrixInteraction.wettingNonwettingInterface().basicParams().pcMax();
717 786240 const auto solidSurface = fluidMatrixInteraction.nonwettingSolidInterface().area(/*Sw=*/0., pcMax);
718 2358720 interfacialArea_[phase0Idx][sPhaseIdx] = solidSurface - areaNS;
719 }
720 else
721 interfacialArea_[phase0Idx][sPhaseIdx] = fluidMatrixInteraction.wettingSolidInterface().area(Sw, pc);
722
723 3144960 interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
724 1572480 interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
725
726 1572480 interfacialArea_[phase1Idx][sPhaseIdx] = areaNS;
727 3144960 interfacialArea_[sPhaseIdx][phase1Idx] = interfacialArea_[phase1Idx][sPhaseIdx];
728 1572480 interfacialArea_[phase1Idx][phase1Idx] = 0.;
729 786240 }
730
731 /*!
732 * \brief The specific interfacial area between two fluid phases [m^2 / m^3]
733 * \note This is _only_ required by the kinetic mass/energy modules
734 */
735 const Scalar interfacialArea(const unsigned int phaseIIdx, const unsigned int phaseJIdx) const
736 {
737 // there is no interfacial area between a phase and itself
738 assert(phaseIIdx not_eq phaseJIdx);
739 28344960 return interfacialArea_[phaseIIdx][phaseJIdx];
740 }
741
742 //! access function Reynolds Number
743 const Scalar reynoldsNumber(const unsigned int phaseIdx) const
744 80640 { return reynoldsNumber_[phaseIdx]; }
745
746 //! access function Prandtl Number
747 const Scalar prandtlNumber(const unsigned int phaseIdx) const
748 80640 { return prandtlNumber_[phaseIdx]; }
749
750 //! access function Nusselt Number
751 const Scalar nusseltNumber(const unsigned int phaseIdx) const
752
4/8
✓ Branch 0 taken 2362080 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2362080 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2362080 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2362080 times.
✗ Branch 7 not taken.
18896640 { return nusseltNumber_[phaseIdx]; }
753
754 //! access function Schmidt Number
755 const Scalar schmidtNumber(const unsigned int phaseIdx) const
756 { return schmidtNumber_[phaseIdx]; }
757
758 //! access function Sherwood Number
759 const Scalar sherwoodNumber(const unsigned int phaseIdx) const
760 9448320 { return sherwoodNumber_[phaseIdx]; }
761
762 //! access function characteristic length
763 const Scalar characteristicLength() const
764 { return characteristicLength_; }
765
766 //! access function pre factor energy transfer
767 const Scalar factorEnergyTransfer() const
768 { return factorEnergyTransfer_; }
769
770 //! access function pre factor mass transfer
771 const Scalar factorMassTransfer() const
772 { return factorMassTransfer_; }
773
774 private:
775 //! dimensionless numbers
776 NumFluidPhasesArray reynoldsNumber_;
777 NumFluidPhasesArray prandtlNumber_;
778 NumFluidPhasesArray nusseltNumber_;
779 NumFluidPhasesArray schmidtNumber_;
780 NumFluidPhasesArray sherwoodNumber_;
781 Scalar characteristicLength_;
782 Scalar factorEnergyTransfer_;
783 Scalar factorMassTransfer_;
784 InterfacialAreasArray interfacialArea_;
785 };
786
787 } // namespace Dumux
788
789 #endif
790