volumevariables.hh 40.6 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 NonEquilibriumModel
22
23
24
25
26
27
28
 * \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.
 */
29

30
31
#ifndef DUMUX_NONEQUILIBRIUM_VOLUME_VARIABLES_HH
#define DUMUX_NONEQUILIBRIUM_VOLUME_VARIABLES_HH
32

Bernd Flemisch's avatar
Bernd Flemisch committed
33
#include <cassert>
34
#include <array>
Bernd Flemisch's avatar
Bernd Flemisch committed
35

36
#include <dumux/common/dimensionlessnumbers.hh>
37
#include <dumux/common/parameters.hh>
38

39
40
namespace Dumux {

Katharina Heck's avatar
Katharina Heck committed
41
/*!
42
 * \ingroup NonEquilibriumModel
Katharina Heck's avatar
Katharina Heck committed
43
44
45
46
 * \brief This class contains the volume variables required for the
 *        modules which require the specific interfacial area between
 *        fluid phases.
 */
47
48
template<class Traits, class EquilibriumVolumeVariables, bool enableChemicalNonEquilibrium,
         bool enableThermalNonEquilibrium, int numEnergyEqFluid>
49
50
class NonEquilibriumVolumeVariablesImplementation;

51
template<class Traits, class EquilibriumVolumeVariables>
52
using NonEquilibriumVolumeVariables =
53
54
55
56
57
      NonEquilibriumVolumeVariablesImplementation<Traits,
                                                  EquilibriumVolumeVariables,
                                                  Traits::ModelTraits::enableChemicalNonEquilibrium(),
                                                  Traits::ModelTraits::enableThermalNonEquilibrium(),
                                                  Traits::ModelTraits::numEnergyEqFluid()>;
58

59
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
60
// specialization for the case of NO kinetic mass but kinetic energy transfer of  two fluids and a solid
61
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
62
template<class Traits, class EquilibriumVolumeVariables>
63
64
65
66
67
class NonEquilibriumVolumeVariablesImplementation<Traits,
                                                  EquilibriumVolumeVariables,
                                                  false/*chemicalNonEquilibrium?*/,
                                                  true/*thermalNonEquilibrium?*/,
                                                  2>
68
: public EquilibriumVolumeVariables
69
{
70
    using ParentType = EquilibriumVolumeVariables;
71
72
73
74
    using ParameterCache = typename Traits::FluidSystem::ParameterCache;
    using Scalar = typename Traits::PrimaryVariables::value_type;

    using ModelTraits = typename Traits::ModelTraits;
75

76
    using FS = typename Traits::FluidSystem;
77
78
79
    static constexpr auto numEnergyEqFluid = ModelTraits::numEnergyEqFluid();
    static constexpr auto numEnergyEqSolid = ModelTraits::numEnergyEqSolid();

80
81
82
    static constexpr auto phase0Idx = FS::phase0Idx;
    static constexpr auto phase1Idx = FS::phase1Idx;
    static constexpr auto sPhaseIdx = FS::numPhases;
83
84

    using DimLessNum = DimensionlessNumbers<Scalar>;
85

86
87
88
89
    using NumFluidPhasesArray = std::array<Scalar, ModelTraits::numFluidPhases()>;
    using InterfacialAreasArray =  std::array<std::array<Scalar, ModelTraits::numFluidPhases()+numEnergyEqSolid>,
                                              ModelTraits::numFluidPhases()+numEnergyEqSolid>;

90
public:
91
     using FluidState = typename Traits::FluidState;
92
     using Indices = typename ModelTraits::Indices;
93
94
95
96
    /*!
     * \brief Update all quantities for a given control volume
     *
     * \param elemSol A vector containing all primary variables connected to the element
97
     * \param problem The object specifying the problem which ought to be simulated
98
99
100
101
102
103
104
105
106
107
108
109
110
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    void update(const ElemSol &elemSol,
                const Problem &problem,
                const Element &element,
                const Scv& scv)
    {
        // Update parent type (also completes the fluid state)
        ParentType::update(elemSol, problem, element, scv);

        ParameterCache paramCache;
111
        paramCache.updateAll(this->fluidState());
112
        updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
113
        updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
114
115
    }

116
    /*!
117
     * \brief Updates dimensionless numbers
118
119
120
121
122
123
124
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param fluidState Container for all the secondary variables concerning the fluids
     * \param paramCache The parameter cache corresponding to the fluid state
     * \param problem The problem to be solved
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
125
     */
126
    template<class ElemSol, class Problem, class Element, class Scv>
127
    void updateDimLessNumbers(const ElemSol& elemSol,
128
129
130
131
132
                              const FluidState& fluidState,
                              const ParameterCache& paramCache,
                              const Problem& problem,
                              const Element& element,
                              const Scv& scv)
133
    {
134
135
        factorEnergyTransfer_ = problem.spatialParams().factorEnergyTransfer(element, scv, elemSol);
        characteristicLength_ = problem.spatialParams().characteristicLength(element, scv, elemSol);
136

137
138
        // set the dimensionless numbers and obtain respective quantities
        const unsigned int vIdxGlobal = scv.dofIndex();
139
        for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
140
        {
141
142
143
144
            const auto darcyMagVelocity = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, vIdxGlobal);
            const auto dynamicViscosity = fluidState.viscosity(phaseIdx);
            const auto density = fluidState.density(phaseIdx);
            const auto kinematicViscosity = dynamicViscosity/density;
145
146

            using FluidSystem = typename Traits::FluidSystem;
147
            const auto heatCapacity = FluidSystem::heatCapacity(fluidState, paramCache, phaseIdx);
148
            const auto thermalConductivity = FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);
149
            const auto porosity = this->porosity();
150

151
            reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_,kinematicViscosity);
152
153
154
155
156
            prandtlNumber_[phaseIdx] = DimLessNum::prandtlNumber(dynamicViscosity, heatCapacity, thermalConductivity);
            nusseltNumber_[phaseIdx] = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
                                                                       prandtlNumber_[phaseIdx],
                                                                       porosity,
                                                                       ModelTraits::nusseltFormulation());
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
     * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param fluidState Container for all the secondary variables concerning the fluids
     * \param paramCache The parameter cache corresponding to the fluid state
     * \param problem The problem to be solved
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    void updateInterfacialArea(const ElemSol& elemSol,
                               const FluidState& fluidState,
                               const ParameterCache& paramCache,
                               const Problem& problem,
                               const Element& element,
                               const Scv& 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
183
        const auto& aWettingNonWettingSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element, scv, elemSol);
184
185
186
187
188

        const Scalar pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
        const Scalar Sw = fluidState.saturation(phase0Idx);

        using AwnSurface = typename Problem::SpatialParams::AwnSurface;
189
        const auto awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc);
190
191
192
193
194
        interfacialArea_[phase0Idx][phase1Idx] = awn;
        interfacialArea_[phase1Idx][phase0Idx] = interfacialArea_[phase0Idx][phase1Idx];
        interfacialArea_[phase0Idx][phase0Idx] = 0.;

        using AnsSurface = typename Problem::SpatialParams::AnsSurface;
195
196
        const auto& aNonWettingSolidSurfaceParams = problem.spatialParams().aNonWettingSolidSurfaceParams(element, scv, elemSol);
        const auto ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, Sw, pc);
197

198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
        // 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.
        static const bool computeAwsFromAnsAndPcMax = getParam<bool>("SpatialParams.ComputeAwsFromAnsAndPcMax", true);
        if (computeAwsFromAnsAndPcMax)
        {
            // I know the solid surface from the pore network. But it is more consistent to use the fit value.
            const Scalar pcMax = aWettingNonWettingSurfaceParams.pcMax();
            const auto solidSurface = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax);
            interfacialArea_[phase0Idx][sPhaseIdx] = solidSurface - ans;
        }
        else
        {
            using AwsSurface = typename Problem::SpatialParams::AwsSurface;
            const auto& aWettingSolidSurfaceParams = problem.spatialParams().aWettingSolidSurfaceParams(element, scv, elemSol);
            interfacialArea_[phase0Idx][sPhaseIdx] = AwsSurface::interfacialArea(aWettingSolidSurfaceParams, materialParams, Sw, pc);
        }
214
215
216

        interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
        interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
217

218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
        interfacialArea_[phase1Idx][sPhaseIdx] = ans;
        interfacialArea_[sPhaseIdx][phase1Idx] = interfacialArea_[phase1Idx][sPhaseIdx];
        interfacialArea_[phase1Idx][phase1Idx] = 0.;
    }

    /*!
     * \brief The specific interfacial area between two fluid phases [m^2 / m^3]
     * \note 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];
    }

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

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

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

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

250
    //! access function pre factor energy transfer
251
252
    const Scalar factorEnergyTransfer() const
    { return factorEnergyTransfer_; }
253
254
255

private:
    //! dimensionless numbers
256
257
258
    NumFluidPhasesArray reynoldsNumber_;
    NumFluidPhasesArray prandtlNumber_;
    NumFluidPhasesArray nusseltNumber_;
259

260
261
    Scalar characteristicLength_;
    Scalar factorEnergyTransfer_;
262
    InterfacialAreasArray interfacialArea_;
263
264
};

265
266
267
268
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
// specialization for the case of NO kinetic mass but kinetic energy transfer of  a fluid mixture and a solid
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
template<class Traits, class EquilibriumVolumeVariables>
269
270
271
272
273
274
class NonEquilibriumVolumeVariablesImplementation<Traits,
                                                  EquilibriumVolumeVariables,
                                                  false/*chemicalNonEquilibrium?*/,
                                                  true/*thermalNonEquilibrium?*/,
                                                  1>
: public EquilibriumVolumeVariables
275
276
277
278
279
280
281
282
283
284
285
286
{
    using ParentType = EquilibriumVolumeVariables;
    using ParameterCache = typename Traits::FluidSystem::ParameterCache;
    using Scalar = typename Traits::PrimaryVariables::value_type;

    using ModelTraits = typename Traits::ModelTraits;
    using FS = typename Traits::FluidSystem;
    static constexpr auto numEnergyEqFluid = ModelTraits::numEnergyEqFluid();
    static constexpr auto numEnergyEqSolid = ModelTraits::numEnergyEqSolid();

    using DimLessNum = DimensionlessNumbers<Scalar>;

287
288
    using NumFluidPhasesArray = std::array<Scalar, ModelTraits::numFluidPhases()>;

289
public:
290
291
292
    using Indices = typename ModelTraits::Indices;
    using FluidState = typename Traits::FluidState;

293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
    /*!
     * \brief Update all quantities for a given control volume
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param problem The object specifying the problem which ought to
     *                be simulated
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    void update(const ElemSol &elemSol,
                const Problem &problem,
                const Element &element,
                const Scv& scv)
    {
        // Update parent type (also completes the fluid state)
        ParentType::update(elemSol, problem, element, scv);

        ParameterCache paramCache;
        paramCache.updateAll(this->fluidState());
313
314
        //only update of DimLessNumbers is necessary here, as interfacial area
        // is easy due to only one fluid with a solid and is directly computed in localresidual
315
        updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
316
        updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
317
318
319
320
321
322
323
324
325
326
327
328
329
330
    }

    /*!
     * \brief Updates dimensionless numbers
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param fluidState Container for all the secondary variables concerning the fluids
     * \param paramCache The parameter cache corresponding to the fluid state
     * \param problem The problem to be solved
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    void updateDimLessNumbers(const ElemSol& elemSol,
331
332
333
334
335
                              const FluidState& fluidState,
                              const ParameterCache& paramCache,
                              const Problem& problem,
                              const Element& element,
                              const Scv& scv)
336
337
338
339
340
341
    {
        factorEnergyTransfer_ = problem.spatialParams().factorEnergyTransfer(element, scv, elemSol);
        characteristicLength_ = problem.spatialParams().characteristicLength(element, scv, elemSol);

        // set the dimensionless numbers and obtain respective quantities
        const unsigned int vIdxGlobal = scv.dofIndex();
342
        for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
343
        {
344
345
346
347
            const auto darcyMagVelocity = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, vIdxGlobal);
            const auto dynamicViscosity = fluidState.viscosity(phaseIdx);
            const auto density = fluidState.density(phaseIdx);
            const auto kinematicViscosity = dynamicViscosity/density;
348
349

            using FluidSystem = typename Traits::FluidSystem;
350
            const auto heatCapacity = FluidSystem::heatCapacity(fluidState, paramCache, phaseIdx);
351
            const auto thermalConductivity = FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);
352
            const auto porosity = this->porosity();
353
354

            reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_,kinematicViscosity);
355
356
357
358
359
            prandtlNumber_[phaseIdx] = DimLessNum::prandtlNumber(dynamicViscosity, heatCapacity, thermalConductivity);
            nusseltNumber_[phaseIdx] = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
                                                                       prandtlNumber_[phaseIdx],
                                                                       porosity,
                                                                       ModelTraits::nusseltFormulation());
360
361
362
        }
    }

363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382

    /*!
     * \brief Updates the volume specific interfacial area [m^2 / m^3] between the solid and the fluid phase.
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param fluidState Container for all the secondary variables concerning the fluids
     * \param paramCache The parameter cache corresponding to the fluid state
     * \param problem The problem to be solved
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    void updateInterfacialArea(const ElemSol& elemSol,
                               const FluidState& fluidState,
                               const ParameterCache& paramCache,
                               const Problem& problem,
                               const Element& element,
                               const Scv& scv)
    {
       using FluidSolidInterfacialAreaFormulation = typename Problem::SpatialParams::FluidSolidInterfacialAreaFormulation;
383
       interfacialArea_ = FluidSolidInterfacialAreaFormulation::fluidSolidInterfacialArea(this->porosity(), characteristicLength());
384
385
    }

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

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

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

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

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

406
407
    const Scalar fluidSolidInterfacialArea() const
    {return interfacialArea_;}
408

409
410
private:
    //! dimensionless numbers
411
412
413
    NumFluidPhasesArray reynoldsNumber_;
    NumFluidPhasesArray prandtlNumber_;
    NumFluidPhasesArray nusseltNumber_;
414
415
416

    Scalar characteristicLength_;
    Scalar factorEnergyTransfer_;
417
    Scalar interfacialArea_ ;
418
419
};

420
////////////////////////////////////////////////////////////////////////////////////////////////////
421
422
// specialization for the case of (only) kinetic mass transfer. Be careful, we do not test this!
////////////////////////////////////////////////////////////////////////////////////////////////////
423
template<class Traits, class EquilibriumVolumeVariables>
424
425
426
427
428
429
class NonEquilibriumVolumeVariablesImplementation<Traits,
                                                  EquilibriumVolumeVariables,
                                                  true/*chemicalNonEquilibrium?*/,
                                                  false/*thermalNonEquilibrium?*/,
                                                  0>
: public EquilibriumVolumeVariables
430
{
431
    using ParentType = EquilibriumVolumeVariables;
432
    using FluidState = typename Traits::FluidState;
433
    using FS = typename Traits::FluidSystem;
434
435
436
437
438
    using ParameterCache = typename Traits::FluidSystem::ParameterCache;
    using Scalar = typename Traits::PrimaryVariables::value_type;

    using ModelTraits = typename Traits::ModelTraits;

439
440
441
442
    static constexpr auto phase0Idx = FS::phase0Idx;
    static constexpr auto phase1Idx = FS::phase1Idx;
    static constexpr auto wCompIdx = FS::comp0Idx;
    static constexpr auto nCompIdx = FS::comp1Idx;
443
444

    using DimLessNum = DimensionlessNumbers<Scalar>;
445
446
447

    using NumFluidPhasesArray = std::array<Scalar, ModelTraits::numFluidPhases()>;

448
public:
449
    using Indices = typename ModelTraits::Indices;
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
    /*!
     * \brief Update all quantities for a given control volume
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param problem The object specifying the problem which ought to
     *                be simulated
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    void update(const ElemSol &elemSol,
                const Problem &problem,
                const Element &element,
                const Scv& scv)
    {
        // Update parent type (also completes the fluid state)
        ParentType::update(elemSol, problem, element, scv);

        ParameterCache paramCache;
469
        paramCache.updateAll(this->fluidState());
470
        updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
471
        updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
472
473
    }

474
475
    /*!
     * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
476
477
478
479
480
481
482
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param fluidState Container for all the secondary variables concerning the fluids
     * \param paramCache The parameter cache corresponding to the fluid state
     * \param problem The problem to be solved
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
483
     */
484
    template<class ElemSol, class Problem, class Element, class Scv>
485
486
487
488
489
490
    void updateDimLessNumbers(const ElemSol& elemSol,
                              const FluidState& fluidState,
                              const ParameterCache& paramCache,
                              const Problem& problem,
                              const Element& element,
                              const Scv& scv)
491
    {
492
493
        factorMassTransfer_ = problem.spatialParams().factorMassTransfer(element, scv, elemSol);
        characteristicLength_ = problem.spatialParams().characteristicLength(element, scv, elemSol);
494

495
        // set the dimensionless numbers and obtain respective quantities.
496
        const unsigned int vIdxGlobal = scv.dofIndex();
497
        for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
498
        {
499
500
501
502
            const auto darcyMagVelocity = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, vIdxGlobal);
            const auto dynamicViscosity = fluidState.viscosity(phaseIdx);
            const auto density = fluidState.density(phaseIdx);
            const auto kinematicViscosity = dynamicViscosity/density;
503
504

            // diffusion coefficient of non-wetting component in wetting phase
505
            using FluidSystem = typename Traits::FluidSystem;
506
507
508
            const auto diffCoeff = FluidSystem::binaryDiffusionCoefficient(fluidState,
                                                                           paramCache,
                                                                           phaseIdx,
509
510
                                                                           wCompIdx,
                                                                           nCompIdx);
511
512

            reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_, kinematicViscosity);
513
            schmidtNumber_[phaseIdx] = DimLessNum::schmidtNumber(dynamicViscosity, density, diffCoeff);
514
515
516
            sherwoodNumber_[phaseIdx] = DimLessNum::sherwoodNumber(reynoldsNumber_[phaseIdx],
                                                                   schmidtNumber_[phaseIdx],
                                                                   ModelTraits::sherwoodFormulation());
517
518
519
        }
    }

520
    /*!
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
     * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param fluidState Container for all the secondary variables concerning the fluids
     * \param paramCache The parameter cache corresponding to the fluid state
     * \param problem The problem to be solved
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    void updateInterfacialArea(const ElemSol& elemSol,
                               const FluidState& fluidState,
                               const ParameterCache& paramCache,
                               const Problem& problem,
                               const Element& element,
                               const Scv& scv)
    {
        // obtain parameters for awnsurface and material law
        const auto& awnSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element, scv, elemSol) ;
540
        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol) ;
541
542
543
544

        const auto Sw = fluidState.saturation(phase0Idx) ;
        const auto pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);

545
546
        // when we only consider chemical non-equilibrium there is only mass transfer between
        // the fluid phases, so in 2p only interfacial area between wetting and non-wetting
547
548
549
550
        using AwnSurface = typename Problem::SpatialParams::AwnSurface;
        interfacialArea_ = AwnSurface::interfacialArea(awnSurfaceParams, materialParams, Sw, pc);
    }

551
552
553
554
555
556
    /*!
     * \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
557
558
        assert( (phaseIIdx == phase1Idx && phaseJIdx == phase0Idx)
                || (phaseIIdx == phase0Idx && phaseJIdx == phase1Idx) );
559
560
561
562
        return interfacialArea_;
    }

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

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

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

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

578
    //! access function pre factor mass transfer
579
580
    const Scalar factorMassTransfer() const
    { return factorMassTransfer_; }
581
582
583
584
585

private:
    Scalar characteristicLength_;
    Scalar factorMassTransfer_;
    Scalar interfacialArea_ ;
586
587
588
    NumFluidPhasesArray sherwoodNumber_;
    NumFluidPhasesArray schmidtNumber_;
    NumFluidPhasesArray reynoldsNumber_;
589
590
};

591
// Specialization where everything is in non-equilibrium.
592
template<class Traits, class EquilibriumVolumeVariables>
593
594
595
596
597
598
class NonEquilibriumVolumeVariablesImplementation<Traits,
                                                  EquilibriumVolumeVariables,
                                                  true/*chemicalNonEquilibrium?*/,
                                                  true/*thermalNonEquilibrium?*/,
                                                  2>
: public EquilibriumVolumeVariables
599
{
600
    using ParentType = EquilibriumVolumeVariables;
601
    using FluidState = typename Traits::FluidState;
602
    using FS = typename Traits::FluidSystem;
603
604
605
606
607
608
609
    using ParameterCache = typename Traits::FluidSystem::ParameterCache;
    using Scalar = typename Traits::PrimaryVariables::value_type;

    using ModelTraits = typename Traits::ModelTraits;
    static constexpr auto numEnergyEqFluid = ModelTraits::numEnergyEqFluid();
    static constexpr auto numEnergyEqSolid = ModelTraits::numEnergyEqSolid();

610
611
612
613
614
    static constexpr auto phase0Idx = FS::phase0Idx;
    static constexpr auto phase1Idx = FS::phase1Idx;
    static constexpr auto sPhaseIdx = FS::numPhases;
    static constexpr auto wCompIdx = FS::comp0Idx;
    static constexpr auto nCompIdx = FS::comp1Idx;
615

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

619
620
621
622
    using NumFluidPhasesArray = std::array<Scalar, ModelTraits::numFluidPhases()>;
    using InterfacialAreasArray =  std::array<std::array<Scalar, ModelTraits::numFluidPhases()+numEnergyEqSolid>,
                                              ModelTraits::numFluidPhases()+numEnergyEqSolid>;

623
public:
624
625
    using Indices = typename ModelTraits::Indices;

626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
    /*!
     * \brief Update all quantities for a given control volume
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param problem The object specifying the problem which ought to
     *                be simulated
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    void update(const ElemSol &elemSol,
                const Problem &problem,
                const Element &element,
                const Scv& scv)
    {
        // Update parent type (also completes the fluid state)
        ParentType::update(elemSol, problem, element, scv);

        ParameterCache paramCache;
645
        paramCache.updateAll(this->fluidState());
646
        updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
647
        updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
648
649
    }

650
651
    /*!
     * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
652
653
654
655
656
657
658
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param fluidState Container for all the secondary variables concerning the fluids
     * \param paramCache The parameter cache corresponding to the fluid state
     * \param problem The problem to be solved
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
659
     */
660
    template<class ElemSol, class Problem, class Element, class Scv>
661
662
663
664
665
666
667
668
669
670
671
672
673
    void updateDimLessNumbers(const ElemSol& elemSol,
                              const FluidState& fluidState,
                              const ParameterCache& paramCache,
                              const Problem& problem,
                              const Element& element,
                              const Scv& scv)
    {
        factorMassTransfer_ = problem.spatialParams().factorMassTransfer(element, scv, elemSol);
        factorEnergyTransfer_ = problem.spatialParams().factorEnergyTransfer(element, scv, elemSol);
        characteristicLength_ = problem.spatialParams().characteristicLength(element, scv, elemSol);

        const auto vIdxGlobal = scv.dofIndex();
        using FluidSystem = typename Traits::FluidSystem;
674
        for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
675
        {
676
677
678
679
680
            const auto darcyMagVelocity = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, vIdxGlobal);
            const auto dynamicViscosity = fluidState.viscosity(phaseIdx);
            const auto density = fluidState.density(phaseIdx);
            const auto kinematicViscosity = dynamicViscosity/density;
            const auto heatCapacity = FluidSystem::heatCapacity(fluidState, paramCache, phaseIdx);
681
682
683
684
685
686
687
688
689
690
691
            const auto thermalConductivity = FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);

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

            reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_, kinematicViscosity);
692
693
694
695
696
697
            prandtlNumber_[phaseIdx] = DimLessNum::prandtlNumber(dynamicViscosity, heatCapacity, thermalConductivity);
            schmidtNumber_[phaseIdx] = DimLessNum::schmidtNumber(dynamicViscosity, density, diffCoeff);
            nusseltNumber_[phaseIdx] = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
                                                                       prandtlNumber_[phaseIdx],
                                                                       porosity,
                                                                       ModelTraits::nusseltFormulation());
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
            // If Diffusion is not enabled, Sherwood is divided by zero
            sherwoodNumber_[phaseIdx] = DimLessNum::sherwoodNumber(reynoldsNumber_[phaseIdx],
                                                                   schmidtNumber_[phaseIdx],
                                                                   ModelTraits::sherwoodFormulation());
        }
    }

     /*!
     * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param fluidState Container for all the secondary variables concerning the fluids
     * \param paramCache The parameter cache corresponding to the fluid state
     * \param problem The problem to be solved
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
     */
    template<class ElemSol, class Problem, class Element, class Scv>
716
717
718
719
720
721
    void updateInterfacialArea(const ElemSol& elemSol,
                               const FluidState& fluidState,
                               const ParameterCache& paramCache,
                               const Problem& problem,
                               const Element& element,
                               const Scv& scv)
722
723
    {
        // obtain (standard) material parameters (needed for the residual saturations)
724
        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
725
726

        //obtain parameters for interfacial area constitutive relations
727
        const auto& aWettingNonWettingSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element, scv, elemSol);
728

729
730
        const Scalar pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
        const Scalar Sw = fluidState.saturation(phase0Idx);
731
732

        using AwnSurface = typename Problem::SpatialParams::AwnSurface;
733
        const auto awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc);
734
735
736
        interfacialArea_[phase0Idx][phase1Idx] = awn;
        interfacialArea_[phase1Idx][phase0Idx] = interfacialArea_[phase0Idx][phase1Idx];
        interfacialArea_[phase0Idx][phase0Idx] = 0.;
737

738
        using AnsSurface = typename Problem::SpatialParams::AnsSurface;
739
740
        const auto& aNonWettingSolidSurfaceParams = problem.spatialParams().aNonWettingSolidSurfaceParams(element, scv, elemSol);
        const auto ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, Sw, pc);
741

742
743
        // 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.
744
745
746
747
748
749
750
751
752
753
754
755
756
757
        static const bool computeAwsFromAnsAndPcMax = getParam<bool>("SpatialParams.ComputeAwsFromAnsAndPcMax", true);
        if (computeAwsFromAnsAndPcMax)
        {
            // I know the solid surface from the pore network. But it is more consistent to use the fit value.
            const Scalar pcMax = aWettingNonWettingSurfaceParams.pcMax();
            const auto solidSurface = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax);
            interfacialArea_[phase0Idx][sPhaseIdx] = solidSurface - ans;
        }
        else
        {
            using AwsSurface = typename Problem::SpatialParams::AwsSurface;
            const auto& aWettingSolidSurfaceParams = problem.spatialParams().aWettingSolidSurfaceParams(element, scv, elemSol);
            interfacialArea_[phase0Idx][sPhaseIdx] = AwsSurface::interfacialArea(aWettingSolidSurfaceParams, materialParams, Sw, pc);
        }
758

759
        interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
760
        interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
761

762
763
764
        interfacialArea_[phase1Idx][sPhaseIdx] = ans;
        interfacialArea_[sPhaseIdx][phase1Idx] = interfacialArea_[phase1Idx][sPhaseIdx];
        interfacialArea_[phase1Idx][phase1Idx] = 0.;
765
766
767
768
    }

    /*!
     * \brief The specific interfacial area between two fluid phases [m^2 / m^3]
769
     * \note This is _only_ required by the kinetic mass/energy modules
770
771
772
773
774
775
776
777
778
     */
    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
779
780
781
    const Scalar reynoldsNumber(const unsigned int phaseIdx) const
    { return reynoldsNumber_[phaseIdx]; }

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

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

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

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

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

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

806
    //! access function pre factor mass transfer
807
808
    const Scalar factorMassTransfer() const
    { return factorMassTransfer_; }
809
810
811

private:
    //! dimensionless numbers
812
813
814
815
816
    NumFluidPhasesArray reynoldsNumber_;
    NumFluidPhasesArray prandtlNumber_;
    NumFluidPhasesArray nusseltNumber_;
    NumFluidPhasesArray schmidtNumber_;
    NumFluidPhasesArray sherwoodNumber_;
817
818
819
    Scalar characteristicLength_;
    Scalar factorEnergyTransfer_;
    Scalar factorMassTransfer_;
820
    InterfacialAreasArray interfacialArea_;
821
};
822

823
824
825
} // namespace Dumux

#endif