volumevariables.hh 24 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
21
 * \ingroup TwoPTwoCModel
22
23
24
25
26
27
 * \brief Contains the quantities which are constant within a
 *        finite volume in the two-phase two-component model.
 */
#ifndef DUMUX_2P2C_VOLUME_VARIABLES_HH
#define DUMUX_2P2C_VOLUME_VARIABLES_HH

28
#include <dumux/material/fluidstates/compositional.hh>
29
30
#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
Timo Koch's avatar
Timo Koch committed
31
#include <dumux/porousmediumflow/volumevariables.hh>
32
#include <dumux/discretization/methods.hh>
33

Timo Koch's avatar
Timo Koch committed
34
35
36
#include "indices.hh" // for formulation

namespace Dumux {
37
38
39
40
41
42
43

/*!
 * \ingroup TwoPTwoCModel
 * \brief Contains the quantities which are constant within a
 *        finite volume in the two-phase two-component model.
 */
template <class TypeTag>
Timo Koch's avatar
Timo Koch committed
44
class TwoPTwoCVolumeVariables : public PorousMediumFlowVolumeVariables<TypeTag>
45
{
Timo Koch's avatar
Timo Koch committed
46
    using ParentType = PorousMediumFlowVolumeVariables<TypeTag>;
47
    using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
48
49

    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
50
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
51
52
    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
53
    using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
54
55
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
    using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
56
57
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
58
    using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
59

60
61
62
    // component indices
    enum
    {
63
64
65
66
67
68
69
        wCompIdx = Indices::wCompIdx,
        nCompIdx = Indices::nCompIdx,
        wPhaseIdx = Indices::wPhaseIdx,
        nPhaseIdx = Indices::nPhaseIdx
    };

    // present phases
70
71
    enum
    {
72
73
74
75
76
77
        wPhaseOnly = Indices::wPhaseOnly,
        nPhaseOnly = Indices::nPhaseOnly,
        bothPhases = Indices::bothPhases
    };

    // formulations
78
79
    enum
    {
80
81
82
83
84
85
        formulation = GET_PROP_VALUE(TypeTag, Formulation),
        pwsn = TwoPTwoCFormulation::pwsn,
        pnsw = TwoPTwoCFormulation::pnsw
    };

    // primary variable indices
86
87
    enum
    {
88
89
90
91
        switchIdx = Indices::switchIdx,
        pressureIdx = Indices::pressureIdx
    };

92
    using Element = typename GridView::template Codim<0>::Entity;
93
94
95
96
97
98
    using PermeabilityType = typename SpatialParams::PermeabilityType;
    using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, FluidSystem>;

    static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
    static constexpr bool useKelvinEquation = GET_PROP_VALUE(TypeTag, UseKelvinEquation);
    static constexpr bool useConstraintSolver = GET_PROP_VALUE(TypeTag, UseConstraintSolver);
99
100
101
    static_assert(useMoles || (!useMoles && useConstraintSolver),
                  "if UseMoles is set false, UseConstraintSolver has to be set to true");

102
103
    using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem, useKelvinEquation>;

104
105
    static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
    static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
106
public:
Timo Koch's avatar
Timo Koch committed
107
    //! The type of the object returned by the fluidState() method
108
    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
109

Timo Koch's avatar
Timo Koch committed
110
111
112
113
    //! Pull member functions of the parent for deriving classes
    using ParentType::temperature;
    using ParentType::enthalpy;

114
    /*!
115
116
117
118
119
120
121
122
     * \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
    */
123
    void update(const ElementSolution &elemSol,
124
125
                const Problem &problem,
                const Element &element,
126
                const SubControlVolume& scv)
127
    {
128
        ParentType::update(elemSol, problem, element, scv);
129

Timo Koch's avatar
Timo Koch committed
130
        Implementation::completeFluidState(elemSol, problem, element, scv, fluidState_);
131
132
133
134

        /////////////
        // calculate the remaining quantities
        /////////////
135
        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
136
137
138
139
140
141

        // 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_);
142
143
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
        {
144
145
146
147
148
149
150
151
152
153
154
            // relative permeabilities
            Scalar kr;
            if (phaseIdx == wPhaseIdx)
                kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx));
            else // ATTENTION: krn requires the wetting phase saturation
                // as parameter!
                kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx));
            relativePermeability_[phaseIdx] = kr;
            Valgrind::CheckDefined(relativePermeability_[phaseIdx]);

            // binary diffusion coefficients
155
156
157
158
159
            diffCoeff_[phaseIdx] = FluidSystem::binaryDiffusionCoefficient(fluidState_,
                                                                           paramCache,
                                                                           phaseIdx,
                                                                           wCompIdx,
                                                                           nCompIdx);
160
161
        }

162
163
164
        // porosity & permeabilty
        porosity_ = problem.spatialParams().porosity(element, scv, elemSol);
        permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
165
166
167
    }

    /*!
168
169
170
171
172
173
174
175
176
     * \brief Complete the fluid state
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param problem The problem
     * \param element The element
     * \param scv The sub control volume
     * \param fluidState The fluid state
     *
     * Set temperature, saturations, capillary pressures, viscosities, densities and enthalpies.
177
     */
178
    static void completeFluidState(const ElementSolution& elemSol,
179
180
                                   const Problem& problem,
                                   const Element& element,
181
182
                                   const SubControlVolume& scv,
                                   FluidState& fluidState)
183
    {
184
        Scalar t = ParentType::temperature(elemSol, problem, element, scv);
185
        fluidState.setTemperature(t);
186

187
        const auto& priVars = ParentType::extractDofPriVars(elemSol, scv);
188
        const auto phasePresence = priVars.state();
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214

        /////////////
        // set the saturations
        /////////////
        Scalar sn;
        if (phasePresence == nPhaseOnly)
            sn = 1.0;
        else if (phasePresence == wPhaseOnly) {
            sn = 0.0;
        }
        else if (phasePresence == bothPhases) {
            if (formulation == pwsn)
                sn = priVars[switchIdx];
            else if (formulation == pnsw)
                sn = 1.0 - priVars[switchIdx];
            else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
        }
        else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
        fluidState.setSaturation(wPhaseIdx, 1 - sn);
        fluidState.setSaturation(nPhaseIdx, sn);

        /////////////
        // set the pressures of the fluid phases
        /////////////

        // calculate capillary pressure
215
        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
        Scalar pc = MaterialLaw::pc(materialParams, 1 - sn);

        if (formulation == pwsn) {
            fluidState.setPressure(wPhaseIdx, priVars[pressureIdx]);
            fluidState.setPressure(nPhaseIdx, priVars[pressureIdx] + pc);
        }
        else if (formulation == pnsw) {
            fluidState.setPressure(nPhaseIdx, priVars[pressureIdx]);
            fluidState.setPressure(wPhaseIdx, priVars[pressureIdx] - pc);
        }
        else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");

        /////////////
        // calculate the phase compositions
        /////////////
        typename FluidSystem::ParameterCache paramCache;

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

        // now comes the tricky part: calculate phase compositions
        if (phasePresence == bothPhases) {
            // both phases are present, phase compositions are a
260
            // result of the nonwetting <-> wetting equilibrium. This is
261
262
263
264
265
266
            // the job of the "MiscibleMultiPhaseComposition"
            // constraint solver
            if(useConstraintSolver) {
                MiscibleMultiPhaseComposition::solve(fluidState,
                                                     paramCache,
                                                     /*setViscosity=*/true,
267
                                                     /*setEnthalpy=*/false);
268
269
270
271
272
273
274
275
276
            }
            // ... 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;

277
278
279
                if (useKelvinEquation)
                    partPressH2O = FluidSystem::kelvinVaporPressure(fluidState, wPhaseIdx, wCompIdx);

280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
                // 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);
            }
        }
303
304
        else if (phasePresence == nPhaseOnly)
        {
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
            // only the nonwetting phase is present, i.e. nonwetting phase
            // composition is stored explicitly.
            if(useMoles)
            {
                fluidState.setMoleFraction(nPhaseIdx, nCompIdx, 1 - priVars[switchIdx]);
                fluidState.setMoleFraction(nPhaseIdx, wCompIdx, priVars[switchIdx]);
            }
            else
            {
                // setMassFraction() has only to be called 1-numComponents times
                fluidState.setMassFraction(nPhaseIdx, wCompIdx, priVars[switchIdx]);
            }

            // calculate the composition of the remaining phases (as
            // well as the densities of all phases). This is the job
            // of the "ComputeFromReferencePhase" constraint solver
            if (useConstraintSolver) {
                ComputeFromReferencePhase::solve(fluidState,
                                                 paramCache,
                                                 nPhaseIdx,
                                                 /*setViscosity=*/true,
326
                                                 /*setEnthalpy=*/false);
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
            }
            // ... 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);
            }
        }
360
361
        else if (phasePresence == wPhaseOnly)
        {
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
            // only the wetting phase is present, i.e. wetting phase
            // composition is stored explicitly.
            if(useMoles) // mole-fraction formulation
            {
                fluidState.setMoleFraction(wPhaseIdx, wCompIdx, 1-priVars[switchIdx]);
                fluidState.setMoleFraction(wPhaseIdx, nCompIdx, priVars[switchIdx]);
            }
            else // mass-fraction formulation
            {
                // setMassFraction() has only to be called 1-numComponents times
                fluidState.setMassFraction(wPhaseIdx, nCompIdx, priVars[switchIdx]);
            }

            // calculate the composition of the remaining phases (as
            // well as the densities of all phases). This is the job
            // of the "ComputeFromReferencePhase" constraint solver
            if (useConstraintSolver) {
                ComputeFromReferencePhase::solve(fluidState,
                                                 paramCache,
                                                 wPhaseIdx,
                                                 /*setViscosity=*/true,
383
                                                 /*setEnthalpy=*/false);
384
385
            }
            // ... or calculated explicitly this way ...
386
387
            else
            {
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
                // 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);
            }
        }

414
415
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
        {
Timo Koch's avatar
Timo Koch committed
416
            // set the viscosity and desity here if constraintsolver is not used
417
418
            if(!useConstraintSolver)
            {
Timo Koch's avatar
Timo Koch committed
419
420
421
                paramCache.updateComposition(fluidState, phaseIdx);
                const Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
                fluidState.setDensity(phaseIdx, rho);
422
                const Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
423
424
425
                fluidState.setViscosity(phaseIdx,mu);
            }
            // compute and set the enthalpy
426
            Scalar h = ParentType::enthalpy(fluidState, paramCache, phaseIdx);
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
            fluidState.setEnthalpy(phaseIdx, h);
        }
   }


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

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

    /*!
     * \brief Returns the mass fraction of a given component in a
     *        given phase within the control volume in \f$[-]\f$.
     *
     * \param phaseIdx The phase index
     * \param compIdx The component index
     */
    Scalar massFraction(const int phaseIdx, const int compIdx) const
    { return fluidState_.massFraction(phaseIdx, compIdx); }

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

    /*!
     * \brief Returns the mass density of a given phase within the
     *        control volume in \f$[kg/m^3]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar density(const int phaseIdx) const
    { return fluidState_.density(phaseIdx); }

    /*!
     * \brief Returns the dynamic viscosity of the fluid within the
     *        control volume in \f$\mathrm{[Pa s]}\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar viscosity(const int phaseIdx) const
    { return fluidState_.viscosity(phaseIdx); }

    /*!
     * \brief Returns the mass density of a given phase within the
     *        control volume in \f$[mol/m^3]\f$.
     *
     * \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
     *        the control volume in \f$[kg/(m*s^2)=N/m^2=Pa]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar pressure(const int phaseIdx) const
    { return fluidState_.pressure(phaseIdx); }

    /*!
     * \brief Returns temperature within the control volume in \f$[K]\f$.
     *
     * 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
     *        the control volume in \f$[-]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar relativePermeability(const int phaseIdx) const
    {
        return relativePermeability_[phaseIdx];
    }

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

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

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

548
549
550
    /*!
     * \brief Returns the average permeability within the control volume in \f$[m^2]\f$.
     */
551
    const PermeabilityType& permeability() const
552
553
    { return permeability_; }

554
555
556
    /*!
     * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
     */
557
    Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
558
559
560
561
562
563
    {
        if(phaseIdx == compIdx)
            DUNE_THROW(Dune::InvalidStateException, "Diffusion coeffiecient called for phaseIdx = compIdx");
        else
            return diffCoeff_[phaseIdx];
    }
564
565
566
567
568


protected:

    Scalar porosity_; //!< Effective porosity within the control volume
569
    PermeabilityType permeability_; //!< Effective permeability within the control volume
570
571
572
573
574
575
576
577
578
579
580
581
    Scalar relativePermeability_[numPhases]; //!< Relative permeability within the control volume
    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); }
};

582
} // end namespace Dumux
583
584

#endif