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 |