volumevariables.hh 18.8 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
23
24
25
26
27
// -*- 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
 *        finite volume in the CO2 model.
 */
#ifndef DUMUX_CO2_VOLUME_VARIABLES_HH
#define DUMUX_CO2_VOLUME_VARIABLES_HH

Bernd Flemisch's avatar
Bernd Flemisch committed
28
29
#include <array>

30
31
32
#include <dune/common/exceptions.hh>

#include <dumux/porousmediumflow/volumevariables.hh>
33
#include <dumux/porousmediumflow/2p/formulation.hh>
Bernd Flemisch's avatar
Bernd Flemisch committed
34
#include <dumux/porousmediumflow/nonisothermal/volumevariables.hh>
Timo Koch's avatar
Timo Koch committed
35
#include <dumux/material/solidstates/updatesolidvolumefractions.hh>
36

37
38
namespace Dumux {

39
/*!
40
 * \ingroup CO2Model
41
42
43
 * \brief Contains the quantities which are are constant within a
 *        finite volume in the CO2 model.
 */
44
45
template <class Traits>
class TwoPTwoCCO2VolumeVariables
46
: public PorousMediumFlowVolumeVariables<Traits>
Katharina Heck's avatar
Katharina Heck committed
47
, public EnergyVolumeVariables<Traits, TwoPTwoCCO2VolumeVariables<Traits> >
48
{
49
50
    using ParentType = PorousMediumFlowVolumeVariables< Traits>;
    using EnergyVolVars = EnergyVolumeVariables<Traits, TwoPTwoCCO2VolumeVariables<Traits> >;
51
52
53

    using Scalar = typename Traits::PrimaryVariables::value_type;
    using ModelTraits = typename Traits::ModelTraits;
Katharina Heck's avatar
Katharina Heck committed
54
    static constexpr int numFluidComps = ParentType::numComponents();
55

56
57
58
    // component indices
    enum
    {
59
60
61
62
        comp0Idx = Traits::FluidSystem::comp0Idx,
        comp1Idx = Traits::FluidSystem::comp1Idx,
        phase0Idx = Traits::FluidSystem::phase0Idx,
        phase1Idx = Traits::FluidSystem::phase1Idx
63
64
    };

65
66
67
    // phase presence indices
    enum
    {
68
69
        firstPhaseOnly = ModelTraits::Indices::firstPhaseOnly,
        secondPhaseOnly = ModelTraits::Indices::secondPhaseOnly,
70
        bothPhases = ModelTraits::Indices::bothPhases
71
72
73
    };

    // primary variable indices
74
75
76
77
    enum
    {
        switchIdx = ModelTraits::Indices::switchIdx,
        pressureIdx = ModelTraits::Indices::pressureIdx
78
79
    };

80
    // formulation
81
    static constexpr auto formulation = ModelTraits::priVarFormulation();
82

83
    // type used for the permeability
84
    using PermeabilityType = typename Traits::PermeabilityType;
85
86
public:
    //! The type of the object returned by the fluidState() method
87
88
89
    using FluidState = typename Traits::FluidState;
    //! The fluid system used here
    using FluidSystem = typename Traits::FluidSystem;
90
91
92
93
94
    //! export type of solid state
    using SolidState = typename Traits::SolidState;
    //! export type of solid system
    using SolidSystem = typename Traits::SolidSystem;

95
96
97
98
99
100
101
102
103

    //! return whether moles or masses are balanced
    static constexpr bool useMoles() { return ModelTraits::useMoles(); }
    //! return the two-phase formulation used here
    static constexpr TwoPFormulation priVarFormulation() { return formulation; }

    // check for permissive combinations
    static_assert(ModelTraits::numPhases() == 2, "NumPhases set in the model is not two!");
    static_assert(ModelTraits::numComponents() == 2, "NumComponents set in the model is not two!");
104
    static_assert((formulation == TwoPFormulation::p0s1 || formulation == TwoPFormulation::p1s0), "Chosen TwoPFormulation not supported!");
105
106
107
108
109
110
111
112
113
114
115
116
117
118

    /*!
     * \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)
    {
        ParentType::update(elemSol, problem, element, scv);
119
        completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
Timo Koch's avatar
Timo Koch committed
120

121
122
123
124
125
126
127
        // 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_);

        using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
        const auto& matParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
128
129
130
131
132
133
134
135
136
137
138

        const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
        const int nPhaseIdx = 1 - wPhaseIdx;

        // relative permeabilities -> require wetting phase saturation as parameter!
        relativePermeability_[wPhaseIdx] = MaterialLaw::krw(matParams, saturation(wPhaseIdx));
        relativePermeability_[nPhaseIdx] = MaterialLaw::krn(matParams, saturation(wPhaseIdx));

        // binary diffusion coefficients
        diffCoeff_[phase0Idx] = FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phase0Idx, comp0Idx, comp1Idx);
        diffCoeff_[phase1Idx] = FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phase1Idx, comp0Idx, comp1Idx);
139
140

        // porosity & permeabilty
Katharina Heck's avatar
Katharina Heck committed
141
        updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
142
        EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
        permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
    }

    /*!
     * \brief Complete the fluid state
     * \note TODO: This is a lot of copy paste from the 2p2c: factor out code!
     *
     * \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.
     */
    template<class ElemSol, class Problem, class Element, class Scv>
159
160
161
162
    void completeFluidState(const ElemSol& elemSol,
                            const Problem& problem,
                            const Element& element,
                            const Scv& scv,
163
164
                            FluidState& fluidState,
                            SolidState& solidState)
165
    {
166
        EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState);
Timo Koch's avatar
Timo Koch committed
167

168
        const auto& priVars = elemSol[scv.localDofIndex()];
Timo Koch's avatar
Timo Koch committed
169
170
        const auto phasePresence = priVars.state();

171
172
173
174
175
        using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
        const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
        fluidState.setWettingPhase(wPhaseIdx);

Timo Koch's avatar
Timo Koch committed
176
        // set the saturations
177
178
179
180
181
182
183
184
185
186
        if (phasePresence == secondPhaseOnly)
        {
            fluidState.setSaturation(phase0Idx, 0.0);
            fluidState.setSaturation(phase1Idx, 1.0);
        }
        else if (phasePresence == firstPhaseOnly)
        {
            fluidState.setSaturation(phase0Idx, 1.0);
            fluidState.setSaturation(phase1Idx, 0.0);
        }
187
188
        else if (phasePresence == bothPhases)
        {
189
190
191
192
193
            if (formulation == TwoPFormulation::p0s1)
            {
                fluidState.setSaturation(phase1Idx, priVars[switchIdx]);
                fluidState.setSaturation(phase0Idx, 1 - priVars[switchIdx]);
            }
194
            else
195
196
197
198
            {
                fluidState.setSaturation(phase0Idx, priVars[switchIdx]);
                fluidState.setSaturation(phase1Idx, 1 - priVars[switchIdx]);
            }
Timo Koch's avatar
Timo Koch committed
199
        }
200
201
        else
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase presence.");
Timo Koch's avatar
Timo Koch committed
202

203
204
205
206
207
208
209
        // set pressures of the fluid phases
        pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
        if (formulation == TwoPFormulation::p0s1)
        {
            fluidState.setPressure(phase0Idx, priVars[pressureIdx]);
            fluidState.setPressure(phase1Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] + pc_
                                                                       : priVars[pressureIdx] - pc_);
Timo Koch's avatar
Timo Koch committed
210
        }
211
212
213
214
215
        else
        {
            fluidState.setPressure(phase1Idx, priVars[pressureIdx]);
            fluidState.setPressure(phase0Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] - pc_
                                                                       : priVars[pressureIdx] + pc_);
Timo Koch's avatar
Timo Koch committed
216
217
218
219
220
221
222
223
224
225
        }

        // calculate the phase compositions
        typename FluidSystem::ParameterCache paramCache;
        // both phases are present
        if (phasePresence == bothPhases)
        {
            //Get the equilibrium mole fractions from the FluidSystem and set them in the fluidState
            //xCO2 = equilibrium mole fraction of CO2 in the liquid phase
            //yH2O = equilibrium mole fraction of H2O in the gas phase
226
227
            const auto xwCO2 = FluidSystem::equilibriumMoleFraction(fluidState, paramCache, phase0Idx);
            const auto xgH2O = FluidSystem::equilibriumMoleFraction(fluidState, paramCache, phase1Idx);
Timo Koch's avatar
Timo Koch committed
228
229
            const auto xwH2O = 1 - xwCO2;
            const auto xgCO2 = 1 - xgH2O;
230
231
232
233
            fluidState.setMoleFraction(phase0Idx, comp0Idx, xwH2O);
            fluidState.setMoleFraction(phase0Idx, comp1Idx, xwCO2);
            fluidState.setMoleFraction(phase1Idx, comp0Idx, xgH2O);
            fluidState.setMoleFraction(phase1Idx, comp1Idx, xgCO2);
Timo Koch's avatar
Timo Koch committed
234
235
236
237
        }

        // only the nonwetting phase is present, i.e. nonwetting phase
        // composition is stored explicitly.
238
        else if (phasePresence == secondPhaseOnly)
Timo Koch's avatar
Timo Koch committed
239
        {
240
            if( useMoles() ) // mole-fraction formulation
241
242
            {
                // set the fluid state
243
244
                fluidState.setMoleFraction(phase1Idx, comp0Idx, priVars[switchIdx]);
                fluidState.setMoleFraction(phase1Idx, comp1Idx, 1-priVars[switchIdx]);
245
                // TODO give values for non-existing wetting phase
246
                const auto xwCO2 = FluidSystem::equilibriumMoleFraction(fluidState, paramCache, phase0Idx);
Timo Koch's avatar
Timo Koch committed
247
                const auto xwH2O = 1 - xwCO2;
248
249
                fluidState.setMoleFraction(phase0Idx, comp1Idx, xwCO2);
                fluidState.setMoleFraction(phase0Idx, comp0Idx, xwH2O);
250
251
252
253
            }
            else // mass-fraction formulation
            {
                // setMassFraction() has only to be called 1-numComponents times
254
                fluidState.setMassFraction(phase1Idx, comp0Idx, priVars[switchIdx]);
255
                // TODO give values for non-existing wetting phase
256
                const auto xwCO2 = FluidSystem::equilibriumMoleFraction(fluidState, paramCache, phase0Idx);
Timo Koch's avatar
Timo Koch committed
257
                const auto xwH2O = 1 - xwCO2;
258
259
                fluidState.setMoleFraction(phase0Idx, comp1Idx, xwCO2);
                fluidState.setMoleFraction(phase0Idx, comp0Idx, xwH2O);
260
            }
Timo Koch's avatar
Timo Koch committed
261
        }
262

Timo Koch's avatar
Timo Koch committed
263
264
        // only the wetting phase is present, i.e. wetting phase
        // composition is stored explicitly.
265
        else if (phasePresence == firstPhaseOnly)
Timo Koch's avatar
Timo Koch committed
266
        {
267
            if( useMoles() ) // mole-fraction formulation
Timo Koch's avatar
Timo Koch committed
268
269
            {
                // convert mass to mole fractions and set the fluid state
270
271
                fluidState.setMoleFraction(phase0Idx, comp0Idx, 1-priVars[switchIdx]);
                fluidState.setMoleFraction(phase0Idx, comp1Idx, priVars[switchIdx]);
Timo Koch's avatar
Timo Koch committed
272
                //  TODO give values for non-existing nonwetting phase
273
                Scalar xnH2O = FluidSystem::equilibriumMoleFraction(fluidState, paramCache, phase1Idx);
Timo Koch's avatar
Timo Koch committed
274
                Scalar xnCO2 = 1 - xnH2O;
275
276
                fluidState.setMoleFraction(phase1Idx, comp1Idx, xnCO2);
                fluidState.setMoleFraction(phase1Idx, comp0Idx, xnH2O);
Timo Koch's avatar
Timo Koch committed
277
278
279
280
            }
            else // mass-fraction formulation
            {
                // setMassFraction() has only to be called 1-numComponents times
281
                fluidState.setMassFraction(phase0Idx, comp1Idx, priVars[switchIdx]);
Timo Koch's avatar
Timo Koch committed
282
                //  TODO give values for non-existing nonwetting phase
283
                Scalar xnH2O = FluidSystem::equilibriumMoleFraction(fluidState, paramCache, phase1Idx);
Timo Koch's avatar
Timo Koch committed
284
                Scalar xnCO2 = 1 - xnH2O;
285
286
                fluidState.setMoleFraction(phase1Idx, comp1Idx, xnCO2);
                fluidState.setMoleFraction(phase1Idx, comp0Idx, xnH2O);
Timo Koch's avatar
Timo Koch committed
287
288
289
            }
        }

290
        for (int phaseIdx = 0; phaseIdx < ModelTraits::numPhases(); ++phaseIdx)
Timo Koch's avatar
Timo Koch committed
291
292
293
294
295
296
297
298
299
        {
            // set the viscosity and desity here if constraintsolver is not used
            paramCache.updateComposition(fluidState, phaseIdx);
            const Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
            fluidState.setDensity(phaseIdx, rho);
            const Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
            fluidState.setViscosity(phaseIdx,mu);

            // compute and set the enthalpy
300
            Scalar h = EnergyVolVars::enthalpy(fluidState, paramCache, phaseIdx);
Timo Koch's avatar
Timo Koch committed
301
302
            fluidState.setEnthalpy(phaseIdx, h);
        }
303
    }
304
305
306
307
308
309
310

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

311
312
313
314
315
316
317
    /*!
     * \brief Returns the phase state for the control volume.
     */
    const SolidState &solidState() const
    { return solidState_; }


318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
    /*!
     * \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
416
    { return fluidState_.pressure(phase1Idx) - fluidState_.pressure(phase0Idx); }
417
418
419
420
421

    /*!
     * \brief Returns the average porosity within the control volume in \f$[-]\f$.
     */
    Scalar porosity() const
422
    { return solidState_.porosity(); }
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442

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

    /*!
     * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
     */
    Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
    {
        if(phaseIdx == compIdx)
            DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for phaseIdx = compIdx");
        else
            return diffCoeff_[phaseIdx];
    }

private:
    FluidState fluidState_;
443
    SolidState solidState_;
444
    Scalar pc_;                     //!< The capillary pressure
445
446
447
448
449
450
451
    PermeabilityType permeability_; //!< Effective permeability within the control volume

    //!< Relative permeability within the control volume
    std::array<Scalar, ModelTraits::numPhases()> relativePermeability_;

    //!< Binary diffusion coefficients for the phases
    std::array<Scalar, ModelTraits::numPhases()> diffCoeff_;
452
453
};

Timo Koch's avatar
Timo Koch committed
454
} // end namespace Dumux
455
456

#endif