volumevariables.hh 23.2 KB
Newer Older
1
2
3
4
5
6
7
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
9
10
11
12
13
14
15
16
17
18
19
20
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
21
 * \ingroup NIModel
22
23
24
 * \brief Base class for the model specific class which provides
 *        access to all volume averaged quantities.
 */
25

26
27
28
#ifndef DUMUX_ENERGY_VOLUME_VARIABLES_HH
#define DUMUX_ENERGY_VOLUME_VARIABLES_HH

29
#include <type_traits>
30
#include <dune/common/std/type_traits.hh>
31

32
#include <dumux/material/solidsystems/1csolid.hh>
33
34
35
36
#include <dumux/porousmediumflow/volumevariables.hh>

namespace Dumux {

37
38
#ifndef DOXYGEN
namespace Detail {
39
// helper structs and functions detecting if the user-defined spatial params class has user-specified functions
40
41
// for solidHeatCapacity, solidDensity, and solidThermalConductivity.

42
43
template <typename T, typename ...Ts>
using SolidHeatCapacityDetector = decltype(std::declval<T>().solidHeatCapacity(std::declval<Ts>()...));
44

45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
template<class T, typename ...Args>
static constexpr bool hasSolidHeatCapacity()
{ return Dune::Std::is_detected<SolidHeatCapacityDetector, T, Args...>::value; }

template <typename T, typename ...Ts>
using SolidDensityDetector = decltype(std::declval<T>().solidDensity(std::declval<Ts>()...));

template<class T, typename ...Args>
static constexpr bool hasSolidDensity()
{ return Dune::Std::is_detected<SolidDensityDetector, T, Args...>::value; }

template <typename T, typename ...Ts>
using SolidThermalConductivityDetector = decltype(std::declval<T>().solidThermalConductivity(std::declval<Ts>()...));

template<class T, typename ...Args>
static constexpr bool hasSolidThermalConductivity()
{ return Dune::Std::is_detected<SolidThermalConductivityDetector, T, Args...>::value; }
62
63
64
65
66
67
68
69
70
71
72

template<class SolidSystem>
struct isInertSolidPhase : public std::false_type {};

template<class Scalar, class Component>
struct isInertSolidPhase<SolidSystems::InertSolidPhase<Scalar, Component>> : public std::true_type {};

} // end namespace Detail
#endif


73
74
75
76
77
// forward declaration
template <class IsothermalTraits, class Impl, bool enableEnergyBalance>
class EnergyVolumeVariablesImplementation;

/*!
Katharina Heck's avatar
Katharina Heck committed
78
 * \ingroup NIModel
79
 * \brief Base class for the model specific class which provides
80
81
82
 *        access to all volume averaged quantities.
 *
 * The volume variables base class is specialized for isothermal and non-isothermal models.
83
84
85
86
87
 */
template<class IsothermalTraits, class Impl>
using EnergyVolumeVariables = EnergyVolumeVariablesImplementation<IsothermalTraits,Impl, IsothermalTraits::ModelTraits::enableEnergyBalance()>;

/*!
Katharina Heck's avatar
Katharina Heck committed
88
 * \ingroup NIModel
89
90
91
92
93
94
95
96
97
98
99
100
101
102
 * \brief The isothermal base class
 */
template<class IsothermalTraits, class Impl>
class EnergyVolumeVariablesImplementation<IsothermalTraits, Impl, false>
{
    using Scalar = typename IsothermalTraits::PrimaryVariables::value_type;

public:
    using FluidState = typename IsothermalTraits::FluidState;
    using SolidState = typename IsothermalTraits::SolidState;
    using FluidSystem = typename IsothermalTraits::FluidSystem;

    //! The temperature is obtained from the problem as a constant for isothermal models
    template<class ElemSol, class Problem, class Element, class Scv>
Dennis Gläser's avatar
Dennis Gläser committed
103
    void updateTemperature(const ElemSol& elemSol,
104
                           const Problem& problem,
Dennis Gläser's avatar
Dennis Gläser committed
105
106
                           const Element& element,
                           const Scv& scv,
107
                           FluidState& fluidState,
Dennis Gläser's avatar
Dennis Gläser committed
108
                           SolidState& solidState)
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
    {
        // retrieve temperature from solution vector, all phases have the same temperature
        Scalar T = problem.temperatureAtPos(scv.dofPosition());
        for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
        {
            fluidState.setTemperature(phaseIdx, T);
        }
        solidState.setTemperature(T);
    }

    template<class ElemSol, class Problem, class Element, class Scv>
    void updateSolidEnergyParams(const ElemSol &elemSol,
                                 const Problem& problem,
                                 const Element &element,
                                 const Scv &scv,
                                 SolidState & solidState)
    {}

    //! The phase enthalpy is zero for isothermal models
    //! This is needed for completing the fluid state
    template<class FluidState, class ParameterCache>
    static Scalar enthalpy(const FluidState& fluidState,
                           const ParameterCache& paramCache,
                           const int phaseIdx)
    {
        return 0;
    }

137
138
139
140
    //! The effective thermal conductivity is zero for isothermal models
    void updateEffectiveThermalConductivity()
    {}

141
142
143
};

//! The non-isothermal implicit volume variables base class
144
145
template<class Traits, class Impl>
class EnergyVolumeVariablesImplementation<Traits, Impl, true>
146
{
147
148
149
150
    using Scalar = typename Traits::PrimaryVariables::value_type;
    using Idx = typename Traits::ModelTraits::Indices;
    using ParentType = PorousMediumFlowVolumeVariables<Traits>;
    using EffCondModel = typename Traits::EffectiveThermalConductivityModel;
151

152
153
154
155
156
    static constexpr int temperatureIdx = Idx::temperatureIdx;
    static constexpr int numEnergyEq = Traits::ModelTraits::numEnergyEq();

    static constexpr bool fullThermalEquilibrium = (numEnergyEq == 1);
    static constexpr bool fluidThermalEquilibrium = (numEnergyEq == 2);
157
158
159

public:
    // export the fluidstate
160
    using FluidState = typename Traits::FluidState;
161
    //! export the underlying fluid system
162
    using FluidSystem = typename Traits::FluidSystem;
163
164
165
166
167
    //! Export the indices
    using Indices = Idx;
    // export the solidstate
    using SolidState = typename Traits::SolidState;
    //! export the underlying solid system
168
    using SolidSystem = typename Traits::SolidSystem;
169
170
171

    //! The temperature is obtained from the problem as a constant for isothermal models
    template<class ElemSol, class Problem, class Element, class Scv>
Dennis Gläser's avatar
Dennis Gläser committed
172
    void updateTemperature(const ElemSol& elemSol,
173
                           const Problem& problem,
Dennis Gläser's avatar
Dennis Gläser committed
174
175
                           const Element& element,
                           const Scv& scv,
176
                           FluidState& fluidState,
Dennis Gläser's avatar
Dennis Gläser committed
177
                           SolidState& solidState)
178
    {
179
        if constexpr (fullThermalEquilibrium)
180
181
        {
            // retrieve temperature from solution vector, all phases have the same temperature
182
            const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx];
183
184
            for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
            {
Dennis Gläser's avatar
Dennis Gläser committed
185
                fluidState.setTemperature(phaseIdx, T);
186
187
188
189
190
191
            }
            solidState.setTemperature(T);
        }

        else
        {
192
193
            // this means we have 1 temp for fluid phase, one for solid
            if constexpr (fluidThermalEquilibrium)
194
            {
195
                const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx];
196
197
198
199
                for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
                {
                    fluidState.setTemperature(phaseIdx, T);
                }
200
            }
201
            // this is for numEnergyEqFluid > 1
202
203
204
205
206
207
208
209
            else
            {
                for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
                {
                    // retrieve temperatures from solution vector, phases might have different temperature
                    const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx + phaseIdx];
                    fluidState.setTemperature(phaseIdx, T);
                }
210
            }
211
            const Scalar solidTemperature = elemSol[scv.localDofIndex()][temperatureIdx+numEnergyEq-1];
212
213
214
215
216
217
218
219
220
221
222
            solidState.setTemperature(solidTemperature);
        }
    }

    template<class ElemSol, class Problem, class Element, class Scv>
    void updateSolidEnergyParams(const ElemSol &elemSol,
                                 const Problem& problem,
                                 const Element &element,
                                 const Scv &scv,
                                 SolidState & solidState)
    {
223
224
        Scalar cs = solidHeatCapacity_(elemSol, problem, element, scv, solidState);
        solidState.setHeatCapacity(cs);
225

226
227
        Scalar rhos = solidDensity_(elemSol, problem, element, scv, solidState);
        solidState.setDensity(rhos);
228

229
230
        Scalar lambdas = solidThermalConductivity_(elemSol, problem, element, scv, solidState);
        solidState.setThermalConductivity(lambdas);
231
232
    }

233
234
235
    // updates the effective thermal conductivity
    void updateEffectiveThermalConductivity()
    {
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
        if constexpr (fullThermalEquilibrium)
        {
            // Full Thermal Equilibirum: One effectiveThermalConductivity value for all phases (solid & fluid).
            lambdaEff_[0] = EffCondModel::effectiveThermalConductivity(asImp_());
        }
        else if constexpr (fluidThermalEquilibrium)
        {
            // Fluid Thermal Equilibrium (Partial Nonequilibrium): One effectiveThermalConductivity for the fluids, one for the solids.
            Scalar fluidLambda = 0.0;
            for (int phaseIdx = 0; phaseIdx < FluidSystem::numPhases; phaseIdx++)
                fluidLambda += fluidThermalConductivity(phaseIdx) * asImp_().saturation(phaseIdx) * asImp_().porosity();

            lambdaEff_[0] = fluidLambda;
            lambdaEff_[numEnergyEq-1] = solidThermalConductivity() * (1.0 - asImp_().porosity());
        }
        else
        {
            // Full Thermal Nonequilibrium: One effectiveThermal Conductivity per phase (solid & fluid).
            for (int phaseIdx = 0; phaseIdx < FluidSystem::numPhases; phaseIdx++)
                lambdaEff_[phaseIdx] = fluidThermalConductivity(phaseIdx) * asImp_().saturation(phaseIdx) * asImp_().porosity();
            lambdaEff_[numEnergyEq-1] = solidThermalConductivity() * (1.0 - asImp_().porosity());
        }
258
259
    }

260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
    /*!
     * \brief Returns the total internal energy of a phase in the
     *        sub-control volume.
     *
     * \param phaseIdx The phase index
     */
    Scalar internalEnergy(const int phaseIdx) const
    { return asImp_().fluidState().internalEnergy(phaseIdx); }

    /*!
     * \brief Returns the total enthalpy of a phase in the sub-control
     *        volume.
     *
     * \param phaseIdx The phase index
     */
    Scalar enthalpy(const int phaseIdx) const
    { return asImp_().fluidState().enthalpy(phaseIdx); }

    /*!
     * \brief Returns the temperature in fluid / solid phase(s)
     *        the sub-control volume.
     */
    Scalar temperatureSolid() const
    { return asImp_().solidState().temperature(); }

285
286
287
288
289
290
291
292
293

    /*!
     * \brief Returns the temperature of a fluid phase assuming thermal nonequilibrium
     *        the sub-control volume.
     * \param phaseIdx The local index of the phases
     */
    Scalar temperatureFluid(const int phaseIdx) const
    { return asImp_().fluidState().temperature(phaseIdx); }

294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
    /*!
     * \brief Returns the total heat capacity \f$\mathrm{[J/(kg K)]}\f$ of the rock matrix in
     *        the sub-control volume.
     */
    Scalar solidHeatCapacity() const
    { return asImp_().solidState().heatCapacity(); }

    /*!
     * \brief Returns the mass density \f$\mathrm{[kg/m^3]}\f$ of the rock matrix in
     *        the sub-control volume.
     */
    Scalar solidDensity() const
    {  return  asImp_().solidState().density(); }

    /*!
309
310
     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
     *        of the solid phase in the sub-control volume.
311
312
313
314
315
     */
    Scalar solidThermalConductivity() const
    { return asImp_().solidState().thermalConductivity(); }

    /*!
316
317
     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
     *        of a fluid phase in the sub-control volume.
318
319
320
321
     */
    Scalar fluidThermalConductivity(const int phaseIdx) const
    { return FluidSystem::thermalConductivity(asImp_().fluidState(), phaseIdx); }

322
323
    /*!
     * \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ in
324
     *        the sub-control volume. Specific to equilibirum models (case fullThermalEquilibrium).
325
     */
326
327
    template< bool enable = fullThermalEquilibrium,
              std::enable_if_t<enable, int> = 0>
328
    Scalar effectiveThermalConductivity() const
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
    { return lambdaEff_[0]; }

    /*!
     * \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of the fluids in
     *        the sub-control volume. Specific to partially nonequilibrium models (case fluidThermalEquilibrium).
     */
    template< bool enable = fluidThermalEquilibrium,
              std::enable_if_t<enable, int> = 0>
    Scalar effectiveFluidThermalConductivity() const
    { return lambdaEff_[0]; }

    /*!
     * \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
     *        of the solid phase in the sub-control volume.
     *        Specific to partially nonequilibrium models (case fluidThermalEquilibrium)
     */
    template< bool enable = fluidThermalEquilibrium,
              std::enable_if_t<enable, int> = 0>
    Scalar effectiveSolidThermalConductivity() const
    { return lambdaEff_[numEnergyEq-1]; }

    /*!
     * \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
     *        per fluid phase in the sub-control volume.
     *        Specific to nonequilibrium models (case full non-equilibrium)
     */
    template< bool enable = (!fullThermalEquilibrium && !fluidThermalEquilibrium),
              std::enable_if_t<enable, int> = 0>
    Scalar effectivePhaseThermalConductivity(const int phaseIdx) const
    { return lambdaEff_[phaseIdx]; }
359

360
361
362
363
364
365
366
367
368
369
370
    //! The phase enthalpy is zero for isothermal models
    //! This is needed for completing the fluid state
    template<class ParameterCache>
    static Scalar enthalpy(const FluidState& fluidState,
                           const ParameterCache& paramCache,
                           const int phaseIdx)
    {
        return FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
    }

protected:
Katharina Heck's avatar
Katharina Heck committed
371
    const Impl &asImp_() const { return *static_cast<const Impl*>(this); }
372
373
    Impl &asImp_() { return *static_cast<Impl*>(this); }

374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
private:
    /*!
     * It has to be decided if the full solid system / solid state interface is used (general option, but more complicated),
     * or the simple nonisothermal spatial params interface (simpler but less general).
     * In the simple nonisothermal spatial params interface the functions solidHeatCapacity, solidDensity, and solidThermalConductivity
     * in the spatial params overwrite the parameters given in the solid system. This only makes sense in combination
     * with the simplest solid system InertSolidPhase, and can be used to quickly change parameters in certain domain regions.
     * For setups with more general solids with several components these functions should not exist. Instead, the solid system
     * determines the values for solidHeatCapacity, solidDensity, and solidThermalConductivity depending on the given composition.
     */

    /*!
     * \name Access functions for the solidsystem / solidstate interface
     */
    // \{

    /*!
391
392
     * \brief Gets the solid heat capacity in an scv.
     *
393
394
395
396
397
398
399
     * \param elemSol the element solution vector
     * \param problem the problem to solve
     * \param element the element (codim-0-entity) the scv belongs to
     * \param scv the sub control volume
     * \param solidState the solid state
     * \note this gets selected if the user uses the solidsystem / solidstate interface
     */
400
401
402
403
404
405
406
    template<class ElemSol, class Problem, class Element, class Scv,
             std::enable_if_t<!Detail::hasSolidHeatCapacity<typename Problem::SpatialParams, Element, Scv, ElemSol, SolidState>(), int> = 0>
    Scalar solidHeatCapacity_(const ElemSol& elemSol,
                              const Problem& problem,
                              const Element& element,
                              const Scv& scv,
                              const SolidState& solidState)
407
408
409
410
411
    {
        return SolidSystem::heatCapacity(solidState);
    }

    /*!
412
413
     * \brief Gets the solid density in an scv.
     *
414
415
416
417
418
419
420
     * \param elemSol the element solution vector
     * \param problem the problem to solve
     * \param element the element (codim-0-entity) the scv belongs to
     * \param scv the sub control volume
     * \param solidState the solid state
     * \note this gets selected if the user uses the solidsystem / solidstate interface
     */
421
422
423
424
425
426
427
    template<class ElemSol, class Problem, class Element, class Scv,
             std::enable_if_t<!Detail::hasSolidDensity<typename Problem::SpatialParams, Element, Scv, ElemSol, SolidState>(), int> = 0>
    Scalar solidDensity_(const ElemSol& elemSol,
                         const Problem& problem,
                         const Element& element,
                         const Scv& scv,
                         const SolidState& solidState)
428
429
430
431
432
    {
        return SolidSystem::density(solidState);
    }

    /*!
433
434
     * \brief Gets the solid's thermal conductivity in an scv.
     *
435
436
437
438
439
440
441
     * \param elemSol the element solution vector
     * \param problem the problem to solve
     * \param element the element (codim-0-entity) the scv belongs to
     * \param scv the sub control volume
     * \param solidState the solid state
     * \note this gets selected if the user uses the solidsystem / solidstate interface
     */
442
443
444
445
446
447
448
    template<class ElemSol, class Problem, class Element, class Scv,
             std::enable_if_t<!Detail::hasSolidThermalConductivity<typename Problem::SpatialParams, Element, Scv, ElemSol, SolidState>(), int> = 0>
    Scalar solidThermalConductivity_(const ElemSol& elemSol,
                                     const Problem& problem,
                                     const Element& element,
                                     const Scv& scv,
                                     const SolidState& solidState)
449
450
451
452
453
454
455
456
457
458
459
460
461
    {
        return SolidSystem::thermalConductivity(solidState);
    }

    // \}

    /*!
     * \name Access functions for the simple nonisothermal spatial params interface in
     *       combination with an InertSolidPhase as solid system
     */
    // \{

    /*!
462
463
     * \brief Gets the solid heat capacity in an scv.
     *
464
465
466
467
468
469
470
471
     * \param elemSol the element solution vector
     * \param problem the problem to solve
     * \param element the element (codim-0-entity) the scv belongs to
     * \param scv the sub control volume
     * \param solidState the solid state
     * \note this gets selected if the user uses the simple spatial params interface in
     *       combination with an InertSolidPhase as solid system
     */
472
473
474
475
476
477
478
    template<class ElemSol, class Problem, class Element, class Scv,
             std::enable_if_t<Detail::hasSolidHeatCapacity<typename Problem::SpatialParams, Element, Scv, ElemSol, SolidState>(), int> = 0>
    Scalar solidHeatCapacity_(const ElemSol& elemSol,
                              const Problem& problem,
                              const Element& element,
                              const Scv& scv,
                              const SolidState& solidState)
479
480
481
482
483
484
485
486
    {
        static_assert(Detail::isInertSolidPhase<SolidSystem>::value,
            "solidHeatCapacity can only be overwritten in the spatial params when the solid system is a simple InertSolidPhase\n"
            "If you select a proper solid system, the solid heat capacity will be computed as stated in the solid system!");
        return problem.spatialParams().solidHeatCapacity(element, scv, elemSol, solidState);
    }

    /*!
487
488
     * \brief Gets the solid density in an scv.
     *
489
490
491
492
493
494
495
496
     * \param elemSol the element solution vector
     * \param problem the problem to solve
     * \param element the element (codim-0-entity) the scv belongs to
     * \param scv the sub control volume
     * \param solidState the solid state
     * \note this gets selected if the user uses the simple spatial params interface in
     *       combination with an InertSolidPhase as solid system
     */
497
498
499
500
501
502
503
    template<class ElemSol, class Problem, class Element, class Scv,
             std::enable_if_t<Detail::hasSolidDensity<typename Problem::SpatialParams, Element, Scv, ElemSol, SolidState>(), int> = 0>
    Scalar solidDensity_(const ElemSol& elemSol,
                         const Problem& problem,
                         const Element& element,
                         const Scv& scv,
                         const SolidState& solidState)
504
505
506
507
508
509
510
511
    {
        static_assert(Detail::isInertSolidPhase<SolidSystem>::value,
            "solidDensity can only be overwritten in the spatial params when the solid system is a simple InertSolidPhase\n"
            "If you select a proper solid system, the solid density will be computed as stated in the solid system!");
        return problem.spatialParams().solidDensity(element, scv, elemSol, solidState);
    }

    /*!
512
513
     * \brief Gets the solid's heat capacity in an scv.
     *
514
515
516
517
518
519
520
521
     * \param elemSol the element solution vector
     * \param problem the problem to solve
     * \param element the element (codim-0-entity) the scv belongs to
     * \param scv the sub control volume
     * \param solidState the solid state
     * \note this gets selected if the user uses the simple spatial params interface in
     *       combination with an InertSolidPhase as solid system
     */
522
523
524
525
526
527
528
    template<class ElemSol, class Problem, class Element, class Scv,
             std::enable_if_t<Detail::hasSolidThermalConductivity<typename Problem::SpatialParams, Element, Scv, ElemSol, SolidState>(), int> = 0>
    Scalar solidThermalConductivity_(const ElemSol& elemSol,
                                     const Problem& problem,
                                     const Element& element,
                                     const Scv& scv,
                                     const SolidState& solidState)
529
530
531
532
533
534
535
    {
        static_assert(Detail::isInertSolidPhase<SolidSystem>::value,
            "solidThermalConductivity can only be overwritten in the spatial params when the solid system is a simple InertSolidPhase\n"
            "If you select a proper solid system, the solid thermal conductivity will be computed as stated in the solid system!");
        return problem.spatialParams().solidThermalConductivity(element, scv, elemSol, solidState);
    }

536
    std::array<Scalar, numEnergyEq> lambdaEff_;
537
538
    // \}

539
540
541
542
543
};

} // end namespace Dumux

#endif