mpncvolumevariablesiakinetic.hh 29.4 KB
Newer Older
1
2
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
3
/*****************************************************************************
4
 *   See the file COPYING for full copying permissions.                      *
5
6
7
8
9
10
11
12
 *                                                                           *
 *   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          *
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
 *   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
 *
 * \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_MPNC_VOLUME_VARIABLES_IA_KINETIC_HH
#define DUMUX_MPNC_VOLUME_VARIABLES_IA_KINETIC_HH

32
#include <dune/common/version.hh>
33
#include <dumux/common/dimensionlessnumbers.hh>
34
35
#include "mpncvolumevariablesia.hh"
#include "mpncpropertieskinetic.hh"
36
37
38
39
40
41
42
43
44

namespace Dumux
{


////////////////////////////////////////////////////////////////////////////////////////////////////
// specialization for the case of kinetic mass AND energy transfer
////////////////////////////////////////////////////////////////////////////////////////////////////
template <class TypeTag, bool enableKinetic >
45
class MPNCVolumeVariablesIA<TypeTag, enableKinetic, /*numEnergyEqs=*/3>
46
47
48
49
50
51
52
53
54
55
{
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
    typedef typename FluidSystem::ParameterCache ParameterCache;
    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
56
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
57
    typedef typename GridView::template Codim<0>::Entity Element;
58
59
    typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;

60
61
62
63
64
65
66
67
68

    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 { dim = GridView::dimension};
Bernd Flemisch's avatar
Bernd Flemisch committed
69
    enum { dimWorld = GridView::dimensionworld};
70
    enum { numEnergyEqs = Indices::numPrimaryEnergyVars};
71
72
    enum { nusseltFormulation = GET_PROP_VALUE(TypeTag, NusseltFormulation)} ;
    enum { sherwoodFormulation = GET_PROP_VALUE(TypeTag, SherwoodFormulation)} ;
73
74
    enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion)} ;

75
76

    typedef DimensionlessNumbers<Scalar> DimLessNum ;
77

78
    typedef Dune::FieldVector<Scalar,dimWorld> GlobalPosition;
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


    typedef typename GET_PROP_TYPE(TypeTag, AwnSurface) AwnSurface;
    typedef typename AwnSurface::Params AwnSurfaceParams;

    typedef typename GET_PROP_TYPE(TypeTag, AwsSurface) AwsSurface;
    typedef typename AwsSurface::Params AwsSurfaceParams;

    typedef typename GET_PROP_TYPE(TypeTag, AnsSurface) AnsSurface;
    typedef typename AnsSurface::Params AnsSurfaceParams;


public:
    /*!
     * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
     */
    void update(const VolumeVariables & volVars,
                const FluidState & fluidState,
                const ParameterCache &paramCache,
                const PrimaryVariables &priVars,
                const Problem &problem,
                const Element & element,
                const FVElementGeometry & fvGeometry,
                const unsigned int scvIdx)
    {
104
        // obtain (standard) material parameters (needed for the residual saturations)
105
        const MaterialLawParams & materialParams  = problem.spatialParams().materialLawParams(element,fvGeometry,scvIdx);
106

107
        //obtain parameters for interfacial area constitutive relations
108
109
110
111
        const AwnSurfaceParams & aWettingNonWettingSurfaceParams
               = problem.spatialParams().aWettingNonWettingSurfaceParams(element,fvGeometry,scvIdx);
        const AnsSurfaceParams & aNonWettingSolidSurfaceParams
               = problem.spatialParams().aNonWettingSolidSurfaceParams(element,fvGeometry,scvIdx);
112
113
114
115

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

116
117
        const Scalar pc = fluidState.pressure(nPhaseIdx) - fluidState.pressure(wPhaseIdx);
        const Scalar Sw = fluidState.saturation(wPhaseIdx);
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
        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
140
        awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc ); // 10.; //
141
142
143
144
145

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

146
        Scalar ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams,Sw, pc ); // 10.; //
147
148
149
150
151
152
//        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
153
       const Scalar pcMax = problem.spatialParams().pcMax(element,
154
155
156
157
                                                            fvGeometry,
                                                            scvIdx);
        // I know the solid surface from the pore network. But it is more consistent to use the fit value.
        // solidSurface_   = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.specificSolidsurface);
158
        solidSurface_   = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax );
159
160
161
162
163
164
165
166
167
168
        Valgrind::CheckDefined(solidSurface_);

//        if (ans > solidSurface_){
//            const  GlobalPosition & globalPos =  fvGeometry.subContVol[scvIdx].global ;
//            std::stringstream positionString ;
//            positionString << " Here:";
//            for(int i=0; i<dim; i++)
//                positionString << " x"<< (i+1) << "="  << globalPos[i] << " "   ;
//            positionString << "\n";
//
169
170
//              std::cout<<"a_{ns} > a_s, set a_{ns}=" << ans <<" to a_{ns}=a_s="<<solidSurface_ << "
//                                    with S_w="<< Sw << " p_c= "<< pc <<  positionString.str() ;
171
//
172
//              ans = solidSurface_ ;
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
//        }

#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
188
189
        const AwsSurfaceParams & aWettingSolidSurfaceParams
               = problem.spatialParams().aWettingSolidSurfaceParams();
190
        Valgrind::CheckDefined(aWettingSolidSurfaceParams);
191
        const Scalar aws = AwsSurface::interfacialArea(aWettingSolidSurfaceParams,materialParams, Sw, pc ); // 10.; //
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
        interfacialArea_[wPhaseIdx][sPhaseIdx] = aws ;
        interfacialArea_[sPhaseIdx][wPhaseIdx] = interfacialArea_[wPhaseIdx][sPhaseIdx];
        interfacialArea_[sPhaseIdx][sPhaseIdx] = 0. ;
#endif

        Valgrind::CheckDefined(interfacialArea_);

        factorMassTransfer_   = problem.spatialParams().factorMassTransfer(element,
                                                                           fvGeometry,
                                                                           scvIdx);

        factorEnergyTransfer_   = problem.spatialParams().factorEnergyTransfer(element,
                                                                               fvGeometry,
                                                                               scvIdx);

        characteristicLength_   = problem.spatialParams().characteristicLength(element,
                                                                               fvGeometry,
                                                                               scvIdx);

        // setting the dimensionless numbers.
        // obtaining the respective quantities.
213
214
215
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
        const unsigned int globalVertexIdx = problem.vertexMapper().subIndex(element, scvIdx, dim);
#else
216
        const unsigned int globalVertexIdx = problem.vertexMapper().map(element, scvIdx, dim);
217
#endif
218
219
220
221
222
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            const Scalar darcyMagVelocity     = problem.model().volumeDarcyMagVelocity(phaseIdx, globalVertexIdx);
            const Scalar dynamicViscosity     = fluidState.viscosity(phaseIdx);
            const Scalar density              = fluidState.density(phaseIdx);
            const Scalar kinematicViscosity   = dynamicViscosity / density;
223
224
            const Scalar heatCapacity         = FluidSystem::heatCapacity(fluidState, paramCache, phaseIdx);
            const Scalar thermalConductivity  = FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);
225

226

227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
            // diffusion coefficient of non-wetting component in wetting phase
            const Scalar diffCoeff = volVars.diffCoeff(phaseIdx, wCompIdx, nCompIdx) ;
            const Scalar porosity = problem.spatialParams().porosity(element,
                                                                   fvGeometry,
                                                                   scvIdx);

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

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

            nusseltNumber_[phaseIdx]    = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
                                                                          prandtlNumber_[phaseIdx],
243
244
                                                                          porosity,
                                                                          nusseltFormulation);
245
246
247
248
249

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

250
            // If Diffusion is not enabled, Sherwood is divided by zero
251
            sherwoodNumber_[phaseIdx]   = DimLessNum::sherwoodNumber(reynoldsNumber_[phaseIdx],
252
253
                                                                      schmidtNumber_[phaseIdx],
                                                                      sherwoodFormulation);
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
283
284
285
286
287
288
289
290
291
292
293
294
295
296
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
        }
    }

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


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
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
// specialization for the case of NO kinetic mass but kinteic energy transfer of a fluid mixture and solid
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <class TypeTag>
class MPNCVolumeVariablesIA<TypeTag, /*enableKinetic=*/false, /*numEnergyEqs=*/2>
{
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
    typedef typename FluidSystem::ParameterCache ParameterCache;
    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
    typedef typename GridView::template Codim<0>::Entity Element;

    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
    enum { dim = GridView::dimension};
    enum { dimWorld = GridView::dimensionworld};
    enum { numEnergyEqs = Indices::numPrimaryEnergyVars};
    enum { nusseltFormulation = GET_PROP_VALUE(TypeTag, NusseltFormulation)} ;

    typedef DimensionlessNumbers<Scalar> DimLessNum ;
public:
    /*!
     * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
     */
    void update(const VolumeVariables & volVars,
                const FluidState & fluidState,
                const ParameterCache &paramCache,
                const PrimaryVariables &priVars,
                const Problem &problem,
                const Element & element,
                const FVElementGeometry & fvGeometry,
                const unsigned int scvIdx)
    {
        factorMassTransfer_   = problem.spatialParams().factorMassTransfer(element,
                                                                           fvGeometry,
                                                                           scvIdx);

        factorEnergyTransfer_   = problem.spatialParams().factorEnergyTransfer(element,
                                                                               fvGeometry,
                                                                               scvIdx);

        characteristicLength_   = problem.spatialParams().characteristicLength(element,
                                                                               fvGeometry,
                                                                               scvIdx);

        // setting the dimensionless numbers.
        // obtaining the respective quantities.
388
389
390
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
        const unsigned int globalVertexIdx = problem.vertexMapper().subIndex(element, scvIdx, dim);
#else
391
        const unsigned int globalVertexIdx = problem.vertexMapper().map(element, scvIdx, dim);
392
#endif
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
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            const Scalar darcyMagVelocity     = problem.model().volumeDarcyMagVelocity(phaseIdx, globalVertexIdx);
            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,
                                                                   fvGeometry,
                                                                   scvIdx);

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

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


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


    //! 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_);
#endif
    }

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


478
479
480
481
////////////////////////////////////////////////////////////////////////////////////////////////////
// specialization for the case of (only) kinetic mass transfer
////////////////////////////////////////////////////////////////////////////////////////////////////
template <class TypeTag>
482
class MPNCVolumeVariablesIA<TypeTag, /*enableKinetic=*/true, /*numEnergyEqs=*/0>
483
484
485
486
487
488
489
490
491
492
493
{
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
    typedef typename FluidSystem::ParameterCache ParameterCache;
    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
    typedef typename GridView::template Codim<0>::Entity Element;
494
    typedef DimensionlessNumbers<Scalar> DimLessNum;
495
496
497
498
499
500
501
502
503
504

    typedef typename GET_PROP_TYPE(TypeTag, AwnSurface) AwnSurface;
    typedef typename AwnSurface::Params AwnSurfaceParams;

    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
    enum { wPhaseIdx = FluidSystem::wPhaseIdx };
    enum { nPhaseIdx = FluidSystem::nPhaseIdx };
    enum { wCompIdx  = FluidSystem::wCompIdx };
    enum { nCompIdx  = FluidSystem::nCompIdx };
    enum { dim       = GridView::dimension};
505
    enum { dimWorld  = GridView::dimensionworld};
506
    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
507
    enum { sherwoodFormulation = GET_PROP_VALUE(TypeTag, SherwoodFormulation)} ;
508
    typedef Dune::FieldVector<Scalar,dimWorld> GlobalPosition;
509
510
    typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;

511
512
513
514
515
516
517
518
519
520
521
522
523
524
525

public:
    /*!
     * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
     */
    void update(const VolumeVariables & volVars,
                const FluidState & fluidState,
                const ParameterCache & paramCache,
                const PrimaryVariables &priVars,
                const Problem &problem,
                const Element & element,
                const FVElementGeometry & fvGeometry,
                const unsigned int scvIdx)
    {
        //obtain parameters for awnsurface
526
527
528
529
        const AwnSurfaceParams & awnSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element,fvGeometry,scvIdx) ;

        // obtain (standard) material parameters (needed for the residual saturations)
        const MaterialLawParams & materialParams  = problem.spatialParams().materialLawParams(element,fvGeometry,scvIdx) ;
530
531
532
533
534
535

        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
536
        interfacialArea_ = AwnSurface::interfacialArea(awnSurfaceParams, materialParams, Sw, pc );
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

        Valgrind::CheckDefined(interfacialArea_);

        factorMassTransfer_   = problem.spatialParams().factorMassTransfer(element,
                                                                           fvGeometry,
                                                                           scvIdx);

        characteristicLength_   = problem.spatialParams().characteristicLength(element,
                                                                               fvGeometry,
                                                                               scvIdx);
        // setting the dimensionless numbers.
        // obtaining the respective quantities.
        int globalVertexIdx = problem.vertexMapper().map(element, scvIdx, dim);
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            const Scalar darcyMagVelocity     = problem.model().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 = volVars.diffCoeff(phaseIdx, wCompIdx, nCompIdx) ;

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

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

            sherwoodNumber_[phaseIdx]   = DimLessNum::sherwoodNumber(reynoldsNumber_[phaseIdx],
568
569
                                                                      schmidtNumber_[phaseIdx],
                                                                      sherwoodFormulation);
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
        }
    }

    /*!
     * \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 and phaseJIdx == wPhaseIdx) or (phaseIIdx == wPhaseIdx and 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] ;
};



} // namespace Dumux

#endif