volumevariables.hh 39.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- 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
Katharina Heck's avatar
Katharina Heck committed
21
 * \ingroup PorousmediumNonEquilibriumModel
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
 * \brief This class contains the volume variables required for the
 *        modules which require the specific interfacial area between
 *        fluid phases.
 *
 * This files contains all specializations which use 'real'
 * interfacial areas.
 */
#ifndef DUMUX_NONEQUILIBRIUM_VOLUME_VARIABLES__HH
#define DUMUX_NONEQUILIBRIUM_VOLUME_VARIABLES__HH

#include <dumux/common/dimensionlessnumbers.hh>

#include <dumux/porousmediumflow/volumevariables.hh>

#include <dumux/material/fluidstates/nonequilibrium.hh>
#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>

namespace Dumux
{
Katharina Heck's avatar
Katharina Heck committed
41
42
43
44
45
46
47
/*!
 * \file
 * \ingroup PorousmediumNonEquilibriumModel
 * \brief This class contains the volume variables required for the
 *        modules which require the specific interfacial area between
 *        fluid phases.
 */
48
49
50
51
52
53
54
55
56
57
58
59
60
61
// forward declaration
template <class TypeTag, bool enableChemicalNonEquilibrium ,bool enableThermalNonEquilibrium>
class NonEquilibriumVolumeVariablesImplementation;

template <class TypeTag>
using NonEquilibriumVolumeVariables =
NonEquilibriumVolumeVariablesImplementation<TypeTag, GET_PROP_VALUE(TypeTag, EnableChemicalNonEquilibrium), GET_PROP_VALUE(TypeTag, EnableThermalNonEquilibrium)>;

template <class TypeTag>
class NonEquilibriumVolumeVariablesImplementation<TypeTag, false/*enableChemicalNonEquilibrium*/, false/*enableThermalNonEquilibrium*/>
{
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
62
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
63
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
    using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
    using Element = typename GridView::template Codim<0>::Entity;
    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
    using ParameterCache = typename FluidSystem::ParameterCache;
    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
public:
    void updateInterfacialArea(const ElementSolution& elemSol,
                               const FluidState & fluidState,
                               const ParameterCache &paramCache,
                               const Problem &problem,
                               const Element & element,
                               const SubControlVolume& scv)
     {}

    void updateTemperatures(const ElementSolution& elemSol,
                             const Problem &problem,
                             const Element& element,
                             const SubControlVolume& scv,
                             FluidState& fluidState)
     {}


    void updateMoleFraction(FluidState & actualFluidState,
                            ParameterCache & paramCache,
                            const PrimaryVariables & priVars)
     { }

};

/////////////////////////////////////////////////////////////////////////////////////////////////////////////
// // specialization for the case of NO kinetic mass but kinetic energy transfer of a fluid mixture and solid
// /////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <class TypeTag>
class NonEquilibriumVolumeVariablesImplementation<TypeTag, /*enableChemicalNonEquilibrium*/false, /*enableThermalNonEquilibrium*/true>
{
    using BaseType = PorousMediumFlowVolumeVariables<TypeTag>;
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
105
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
106
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
107
108
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
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
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
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
247
248
249
250
251
252
253
254
255
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
    using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
    using Element = typename GridView::template Codim<0>::Entity;
    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
    using ParameterCache = typename FluidSystem::ParameterCache;
    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);

    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
    enum { nusseltFormulation = GET_PROP_VALUE(TypeTag, NusseltFormulation)} ;
    enum { numEnergyEqFluid = GET_PROP_VALUE(TypeTag, NumEnergyEqFluid) };
    enum { numEnergyEqSolid = GET_PROP_VALUE(TypeTag, NumEnergyEqSolid) };
    enum { temperature0Idx = Indices::temperature0Idx };

    static_assert((numEnergyEqFluid < 2),
                   "This model is a specialization for a energy transfer of a fluid  mixture and a solid");

    using DimLessNum = DimensionlessNumbers<Scalar>;
public:
    /*!
     * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
     */
    void updateInterfacialArea(const ElementSolution& elemSol,
                               const FluidState & fluidState,
                               const ParameterCache &paramCache,
                               const Problem &problem,
                               const Element & element,
                               const SubControlVolume& scv)
    {
        factorMassTransfer_   = problem.spatialParams().factorMassTransfer(element, scv, elemSol);
        factorEnergyTransfer_   = problem.spatialParams().factorEnergyTransfer(element, scv, elemSol);
        characteristicLength_   = problem.spatialParams().characteristicLength(element, scv, elemSol);

        // setting the dimensionless numbers.
        // obtaining the respective quantities.
         const unsigned int vIdxGlobal = scv.dofIndex();
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
        {
            const Scalar darcyMagVelocity     = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, vIdxGlobal);
            const Scalar dynamicViscosity     = fluidState.viscosity(phaseIdx);
            const Scalar density              = fluidState.density(phaseIdx);
            const Scalar kinematicViscosity   = dynamicViscosity / density;
            const Scalar heatCapacity         = FluidSystem::heatCapacity(fluidState,
                                                                          paramCache,
                                                                          phaseIdx);
            const Scalar thermalConductivity  = FluidSystem::thermalConductivity(fluidState,
                                                                           paramCache,
                                                                           phaseIdx);
            const Scalar porosity = problem.spatialParams().porosity(element,
                                                                   scv,
                                                                   elemSol);

            reynoldsNumber_[phaseIdx]   = DimLessNum::reynoldsNumber(darcyMagVelocity,
                                                                     characteristicLength_,
                                                                     kinematicViscosity);

            prandtlNumber_[phaseIdx]    = DimLessNum::prandtlNumber(dynamicViscosity,
                                                                    heatCapacity,
                                                                    thermalConductivity);


            nusseltNumber_[phaseIdx]    = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
                                                                          prandtlNumber_[phaseIdx],
                                                                          porosity,
                                                                          nusseltFormulation);
        }
    }

    /*!
     * \brief Update all quantities for a given control volume
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
     */
    void updateTemperatures(const ElementSolution& elemSol,
                            const Problem &problem,
                            const Element& element,
                            const SubControlVolume& scv,
                            FluidState& fluidState)
    {
        if (numEnergyEqFluid >1)
        for(int phaseIdx=0; phaseIdx < numEnergyEqFluid; ++phaseIdx)
        {
            // retrieve temperature from solution vector
            const Scalar T = BaseType::extractDofPriVars(elemSol, scv)[temperature0Idx + phaseIdx];
            fluidState.setTemperature(phaseIdx, T);
        }
        else
        {
            const Scalar T = BaseType::extractDofPriVars(elemSol, scv)[temperature0Idx];
            fluidState.setTemperature(T);
        }
        for(int solidPhaseIdx = numEnergyEqFluid; solidPhaseIdx < numEnergyEqFluid+numEnergyEqSolid; ++solidPhaseIdx)
        {
            temperatureSolid_ = BaseType::extractDofPriVars(elemSol, scv)[temperature0Idx + solidPhaseIdx];

        }
    }

    void updateMoleFraction(FluidState & actualFluidState,
                            ParameterCache & paramCache,
                            const PrimaryVariables & priVars)
    { }


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

    //! access function Reynolds Number
    const Scalar reynoldsNumber(const unsigned int phaseIdx) const
    { return reynoldsNumber_[phaseIdx]; }

    //! access function Prandtl Number
    const Scalar prandtlNumber(const unsigned int phaseIdx) const
    { return prandtlNumber_[phaseIdx]; }

    //! access function Nusselt Number
    const Scalar nusseltNumber(const unsigned int phaseIdx) const
    { return nusseltNumber_[phaseIdx]; }

    //! access function characteristic length
    const Scalar characteristicLength() const
    { return characteristicLength_; }

    //! access function pre factor energy transfer
    const Scalar factorEnergyTransfer() const
    { return factorEnergyTransfer_; }

    //! access function pre factor mass transfer
    const Scalar factorMassTransfer() const
    { return factorMassTransfer_; }

    /*!
     * \brief If running in valgrind this makes sure that all
     *        quantities in the volume variables are defined.
     */
    void checkDefined() const
    {
#if !defined NDEBUG && HAVE_VALGRIND
        Valgrind::CheckDefined(reynoldsNumber_);
        Valgrind::CheckDefined(prandtlNumber_);
        Valgrind::CheckDefined(nusseltNumber_);
        Valgrind::CheckDefined(characteristicLength_);
        Valgrind::CheckDefined(factorEnergyTransfer_);
        Valgrind::CheckDefined(factorMassTransfer_);
        Valgrind::CheckDefined(temperatureSolid_);
#endif
    }

private:
    //! dimensionless numbers
    Scalar reynoldsNumber_[numPhases];
    Scalar prandtlNumber_[numPhases];
    Scalar nusseltNumber_[numPhases];
    Scalar characteristicLength_;
    Scalar factorEnergyTransfer_;
    Scalar factorMassTransfer_;
    Scalar solidSurface_ ;
    Scalar temperatureSolid_;
};

////////////////////////////////////////////////////////////////////////////////////////////////////
// // specialization for the case of (only) kinetic mass transfer. Be careful, we do not test this!
// ////////////////////////////////////////////////////////////////////////////////////////////////////
template <class TypeTag>
class NonEquilibriumVolumeVariablesImplementation<TypeTag, true/*enableChemicalNonEquilibrium*/, false/*enableThermalNonEquilibrium*/>
{
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
283
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
284
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
    using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
    using Element = typename GridView::template Codim<0>::Entity;
    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
    using ParameterCache = typename FluidSystem::ParameterCache;
    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);

    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
    enum { wPhaseIdx = FluidSystem::wPhaseIdx };
    enum { nPhaseIdx = FluidSystem::nPhaseIdx };
    enum { nCompIdx = FluidSystem::nCompIdx } ;
    enum { wCompIdx = FluidSystem::wCompIdx } ;
    enum { dim = GridView::dimension};
    enum { sherwoodFormulation = GET_PROP_VALUE(TypeTag, SherwoodFormulation)} ;
    enum { moleFrac00Idx = Indices::conti0EqIdx };

    using AwnSurface = typename GET_PROP_TYPE(TypeTag, AwnSurface);
    using AwnSurfaceParams = typename  AwnSurface::Params;

    using DimLessNum = DimensionlessNumbers<Scalar>;
305
    using ConstraintSolver = MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
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
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
public:
    /*!
     * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
     */
    void updateInterfacialArea(const ElementSolution& elemSol,
                               const FluidState & fluidState,
                               const ParameterCache &paramCache,
                               const Problem &problem,
                               const Element & element,
                               const SubControlVolume& scv)
    {
        //obtain parameters for awnsurface
        const AwnSurfaceParams & awnSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element, scv, elemSol) ;

        // obtain (standard) material parameters (needed for the residual saturations)
        const auto &materialParams  = problem.spatialParams().materialLawParams(element, scv, elemSol) ;

        Valgrind::CheckDefined(awnSurfaceParams);
        const Scalar Sw = fluidState.saturation(wPhaseIdx) ;
        const Scalar pc = fluidState.pressure(nPhaseIdx) - fluidState.pressure(wPhaseIdx);

        // so far there is only a model for kinetic mass transfer between fluid phases
        interfacialArea_ = AwnSurface::interfacialArea(awnSurfaceParams, materialParams, Sw, pc );

        Valgrind::CheckDefined(interfacialArea_);

        factorMassTransfer_   = problem.spatialParams().factorMassTransfer(element, scv, elemSol);

        characteristicLength_   = problem.spatialParams().characteristicLength(element, scv, elemSol);
        // setting the dimensionless numbers.
        // obtaining the respective quantities.
        int globalVertexIdx = problem.vertexMapper().subIndex(element, scv, dim);
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
        {
            const Scalar darcyMagVelocity     = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, globalVertexIdx);
            const Scalar dynamicViscosity     = fluidState.viscosity(phaseIdx);
            const Scalar density              = fluidState.density(phaseIdx);
            const Scalar kinematicViscosity   = dynamicViscosity / density;

            // diffusion coefficient of non-wetting component in wetting phase
            const Scalar diffCoeff = FluidSystem::binaryDiffusionCoefficient(fluidState, paramCache, phaseIdx, wCompIdx, nCompIdx);

            reynoldsNumber_[phaseIdx]   = DimLessNum::reynoldsNumber(darcyMagVelocity,
                                                                     characteristicLength_,
                                                                     kinematicViscosity);

            schmidtNumber_[phaseIdx]    = DimLessNum::schmidtNumber(dynamicViscosity,
                                                                    density,
                                                                    diffCoeff);
            sherwoodNumber_[phaseIdx]   = DimLessNum::sherwoodNumber(reynoldsNumber_[phaseIdx],
                                                                     schmidtNumber_[phaseIdx],
                                                                     sherwoodFormulation);
        }
    }

     /*!
     * \brief Update composition of all phases in the mutable
     *        parameters from the primary variables.
     *
     *        \param actualFluidState Container for all the secondary variables concerning the fluids
     *        \param paramCache Container for cache parameters
     *        \param priVars The primary Variables
     *        \param *hint the volume variables, usable for initial guess of composition
     */
    void updateMoleFraction(FluidState & actualFluidState,
                            ParameterCache & paramCache,
                            const PrimaryVariables & priVars)
    {
        // setting the mole fractions of the fluid state
        for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx)
        {
                // set the component mole fractions
                for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
                    actualFluidState.setMoleFraction(phaseIdx,
                           compIdx,
                           priVars[moleFrac00Idx +
                                   phaseIdx*numComponents +
                                   compIdx]);
                }
        }

//          // For using the ... other way of calculating equilibrium
//          THIS IS ONLY FOR silencing Valgrind but is not used in this model
            for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx)
                for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
                    const Scalar phi = FluidSystem::fugacityCoefficient(actualFluidState,
                                                                        paramCache,
                                                                        phaseIdx,
                                                                        compIdx);
                    actualFluidState.setFugacityCoefficient(phaseIdx,
                                                      compIdx,
                                                      phi);
            }

            FluidState equilFluidState; // the fluidState *on the interface* i.e. chemical equilibrium
            equilFluidState.assign(actualFluidState) ;
            ConstraintSolver::solve(equilFluidState,
                                    paramCache,
                                    /*setViscosity=*/false,
                                    /*setEnthalpy=*/false) ;

            // Setting the equilibrium composition (in a kinetic model not necessarily the same as the actual mole fraction)
            for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
                for (int compIdx=0; compIdx< numComponents; ++ compIdx){
                    xEquil_[phaseIdx][compIdx] = equilFluidState.moleFraction(phaseIdx, compIdx);
                }
            }

            // compute densities of all phases
            for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
                const Scalar rho = FluidSystem::density(actualFluidState, paramCache, phaseIdx);
                actualFluidState.setDensity(phaseIdx, rho);
            }

        }

    /*!
     * \brief The mole fraction we would have in the case of chemical equilibrium /
     *        on the interface.
     *
     *     \param phaseIdx The index of the fluid phase
     *     \param compIdx The local index of the component
     */
    const Scalar xEquil(const unsigned int phaseIdx, const unsigned int compIdx) const
    {
        return xEquil_[phaseIdx][compIdx] ;
    }

    /*!
     * \brief The specific interfacial area between two fluid phases [m^2 / m^3]
     */
    const Scalar interfacialArea(const unsigned int phaseIIdx, const unsigned int phaseJIdx) const
    {
        // so far there is only a model for kinetic mass transfer between fluid phases
        assert((phaseIIdx == nPhaseIdx && phaseJIdx == wPhaseIdx)
              || (phaseIIdx == wPhaseIdx && phaseJIdx == nPhaseIdx));
        return interfacialArea_;
    }

    //! access function Reynolds Number
    const Scalar reynoldsNumber(const unsigned int phaseIdx) const
    { return reynoldsNumber_[phaseIdx]; }

    //! access function Schmidt Number
    const Scalar schmidtNumber(const unsigned int phaseIdx) const
    { return schmidtNumber_[phaseIdx]; }

    //! access function Sherwood Number
    const Scalar sherwoodNumber(const unsigned int phaseIdx) const
    { return sherwoodNumber_[phaseIdx]; }

    //! access function characteristic length
    const Scalar characteristicLength() const
    { return characteristicLength_; }

    //! access function pre factor mass transfer
    const Scalar factorMassTransfer() const
    { return factorMassTransfer_; }

    /*!
     * \brief If running in valgrind this makes sure that all
     *        quantities in the volume variables are defined.
     */
    void checkDefined() const
    {
        Valgrind::CheckDefined(interfacialArea_);
        Valgrind::CheckDefined(characteristicLength_);
        Valgrind::CheckDefined(factorMassTransfer_);
        Valgrind::CheckDefined(reynoldsNumber_);
        Valgrind::CheckDefined(schmidtNumber_);
        Valgrind::CheckDefined(sherwoodNumber_);
    }

private:
    Scalar characteristicLength_;
    Scalar factorMassTransfer_;
    Scalar solidSurface_ ;
    Scalar interfacialArea_ ;
    Scalar sherwoodNumber_[numPhases] ;
    Scalar schmidtNumber_[numPhases] ;
    Scalar reynoldsNumber_[numPhases] ;
    Scalar xEquil_[numPhases][numComponents];
};

//this is a specialization where everything is in non-equilibrium. we have to do our own stuff for the interfacial area but can use the rest from the others
template <class TypeTag>
class NonEquilibriumVolumeVariablesImplementation<TypeTag, true/*enableChemicalNonEquilibrium*/, true/*enableThermalNonEquilibrium*/>
{
    using BaseType = PorousMediumFlowVolumeVariables<TypeTag>;
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
500
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
501
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
    using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
    using Element = typename GridView::template Codim<0>::Entity;
    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
    using ParameterCache = typename FluidSystem::ParameterCache;
    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);

    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
    enum { wPhaseIdx = FluidSystem::wPhaseIdx };
    enum { nPhaseIdx = FluidSystem::nPhaseIdx };
    enum { sPhaseIdx = FluidSystem::sPhaseIdx };
    enum { nCompIdx = FluidSystem::nCompIdx } ;
    enum { wCompIdx = FluidSystem::wCompIdx } ;
    enum { nusseltFormulation = GET_PROP_VALUE(TypeTag, NusseltFormulation)} ;
    enum { sherwoodFormulation = GET_PROP_VALUE(TypeTag, SherwoodFormulation)} ;
    enum { numEnergyEqFluid = GET_PROP_VALUE(TypeTag, NumEnergyEqFluid) };
    enum { numEnergyEqSolid = GET_PROP_VALUE(TypeTag, NumEnergyEqSolid) };
    enum { temperature0Idx = Indices::temperature0Idx };
    enum { moleFrac00Idx = Indices::conti0EqIdx };

    static_assert((numEnergyEqFluid > 1),
                   "This model only deals with energy transfer between two fluids and one solid phase");
    using DimLessNum = DimensionlessNumbers<Scalar>;

    using AwnSurface = typename GET_PROP_TYPE(TypeTag, AwnSurface);
    using AwnSurfaceParams = typename  AwnSurface::Params;
    using AwsSurface = typename GET_PROP_TYPE(TypeTag, AwsSurface);
    using AwsSurfaceParams = typename  AwsSurface::Params;
    using AnsSurface = typename GET_PROP_TYPE(TypeTag, AnsSurface);
    using AnsSurfaceParams = typename  AnsSurface::Params;

533
    using ConstraintSolver = MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843


public:
    /*!
     * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
     */
    void updateInterfacialArea(const ElementSolution& elemSol,
                                const FluidState & fluidState,
                                const ParameterCache &paramCache,
                                const Problem &problem,
                                const Element & element,
                                const SubControlVolume& scv)
    {
        // obtain (standard) material parameters (needed for the residual saturations)
       const auto& materialParams =
            problem.spatialParams().materialLawParams(element, scv, elemSol);

        //obtain parameters for interfacial area constitutive relations
        const AwnSurfaceParams & aWettingNonWettingSurfaceParams
               = problem.spatialParams().aWettingNonWettingSurfaceParams(element, scv, elemSol);
        const AnsSurfaceParams & aNonWettingSolidSurfaceParams
               = problem.spatialParams().aNonWettingSolidSurfaceParams(element, scv, elemSol);

        Valgrind::CheckDefined(aWettingNonWettingSurfaceParams);
        Valgrind::CheckDefined(aNonWettingSolidSurfaceParams);

        const Scalar pc = fluidState.pressure(nPhaseIdx) - fluidState.pressure(wPhaseIdx);
        const Scalar Sw = fluidState.saturation(wPhaseIdx);
        Valgrind::CheckDefined(Sw);
        Valgrind::CheckDefined(pc);
        Scalar awn;

#define AwnRegul 0
        // This regularizes the interfacial area between the fluid phases.
        // This makes sure, that
        // a) some saturation cannot be lost: Never leave two phase region.
        // b) We cannot leave the fit region: no crazy (e.g. negative) values possible
//        const Scalar Swr =  aWettingNonWettingSurfaceParams.Swr() ;
//        const Scalar Snr =  aWettingNonWettingSurfaceParams.Snr() ;

        // this just leads to a stalling newton error as soon as this kicks in.
        // May be a spline or sth like this would help, but I do not which derivatives
        // to specify.
#if AwnRegul
        if(Sw < 5e-3 ) // or Sw > (1.-1e-5 )
        {
            awn = 0. ; // 10.; //
        }
        else
#endif
        awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc ); // 10.; //

        interfacialArea_[wPhaseIdx][nPhaseIdx] = awn ; //10. ;//
        interfacialArea_[nPhaseIdx][wPhaseIdx] = interfacialArea_[wPhaseIdx][nPhaseIdx];
        interfacialArea_[wPhaseIdx][wPhaseIdx] = 0. ;

        Scalar ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams,Sw, pc ); // 10.; //
//        if (ans <0 )
//            ans = 0 ;

// Switch for using a a_{wn} relations that has some "maximum capillary pressure" as parameter.
// That value is obtained by regularization of the pc(Sw) function.
#if USE_PCMAX
       const Scalar pcMax = problem.spatialParams().pcMax(element,
                                                          scv,
                                                          elemSol);
        // I know the solid surface from the pore network. But it is more consistent to use the fit value.
        solidSurface_   = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax );
        Valgrind::CheckDefined(solidSurface_);
#endif
        interfacialArea_[nPhaseIdx][sPhaseIdx] = ans ; //10. ; //
        interfacialArea_[sPhaseIdx][nPhaseIdx] = interfacialArea_[nPhaseIdx][sPhaseIdx];
        interfacialArea_[nPhaseIdx][nPhaseIdx] = 0. ;

#if USE_PCMAX
        const Scalar aws = solidSurface_ - ans ;
        interfacialArea_[wPhaseIdx][sPhaseIdx] = aws ; //10. ; //
        interfacialArea_[sPhaseIdx][wPhaseIdx] = interfacialArea_[wPhaseIdx][sPhaseIdx];
        interfacialArea_[sPhaseIdx][sPhaseIdx] = 0. ;
#else
        const AwsSurfaceParams & aWettingSolidSurfaceParams
               = problem.spatialParams().aWettingSolidSurfaceParams();
        Valgrind::CheckDefined(aWettingSolidSurfaceParams);
        const Scalar aws = AwsSurface::interfacialArea(aWettingSolidSurfaceParams,materialParams, Sw, pc ); // 10.; //
        interfacialArea_[wPhaseIdx][sPhaseIdx] = aws ;
        interfacialArea_[sPhaseIdx][wPhaseIdx] = interfacialArea_[wPhaseIdx][sPhaseIdx];
        interfacialArea_[sPhaseIdx][sPhaseIdx] = 0. ;
#endif

        Valgrind::CheckDefined(interfacialArea_);

        factorMassTransfer_   = problem.spatialParams().factorMassTransfer(element, scv, elemSol);

        factorEnergyTransfer_   = problem.spatialParams().factorEnergyTransfer(element, scv, elemSol);

        characteristicLength_   = problem.spatialParams().characteristicLength(element, scv, elemSol);


        const unsigned int vIdxGlobal = scv.dofIndex();
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            const Scalar darcyMagVelocity     = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, vIdxGlobal);
            const Scalar dynamicViscosity     = fluidState.viscosity(phaseIdx);
            const Scalar density              = fluidState.density(phaseIdx);
            const Scalar kinematicViscosity   = dynamicViscosity / density;
            const Scalar heatCapacity         = FluidSystem::heatCapacity(fluidState, paramCache, phaseIdx);
            const Scalar thermalConductivity  = FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);


            // diffusion coefficient of non-wetting component in wetting phase
            const Scalar diffCoeff =  FluidSystem::binaryDiffusionCoefficient(fluidState, paramCache, phaseIdx, wCompIdx, nCompIdx);
            const Scalar porosity = problem.spatialParams().porosity(element,
                                                                    scv,
                                                                   elemSol);

            reynoldsNumber_[phaseIdx]   = DimLessNum::reynoldsNumber(darcyMagVelocity,
                                                                     characteristicLength_,
                                                                     kinematicViscosity);

            prandtlNumber_[phaseIdx]    = DimLessNum::prandtlNumber(dynamicViscosity,
                                                                    heatCapacity,
                                                                    thermalConductivity);

            nusseltNumber_[phaseIdx]    = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
                                                                          prandtlNumber_[phaseIdx],
                                                                          porosity,
                                                                          nusseltFormulation);

            schmidtNumber_[phaseIdx]    = DimLessNum::schmidtNumber(dynamicViscosity,
                                                                    density,
                                                                    diffCoeff);

            // If Diffusion is not enabled, Sherwood is divided by zero
            sherwoodNumber_[phaseIdx]   = DimLessNum::sherwoodNumber(reynoldsNumber_[phaseIdx],
                                                                      schmidtNumber_[phaseIdx],
                                                                      sherwoodFormulation);
        }
    }

    void updateMoleFraction(FluidState & actualFluidState,
                            ParameterCache & paramCache,
                            const PrimaryVariables & priVars)
    {
        // setting the mole fractions of the fluid state
        for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
                // set the component mole fractions
                for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
                    actualFluidState.setMoleFraction(phaseIdx,
                           compIdx,
                           priVars[moleFrac00Idx +
                                   phaseIdx*numComponents +
                                   compIdx]);
                }
            }

//            // For using the ... other way of calculating equilibrium
//             THIS IS ONLY FOR silencing Valgrind but is not used in this model
            for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx)
                for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
                    const Scalar phi = FluidSystem::fugacityCoefficient(actualFluidState,
                                                                        paramCache,
                                                                        phaseIdx,
                                                                        compIdx);
                    actualFluidState.setFugacityCoefficient(phaseIdx,
                                                      compIdx,
                                                      phi);
            }

            FluidState equilFluidState; // the fluidState *on the interface* i.e. chemical equilibrium
            equilFluidState.assign(actualFluidState) ;
            ConstraintSolver::solve(equilFluidState,
                                    paramCache,
                                    /*setViscosity=*/false,
                                    /*setEnthalpy=*/false) ;

            // Setting the equilibrium composition (in a kinetic model not necessarily the same as the actual mole fraction)
            for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
                for (int compIdx=0; compIdx< numComponents; ++ compIdx){
                    xEquil_[phaseIdx][compIdx] = equilFluidState.moleFraction(phaseIdx, compIdx);
                }
            }

            // compute densities of all phases
            for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
                const Scalar rho = FluidSystem::density(actualFluidState, paramCache, phaseIdx);
                actualFluidState.setDensity(phaseIdx, rho);
            }

     }

    void updateTemperatures(const ElementSolution& elemSol,
                            const Problem &problem,
                            const Element& element,
                            const SubControlVolume& scv,
                            FluidState& fluidState)
    {
        for(int phaseIdx=0; phaseIdx < numEnergyEqFluid; ++phaseIdx)
        {
            // retrieve temperature from solution vector
            const Scalar T = BaseType::extractDofPriVars(elemSol, scv)[temperature0Idx + phaseIdx];
            fluidState.setTemperature(phaseIdx, T);
        }
        for(int solidPhaseIdx = numEnergyEqFluid; solidPhaseIdx < numEnergyEqFluid+numEnergyEqSolid; ++solidPhaseIdx)
        {
            temperatureSolid_ = BaseType::extractDofPriVars(elemSol, scv)[temperature0Idx + solidPhaseIdx];
        }
    }

    /*!
     * \brief Returns the temperature in fluid / solid phase(s)
     *        the sub-control volume.
     * \param phaseIdx The local index of the phases
     */
    const Scalar temperatureSolid() const
    { return temperatureSolid_; }
    /*!
     * \brief The mole fraction we would have in the case of chemical equilibrium /
     *        on the interface.
     *
     *     \param phaseIdx The index of the fluid phase
     *     \param compIdx The local index of the component
     */
    const Scalar xEquil(const unsigned int phaseIdx, const unsigned int compIdx) const
    {
        return xEquil_[phaseIdx][compIdx] ;
    }

    /*!
     * \brief The specific interfacial area between two fluid phases [m^2 / m^3]
     *
     * This is _only_ required by the kinetic mass/energy modules
     *
     */
    const Scalar interfacialArea(const unsigned int phaseIIdx, const unsigned int phaseJIdx) const
    {
        // there is no interfacial area between a phase and itself
        assert(phaseIIdx not_eq phaseJIdx);
        return interfacialArea_[phaseIIdx][phaseJIdx];
    }

    //! access function Reynolds Number
    const Scalar reynoldsNumber(const unsigned int phaseIdx) const
    { return reynoldsNumber_[phaseIdx]; }

    //! access function Prandtl Number
    const Scalar prandtlNumber(const unsigned int phaseIdx) const
    { return prandtlNumber_[phaseIdx]; }

    //! access function Nusselt Number
    const Scalar nusseltNumber(const unsigned int phaseIdx) const
    { return nusseltNumber_[phaseIdx]; }

    //! access function Schmidt Number
    const Scalar schmidtNumber(const unsigned int phaseIdx) const
    { return schmidtNumber_[phaseIdx]; }

    //! access function Sherwood Number
    const Scalar sherwoodNumber(const unsigned int phaseIdx) const
    { return sherwoodNumber_[phaseIdx]; }

    //! access function characteristic length
    const Scalar characteristicLength() const
    { return characteristicLength_; }

    //! access function pre factor energy transfer
    const Scalar factorEnergyTransfer() const
    { return factorEnergyTransfer_; }

    //! access function pre factor mass transfer
    const Scalar factorMassTransfer() const
    { return factorMassTransfer_; }

    /*!
     * \brief If running in valgrind this makes sure that all
     *        quantities in the volume variables are defined.
     */
    void checkDefined() const
    {
#if !defined NDEBUG && HAVE_VALGRIND
        Valgrind::CheckDefined(reynoldsNumber_);
        Valgrind::CheckDefined(prandtlNumber_);
        Valgrind::CheckDefined(nusseltNumber_);
        Valgrind::CheckDefined(schmidtNumber_);
        Valgrind::CheckDefined(sherwoodNumber_);
        Valgrind::CheckDefined(characteristicLength_);
        Valgrind::CheckDefined(factorEnergyTransfer_);
        Valgrind::CheckDefined(factorMassTransfer_);
        Valgrind::CheckDefined(interfacialArea_);
#endif
    }

private:
    //! dimensionless numbers
    Scalar reynoldsNumber_[numPhases];
    Scalar prandtlNumber_[numPhases];
    Scalar nusseltNumber_[numPhases];
    Scalar schmidtNumber_[numPhases];
    Scalar sherwoodNumber_[numPhases];
    Scalar characteristicLength_;
    Scalar factorEnergyTransfer_;
    Scalar factorMassTransfer_;
    Scalar solidSurface_ ;
    Scalar interfacialArea_[numPhases+1][numPhases+1];
    Scalar xEquil_[numPhases][numComponents];
    Scalar temperatureSolid_;

};
//
} // namespace Dumux

#endif