2p2cvolumevariables.hh 27.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
// -*- 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
 *
 * \brief Contains the quantities which are constant within a
Beatrix Becker's avatar
Beatrix Becker committed
23
 *        finite volume in the two-phase two-component model.
24
25
26
27
 */
#ifndef DUMUX_2P2C_VOLUME_VARIABLES_HH
#define DUMUX_2P2C_VOLUME_VARIABLES_HH

28
#include <dumux/implicit/common/implicitmodel.hh>
29
30
31
#include <dumux/material/fluidstates/compositionalfluidstate.hh>
#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
32
33
#include "2p2cproperties.hh"
#include "2p2cindices.hh"
34
35
36
37
38
39

namespace Dumux
{

/*!
 * \ingroup TwoPTwoCModel
40
 * \ingroup ImplicitVolumeVariables
Beatrix Becker's avatar
Beatrix Becker committed
41
42
 * \brief Contains the quantities which are constant within a
 *        finite volume in the two-phase two-component model.
43
44
 */
template <class TypeTag>
45
class TwoPTwoCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
46
{
47
    typedef ImplicitVolumeVariables<TypeTag> ParentType;
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) Implementation;
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
    typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
    enum {
        numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
        numComponents = GET_PROP_VALUE(TypeTag, NumComponents)
    };

    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
    enum {
        wCompIdx = Indices::wCompIdx,
        nCompIdx = Indices::nCompIdx,
        wPhaseIdx = Indices::wPhaseIdx,
        nPhaseIdx = Indices::nPhaseIdx
    };

    // present phases
    enum {
        wPhaseOnly = Indices::wPhaseOnly,
        nPhaseOnly = Indices::nPhaseOnly,
        bothPhases = Indices::bothPhases
    };

    // formulations
    enum {
        formulation = GET_PROP_VALUE(TypeTag, Formulation),
78
79
        pwsn = TwoPTwoCFormulation::pwsn,
        pnsw = TwoPTwoCFormulation::pnsw
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
    };

    // primary variable indices
    enum {
        switchIdx = Indices::switchIdx,
        pressureIdx = Indices::pressureIdx
    };

    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GridView::template Codim<0>::Entity Element;
    enum { dim = GridView::dimension};

    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem> MiscibleMultiPhaseComposition;
95
    static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
96
97
98
    static const bool useConstraintSolver = GET_PROP_VALUE(TypeTag, UseConstraintSolver);
    static_assert(useMoles || (!useMoles && useConstraintSolver),
                  "if UseMoles is set false, UseConstraintSolver has to be set to true");
99
100
    typedef Dumux::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;

101
    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
102
    enum { dofCodim = isBox ? dim : 0 };
103

104
105
106
107
108
public:
    //! The type of the object returned by the fluidState() method
    typedef Dumux::CompositionalFluidState<Scalar, FluidSystem> FluidState;

    /*!
109
     * \copydoc ImplicitVolumeVariables::update
110
111
112
113
114
115
116
117
118
119
120
121
122
123
     */
    void update(const PrimaryVariables &priVars,
                const Problem &problem,
                const Element &element,
                const FVElementGeometry &fvGeometry,
                const int scvIdx,
                const bool isOldSol)
    {
        ParentType::update(priVars,
                           problem,
                           element,
                           fvGeometry,
                           scvIdx,
                           isOldSol);
124
        
125
        completeFluidState(priVars, problem, element, fvGeometry, scvIdx, fluidState_, isOldSol);
126
        
127
128
129
130
        /////////////
        // calculate the remaining quantities
        /////////////
        const MaterialLawParams &materialParams =
131
132
        problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
        
133
134
135
136
137
138
        // Second instance of a parameter cache.
        // Could be avoided if diffusion coefficients also
        // became part of the fluid state.
        typename FluidSystem::ParameterCache paramCache;
        paramCache.updateAll(fluidState_);
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
139
140
            
            
141
142
143
144
145
            // relative permeabilities
            Scalar kr;
            if (phaseIdx == wPhaseIdx)
                kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx));
            else // ATTENTION: krn requires the wetting phase saturation
146
                // as parameter!
147
148
149
                kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx));
            relativePermeability_[phaseIdx] = kr;
            Valgrind::CheckDefined(relativePermeability_[phaseIdx]);
150
            
151
152
            // binary diffusion coefficients
            diffCoeff_[phaseIdx] =
153
154
155
156
157
            FluidSystem::binaryDiffusionCoefficient(fluidState_,
                                                    paramCache,
                                                    phaseIdx,
                                                    wCompIdx,
                                                    nCompIdx);
158
159
            Valgrind::CheckDefined(diffCoeff_[phaseIdx]);
        }
160
        
161
162
163
164
165
        // porosity
        porosity_ = problem.spatialParams().porosity(element,
                                                     fvGeometry,
                                                     scvIdx);
        Valgrind::CheckDefined(porosity_);
166
        
167
168
169
170
171
        // energy related quantities not contained in the fluid state
        asImp_().updateEnergy_(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
    }

    /*!
Bernd Flemisch's avatar
Bernd Flemisch committed
172
     * \copydoc ImplicitModel::completeFluidState
173
174
175
176
177
178
179
180
181
182
183
184
185
     * \param isOldSol Specifies whether this is the previous solution or the current one
     */
    static void completeFluidState(const PrimaryVariables& priVars,
                                   const Problem& problem,
                                   const Element& element,
                                   const FVElementGeometry& fvGeometry,
                                   int scvIdx,
                                   FluidState& fluidState,
                                   bool isOldSol = false)
    {
        Scalar t = Implementation::temperature_(priVars, problem, element,
                                                fvGeometry, scvIdx);
        fluidState.setTemperature(t);
186
        
187
        int dofIdxGlobal = problem.model().dofMapper().map(element, scvIdx, dofCodim);
188
        
189
        int phasePresence = problem.model().phasePresence(dofIdxGlobal, isOldSol);
190
        
191
192
193
        /////////////
        // set the saturations
        /////////////
194
        Scalar sn;
195
        if (phasePresence == nPhaseOnly)
196
            sn = 1.0;
197
        else if (phasePresence == wPhaseOnly) {
198
            sn = 0.0;
199
200
        }
        else if (phasePresence == bothPhases) {
201
            if (formulation == pwsn)
202
                sn = priVars[switchIdx];
203
            else if (formulation == pnsw)
204
                sn = 1.0 - priVars[switchIdx];
205
206
207
            else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
        }
        else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
208
209
        fluidState.setSaturation(wPhaseIdx, 1 - sn);
        fluidState.setSaturation(nPhaseIdx, sn);
210
        
211
212
213
        /////////////
        // set the pressures of the fluid phases
        /////////////
214
        
215
216
        // calculate capillary pressure
        const MaterialLawParams &materialParams =
217
        problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
218
        Scalar pc = MaterialLaw::pc(materialParams, 1 - sn);
219
        
220
        if (formulation == pwsn) {
221
            fluidState.setPressure(wPhaseIdx, priVars[pressureIdx]);
222
            fluidState.setPressure(nPhaseIdx, priVars[pressureIdx] + pc);
223
        }
224
        else if (formulation == pnsw) {
225
            fluidState.setPressure(nPhaseIdx, priVars[pressureIdx]);
226
            fluidState.setPressure(wPhaseIdx, priVars[pressureIdx] - pc);
227
228
        }
        else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
229
        
230
231
232
233
        /////////////
        // calculate the phase compositions
        /////////////
        typename FluidSystem::ParameterCache paramCache;
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
        
        //get the phase pressures and set the fugacity coefficients here if constraintsolver is not used
        Scalar pn = 0;
        Scalar pw = 0;
        
        if(!useConstraintSolver) {
            if (formulation == pwsn) {
                pw = priVars[pressureIdx];
                pn = pw + pc;
            }
            else {
                pn = priVars[pressureIdx];
                pw = pn - pc;
            }
            
            for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
                assert(FluidSystem::isIdealMixture(phaseIdx));
                
                for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
                    Scalar phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
                    fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
                }
            }
        }
        
259
260
261
262
263
264
        // now comes the tricky part: calculate phase compositions
        if (phasePresence == bothPhases) {
            // both phases are present, phase compositions are a
            // result of the the nonwetting <-> wetting equilibrium. This is
            // the job of the "MiscibleMultiPhaseComposition"
            // constraint solver
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
            if(useConstraintSolver) {
                MiscibleMultiPhaseComposition::solve(fluidState,
                                                     paramCache,
                                                     /*setViscosity=*/true,
                                                     /*setInternalEnergy=*/false);
            }  
            // ... or calculated explicitly this way ...
            else {  
                //get the partial pressure of the main component of the the wetting phase ("H20") within the nonwetting (gas) phase == vapor pressure due to equilibrium
                //note that in this case the fugacityCoefficient * pw is the vapor pressure (see implementation in respective fluidsystem)
                Scalar partPressH2O = FluidSystem::fugacityCoefficient(fluidState, 
                                                                       wPhaseIdx, 
                                                                       wCompIdx) * pw;
                                                                       
                // get the partial pressure of the main component of the the nonwetting (gas) phase ("Air") 
                Scalar partPressAir = pn - partPressH2O;
                
                //calculate the mole fractions of the components within the nonwetting phase
                Scalar xnn = partPressAir/pn;
                Scalar xnw = partPressH2O/pn;
                
                // calculate the mole fractions of the components within the wetting phase
                //note that in this case the fugacityCoefficient * pw is the Henry Coefficient (see implementation in respective fluidsystem)
                Scalar xwn = partPressAir
                  / (FluidSystem::fugacityCoefficient(fluidState,
                                                      wPhaseIdx,nCompIdx) 
                  * pw);
                
                Scalar xww = 1.0 -xwn;
                
                //set all mole fractions
                fluidState.setMoleFraction(wPhaseIdx, wCompIdx, xww);
                fluidState.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
                fluidState.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
                fluidState.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
                
                paramCache.updateComposition(fluidState, wPhaseIdx); 
                paramCache.updateComposition(fluidState, nPhaseIdx);
                
                //set the phase densities
                Scalar rhoW = FluidSystem::density(fluidState, paramCache, wPhaseIdx);
                Scalar rhoN = FluidSystem::density(fluidState, paramCache, nPhaseIdx);
                
                fluidState.setDensity(wPhaseIdx, rhoW);
                fluidState.setDensity(nPhaseIdx, rhoN);     
            }  
311
312
313
314
        }
        else if (phasePresence == nPhaseOnly) {
            // only the nonwetting phase is present, i.e. nonwetting phase
            // composition is stored explicitly.
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
            
            if(!useMoles) //mass-fraction formulation
            {
                // extract _mass_ fractions in the nonwetting phase
                Scalar massFractionN[numComponents];
                massFractionN[wCompIdx] = priVars[switchIdx];
                massFractionN[nCompIdx] = 1 - massFractionN[wCompIdx];
                
                // calculate average molar mass of the nonwetting phase
                Scalar M1 = FluidSystem::molarMass(wCompIdx);
                Scalar M2 = FluidSystem::molarMass(nCompIdx);
                Scalar X2 = massFractionN[nCompIdx];
                Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
                
                // convert mass to mole fractions and set the fluid state
                fluidState.setMoleFraction(nPhaseIdx, wCompIdx, massFractionN[wCompIdx]*avgMolarMass/M1);
                fluidState.setMoleFraction(nPhaseIdx, nCompIdx, massFractionN[nCompIdx]*avgMolarMass/M2);
            }
            else //mole-fraction formulation
            {
335
336
337
                Scalar moleFractionN[numComponents];
                moleFractionN[wCompIdx] = priVars[switchIdx];
                moleFractionN[nCompIdx] = 1 - moleFractionN[wCompIdx];
338
                
339
340
341
                // set the fluid state
                fluidState.setMoleFraction(nPhaseIdx, wCompIdx, moleFractionN[wCompIdx]);
                fluidState.setMoleFraction(nPhaseIdx, nCompIdx, moleFractionN[nCompIdx]);
342
            }
343
            // calculate the composition of the remaining phases (as
Beatrix Becker's avatar
Beatrix Becker committed
344
            // well as the densities of all phases). This is the job
345
            // of the "ComputeFromReferencePhase" constraint solver
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
            if (useConstraintSolver) {
                ComputeFromReferencePhase::solve(fluidState,
                                                 paramCache,
                                                 nPhaseIdx,
                                                 /*setViscosity=*/true,
                                                 /*setInternalEnergy=*/false);
            }
            // ... or calculated explicitly this way ...
            else {    
                // note that the water phase is actually not existing!
                // thus, this is used as phase switch criterion
                Scalar xnw = priVars[switchIdx];
                Scalar xnn = 1.0 -xnw;
                
                //first, xww: 
                // xnw * pn = "actual" (hypothetical) vapor pressure 
                // fugacityCoefficient * pw = vapor pressure given by thermodynamic conditions
                // Here, xww is not actually the mole fraction of water in the wetting phase
                // xww is only the ratio of "actual" vapor pressure / "thermodynamic" vapor pressure  
                // If xww > 1 : gas is over-saturated with water vapor, 
                // condensation takes place (see switch criterion in model) 
                Scalar xww = xnw * pn 
                  / (FluidSystem::fugacityCoefficient(fluidState,
                                                      wPhaseIdx,wCompIdx)
                     * pw);
                
                // now, xwn: 	     
                //partialPressure / xwn = Henry
                //partialPressure = xnn * pn	
                //xwn = xnn * pn / Henry
                // Henry = fugacityCoefficient * pw
                Scalar xwn = xnn * pn / (FluidSystem::fugacityCoefficient(fluidState,
                                                                          wPhaseIdx,nCompIdx)
                                         * pw); 
                
                fluidState.setMoleFraction(wPhaseIdx, wCompIdx, xww);
                fluidState.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
                
                paramCache.updateComposition(fluidState, wPhaseIdx); 
                paramCache.updateComposition(fluidState, nPhaseIdx);
                
                Scalar rhoW = FluidSystem::density(fluidState, paramCache, wPhaseIdx);
                Scalar rhoN = FluidSystem::density(fluidState, paramCache, nPhaseIdx);
                
                fluidState.setDensity(wPhaseIdx, rhoW);
                fluidState.setDensity(nPhaseIdx, rhoN);
            } 
393
394
395
396
        }
        else if (phasePresence == wPhaseOnly) {
            // only the wetting phase is present, i.e. wetting phase
            // composition is stored explicitly.
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
            if(!useMoles) //mass-fraction formulation
            {
                // extract _mass_ fractions in the nonwetting phase
                Scalar massFractionW[numComponents];
                massFractionW[nCompIdx] = priVars[switchIdx];
                massFractionW[wCompIdx] = 1 - massFractionW[nCompIdx];
                
                // calculate average molar mass of the nonwetting phase
                Scalar M1 = FluidSystem::molarMass(wCompIdx);
                Scalar M2 = FluidSystem::molarMass(nCompIdx);
                Scalar X2 = massFractionW[nCompIdx];
                Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
                
                // convert mass to mole fractions and set the fluid state
                fluidState.setMoleFraction(wPhaseIdx, wCompIdx, massFractionW[wCompIdx]*avgMolarMass/M1);
                fluidState.setMoleFraction(wPhaseIdx, nCompIdx, massFractionW[nCompIdx]*avgMolarMass/M2);
            }
            else //mole-fraction formulation
            {
416
417
418
                Scalar moleFractionW[numComponents];
                moleFractionW[nCompIdx] = priVars[switchIdx];
                moleFractionW[wCompIdx] = 1 - moleFractionW[nCompIdx];
419
                
420
421
422
                // set the fluid state
                fluidState.setMoleFraction(wPhaseIdx, wCompIdx, moleFractionW[wCompIdx]);
                fluidState.setMoleFraction(wPhaseIdx, nCompIdx, moleFractionW[nCompIdx]);
423
            }
424
            // calculate the composition of the remaining phases (as
Beatrix Becker's avatar
Beatrix Becker committed
425
            // well as the densities of all phases). This is the job
426
            // of the "ComputeFromReferencePhase" constraint solver
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
            if (useConstraintSolver) {
                ComputeFromReferencePhase::solve(fluidState,
                                                 paramCache,
                                                 wPhaseIdx,
                                                 /*setViscosity=*/true,
                                                 /*setInternalEnergy=*/false);
            }
            // ... or calculated explicitly this way ...
            else {    
                // note that the gas phase is actually not existing!
                // thus, this is used as phase switch criterion
                Scalar xwn = priVars[switchIdx]; 
                
                //first, xnw:
                //psteam = xnw * pn = partial pressure of water in gas phase
                //psteam = fugacityCoefficient * pw
                Scalar xnw = (FluidSystem::fugacityCoefficient(fluidState,
                                                               wPhaseIdx,wCompIdx)
                              * pw) / pn ;
                
                //now, xnn: 	     
                // xwn = partialPressure / Henry
                // partialPressure = pn * xnn
                // xwn = pn * xnn / Henry
                // xnn = xwn * Henry / pn
                // Henry = fugacityCoefficient * pw
                Scalar xnn = xwn * (FluidSystem::fugacityCoefficient(fluidState,
                                                                     wPhaseIdx,nCompIdx)
                                    * pw) / pn ;
                
                fluidState.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
                fluidState.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
                
                paramCache.updateComposition(fluidState, wPhaseIdx); 
                paramCache.updateComposition(fluidState, nPhaseIdx);
                
                Scalar rhoW = FluidSystem::density(fluidState, paramCache, wPhaseIdx);
                Scalar rhoN = FluidSystem::density(fluidState, paramCache, nPhaseIdx);
                
                fluidState.setDensity(wPhaseIdx, rhoW);
                fluidState.setDensity(nPhaseIdx, rhoN);
            }
469
        }
470
        
471
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
472
473
474
475
476
477
478
479
            //set the viscosity here if constraintsolver is not used
            if(!useConstraintSolver) {
                const Scalar mu =
                FluidSystem::viscosity(fluidState,
                                       paramCache,
                                       phaseIdx);
                fluidState.setViscosity(phaseIdx,mu);
            }	  
480
481
482
483
            // compute and set the enthalpy
            Scalar h = Implementation::enthalpy_(fluidState, paramCache, phaseIdx);
            fluidState.setEnthalpy(phaseIdx, h);
        }
484
485
   }
    
486
487

    /*!
488
     * \brief Returns the phase state within the control volume.
489
490
491
492
493
     */
    const FluidState &fluidState() const
    { return fluidState_; }

    /*!
494
     * \brief Returns the saturation of a given phase within
495
     *        the control volume in \f$[-]\f$.
496
497
498
499
500
501
     *
     * \param phaseIdx The phase index
     */
    Scalar saturation(const int phaseIdx) const
    { return fluidState_.saturation(phaseIdx); }

502
503
    /*!
     * \brief Returns the mass fraction of a given component in a
504
     *        given phase within the control volume in \f$[-]\f$.
505
506
507
508
509
510
511
512
     *
     * \param phaseIdx The phase index
     * \param compIdx The component index
     */
    Scalar massFraction(const int phaseIdx, const int compIdx) const
    { return fluidState_.massFraction(phaseIdx, compIdx); }

    /*!
513
     * \brief Returns the mole fraction of a given component in a
514
     *        given phase within the control volume in \f$[-]\f$.
515
516
517
518
519
520
521
     *
     * \param phaseIdx The phase index
     * \param compIdx The component index
     */
    Scalar moleFraction(const int phaseIdx, const int compIdx) const
    { return fluidState_.moleFraction(phaseIdx, compIdx); }

522
523
    /*!
     * \brief Returns the mass density of a given phase within the
524
     *        control volume in \f$[kg/m^3]\f$.
525
526
527
528
529
530
531
532
     *
     * \param phaseIdx The phase index
     */
    Scalar density(const int phaseIdx) const
    { return fluidState_.density(phaseIdx); }

    /*!
     * \brief Returns the mass density of a given phase within the
533
     *        control volume in \f$[mol/m^3]\f$.
534
535
536
537
538
539
540
541
     *
     * \param phaseIdx The phase index
     */
    Scalar molarDensity(const int phaseIdx) const
    { return fluidState_.density(phaseIdx) / fluidState_.averageMolarMass(phaseIdx); }

    /*!
     * \brief Returns the effective pressure of a given phase within
542
     *        the control volume in \f$[kg/(m*s^2)=N/m^2=Pa]\f$.
543
544
545
546
547
548
549
     *
     * \param phaseIdx The phase index
     */
    Scalar pressure(const int phaseIdx) const
    { return fluidState_.pressure(phaseIdx); }

    /*!
550
     * \brief Returns temperature within the control volume in \f$[K]\f$.
551
552
553
554
555
556
557
558
559
560
     *
     * Note that we assume thermodynamic equilibrium, i.e. the
     * temperature of the rock matrix and of all fluid phases are
     * identical.
     */
    Scalar temperature() const
    { return fluidState_.temperature(/*phaseIdx=*/0); }

    /*!
     * \brief Returns the relative permeability of a given phase within
561
     *        the control volume in \f$[-]\f$.
562
563
564
565
566
567
568
569
570
571
     *
     * \param phaseIdx The phase index
     */
    Scalar relativePermeability(const int phaseIdx) const
    {
        return relativePermeability_[phaseIdx];
    }

    /*!
     * \brief Returns the effective mobility of a given phase within
572
     *        the control volume in \f$[s*m/kg]\f$.
573
574
575
576
577
578
579
580
581
     *
     * \param phaseIdx The phase index
     */
    Scalar mobility(const int phaseIdx) const
    {
        return relativePermeability_[phaseIdx]/fluidState_.viscosity(phaseIdx);
    }

    /*!
582
583
     * \brief Returns the effective capillary pressure within the control volume
     *        in \f$[kg/(m*s^2)=N/m^2=Pa]\f$.
584
585
586
587
588
     */
    Scalar capillaryPressure() const
    { return fluidState_.pressure(nPhaseIdx) - fluidState_.pressure(wPhaseIdx); }

    /*!
589
     * \brief Returns the average porosity within the control volume in \f$[-]\f$.
590
591
592
593
594
     */
    Scalar porosity() const
    { return porosity_; }

    /*!
595
     * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
596
     */
597
    Scalar diffCoeff(const int phaseIdx) const
598
599
600
601
602
603
604
605
606
607
    { return diffCoeff_[phaseIdx]; }


protected:
    static Scalar temperature_(const PrimaryVariables &priVars,
                               const Problem& problem,
                               const Element &element,
                               const FVElementGeometry &fvGeometry,
                               int scvIdx)
    {
608
        return problem.temperatureAtPos(fvGeometry.subContVol[scvIdx].global);
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
    }

    template<class ParameterCache>
    static Scalar enthalpy_(const FluidState& fluidState,
                            const ParameterCache& paramCache,
                            const int phaseIdx)
    {
        return 0;
    }

    /*!
     * \brief Called by update() to compute the energy related quantities
     */
    void updateEnergy_(const PrimaryVariables &sol,
                       const Problem &problem,
                       const Element &element,
                       const FVElementGeometry &fvGeometry,
626
                       const int scvIdx,
627
628
629
                       bool isOldSol)
    { }

630
631
    Scalar porosity_; //!< Effective porosity within the control volume
    Scalar relativePermeability_[numPhases]; //!< Relative permeability within the control volume
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
    Scalar diffCoeff_[numPhases]; //!< Binary diffusion coefficients for the phases
    FluidState fluidState_;

private:
    Implementation &asImp_()
    { return *static_cast<Implementation*>(this); }

    const Implementation &asImp_() const
    { return *static_cast<const Implementation*>(this); }


};

} // end namespace

#endif