volumevariables.hh 19.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
// -*- 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    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (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
 * \ingroup
 * \brief Base class for the model specific class which provides
 *        access to all volume averaged quantities.
 */
#ifndef DUMUX_ENERGY_VOLUME_VARIABLES_HH
#define DUMUX_ENERGY_VOLUME_VARIABLES_HH

28
29
30
#include <type_traits>

#include <dumux/material/solidsystems/inertsolidphase.hh>
31
32
33
34
#include <dumux/porousmediumflow/volumevariables.hh>

namespace Dumux {

35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#ifndef DOXYGEN
namespace Detail {
// helper struct detecting if the user-defined spatial params class has user-specified functions
// for solidHeatCapacity, solidDensity, and solidThermalConductivity.
// for g++ > 5.3, this can be replaced by a lambda
template<class E, class SCV, class ES, class SS>
struct hasSolidHeatCapacity
{
    template<class SpatialParams>
    auto operator()(const SpatialParams& a)
    -> decltype(a.solidHeatCapacity(std::declval<E>(), std::declval<SCV>(), std::declval<ES>(), std::declval<SS>()))
    {}
};

template<class E, class SCV, class ES, class SS>
struct hasSolidDensity
{
    template<class SpatialParams>
    auto operator()(const SpatialParams& a)
    -> decltype(a.solidDensity(std::declval<E>(), std::declval<SCV>(), std::declval<ES>(), std::declval<SS>()))
    {}
};

template<class E, class SCV, class ES, class SS>
struct hasSolidThermalConductivity
{
    template<class SpatialParams>
    auto operator()(const SpatialParams& a)
    -> decltype(a.solidThermalConductivity(std::declval<E>(), std::declval<SCV>(), std::declval<ES>(), std::declval<SS>()))
    {}
};

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


77
78
79
80
81
// forward declaration
template <class IsothermalTraits, class Impl, bool enableEnergyBalance>
class EnergyVolumeVariablesImplementation;

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

/*!
Katharina Heck's avatar
Katharina Heck committed
91
 * \ingroup NIModel
92
93
94
95
96
97
98
99
100
101
102
103
104
105
 * \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
106
    void updateTemperature(const ElemSol& elemSol,
107
                           const Problem& problem,
Dennis Gläser's avatar
Dennis Gläser committed
108
109
                           const Element& element,
                           const Scv& scv,
110
                           FluidState& fluidState,
Dennis Gläser's avatar
Dennis Gläser committed
111
                           SolidState& solidState)
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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
    {
        // 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;
    }

};

//! The non-isothermal implicit volume variables base class
template<class IsothermalTraits, class Impl>
class EnergyVolumeVariablesImplementation<IsothermalTraits, Impl, true>
{
    using Scalar = typename IsothermalTraits::PrimaryVariables::value_type;
    using Idx = typename IsothermalTraits::ModelTraits::Indices;
    using ParentType = PorousMediumFlowVolumeVariables<IsothermalTraits>;

    static const int temperatureIdx = Idx::temperatureIdx;
    static const int numEnergyEq = IsothermalTraits::ModelTraits::numEnergyEq();

public:
    // export the fluidstate
    using FluidState = typename IsothermalTraits::FluidState;
    using SolidState = typename IsothermalTraits::SolidState;
    //! export the underlying fluid system
    using FluidSystem = typename IsothermalTraits::FluidSystem;
    using SolidSystem = typename IsothermalTraits::SolidSystem;

    //! 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
163
    void updateTemperature(const ElemSol& elemSol,
164
                           const Problem& problem,
Dennis Gläser's avatar
Dennis Gläser committed
165
166
                           const Element& element,
                           const Scv& scv,
167
                           FluidState& fluidState,
Dennis Gläser's avatar
Dennis Gläser committed
168
                           SolidState& solidState)
169
170
171
172
    {
        if (numEnergyEq == 1)
        {
            // retrieve temperature from solution vector, all phases have the same temperature
173
            const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx];
174
175
            for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
            {
Dennis Gläser's avatar
Dennis Gläser committed
176
                fluidState.setTemperature(phaseIdx, T);
177
178
179
180
181
182
            }
            solidState.setTemperature(T);
        }

        else
        {
183
184
            //if numEnergyEq == 2 this means we have 1 temp for fluid phase, one for solid
            if (numEnergyEq == 2)
185
            {
186
187
188
189
190
191
192
193
194
195
196
197
                const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx];
                fluidState.setTemperature(T);
            }
            //this is for numEnergyEqFluid > 1
            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);
                }
198
            }
199
            const Scalar solidTemperature = elemSol[scv.localDofIndex()][temperatureIdx+numEnergyEq-1];
200
201
202
203
204
205
206
207
208
209
210
            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)
    {
211
212
        Scalar cs = solidHeatCapacity_(elemSol, problem, element, scv, solidState);
        solidState.setHeatCapacity(cs);
213

214
215
        Scalar rhos = solidDensity_(elemSol, problem, element, scv, solidState);
        solidState.setDensity(rhos);
216

217
218
        Scalar lambdas = solidThermalConductivity_(elemSol, problem, element, scv, solidState);
        solidState.setThermalConductivity(lambdas);
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
    }

    /*!
     * \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.
     * \param phaseIdx The local index of the phases
     */
    Scalar temperatureSolid() const
    { return asImp_().solidState().temperature(); }

247
248
249
250
251
252
253
254
255

    /*!
     * \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); }

256
257
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
285
286
287
288
289
290
291
292
293
    /*!
     * \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(); }

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

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

    //! 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
294
    const Impl &asImp_() const { return *static_cast<const Impl*>(this); }
295
296
    Impl &asImp_() { return *static_cast<Impl*>(this); }

297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
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
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
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
     */
    // \{

    /*!
     * \brief get the solid heat capacity in an scv
     * \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
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    auto solidHeatCapacity_(const ElemSol& elemSol,
                            const Problem& problem,
                            const Element& element,
                            const Scv& scv,
                            const SolidState& solidState)
    -> typename std::enable_if_t<!decltype(
            isValid(Detail::hasSolidHeatCapacity<Element, Scv, ElemSol, SolidState>())(problem.spatialParams())
                                            )::value, Scalar>
    {
        return SolidSystem::heatCapacity(solidState);
    }

    /*!
     * \brief get the solid density in an scv
     * \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
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    auto solidDensity_(const ElemSol& elemSol,
                       const Problem& problem,
                       const Element& element,
                       const Scv& scv,
                       const SolidState& solidState)
    -> typename std::enable_if_t<!decltype(
            isValid(Detail::hasSolidDensity<Element, Scv, ElemSol, SolidState>())(problem.spatialParams())
                                            )::value, Scalar>
    {
        return SolidSystem::density(solidState);
    }

    /*!
     * \brief get the solid's thermal conductivity in an scv
     * \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
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    auto solidThermalConductivity_(const ElemSol& elemSol,
                                   const Problem& problem,
                                   const Element& element,
                                   const Scv& scv,
                                   const SolidState& solidState)
    -> typename std::enable_if_t<!decltype(
            isValid(Detail::hasSolidThermalConductivity<Element, Scv, ElemSol, SolidState>())(problem.spatialParams())
                                            )::value, Scalar>
    {
        return SolidSystem::thermalConductivity(solidState);
    }

    // \}

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

    /*!
     * \brief get the solid heat capacity in an scv
     * \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
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    auto solidHeatCapacity_(const ElemSol& elemSol,
                            const Problem& problem,
                            const Element& element,
                            const Scv& scv,
                            const SolidState& solidState)
    -> typename std::enable_if_t<decltype(
            isValid(Detail::hasSolidHeatCapacity<Element, Scv, ElemSol, SolidState>())(problem.spatialParams())
                                            )::value, Scalar>
    {
        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);
    }

    /*!
     * \brief get the solid density in an scv
     * \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
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    auto solidDensity_(const ElemSol& elemSol,
                       const Problem& problem,
                       const Element& element,
                       const Scv& scv,
                       const SolidState& solidState)
    -> typename std::enable_if_t<decltype(
            isValid(Detail::hasSolidDensity<Element, Scv, ElemSol, SolidState>())(problem.spatialParams())
                                            )::value, Scalar>
    {
        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);
    }

    /*!
     * \brief get the solid's heat capacity in an scv
     * \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
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    auto solidThermalConductivity_(const ElemSol& elemSol,
                                   const Problem& problem,
                                   const Element& element,
                                   const Scv& scv,
                                   const SolidState& solidState)
    -> typename std::enable_if_t<decltype(
            isValid(Detail::hasSolidThermalConductivity<Element, Scv, ElemSol, SolidState>())(problem.spatialParams())
                                            )::value, Scalar>
    {
        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);
    }

    // \}

467
468
469
470
471
};

} // end namespace Dumux

#endif