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

26
27
28
#ifndef DUMUX_2P_VOLUME_VARIABLES_HH
#define DUMUX_2P_VOLUME_VARIABLES_HH

29
#include <dumux/porousmediumflow/volumevariables.hh>
30
#include <dumux/porousmediumflow/nonisothermal/volumevariables.hh>
Timo Koch's avatar
Timo Koch committed
31
32
#include <dumux/material/solidstates/updatesolidvolumefractions.hh>
#include <dumux/porousmediumflow/2p/formulation.hh>
33

34
35
namespace Dumux {

36
37
38
39
40
/*!
 * \ingroup TwoPModel
 * \brief Contains the quantities which are are constant within a
 *        finite volume in the two-phase model.
 */
41
42
template <class Traits>
class TwoPVolumeVariables
43
: public PorousMediumFlowVolumeVariables<Traits>
Katharina Heck's avatar
Katharina Heck committed
44
, public EnergyVolumeVariables<Traits, TwoPVolumeVariables<Traits> >
45
{
46
47
    using ParentType = PorousMediumFlowVolumeVariables<Traits>;
    using EnergyVolVars = EnergyVolumeVariables<Traits, TwoPVolumeVariables<Traits> >;
48
49
50
51
    using PermeabilityType = typename Traits::PermeabilityType;
    using ModelTraits = typename Traits::ModelTraits;
    using Indices = typename ModelTraits::Indices;
    using Scalar = typename Traits::PrimaryVariables::value_type;
52
    using FS = typename Traits::FluidSystem;
53
    static constexpr int numFluidComps = ParentType::numFluidComponents();
54
55
    enum
    {
56
57
        pressureIdx = Indices::pressureIdx,
        saturationIdx = Indices::saturationIdx,
58
59
60

        phase0Idx = FS::phase0Idx,
        phase1Idx = FS::phase1Idx
61
62
    };

63
    static constexpr auto formulation = ModelTraits::priVarFormulation();
64

65
public:
66
    //! Export type of fluid system
67
    using FluidSystem = typename Traits::FluidSystem;
68
    //! Export type of fluid state
69
    using FluidState = typename Traits::FluidState;
70
    //! Export type of solid state
71
    using SolidState = typename Traits::SolidState;
72
    //! Export type of solid system
73
    using SolidSystem = typename Traits::SolidSystem;
74
75

    /*!
76
     * \brief Updates all quantities for a given control volume.
77
78
79
80
81
82
83
     *
     * \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
    */
84
85
    template<class ElemSol, class Problem, class Element, class Scv>
    void update(const ElemSol &elemSol,
86
87
                const Problem &problem,
                const Element &element,
88
                const Scv& scv)
89
    {
90
        ParentType::update(elemSol, problem, element, scv);
91

92
        completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
93

94
        using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
95
        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
96

97
98
        const int wPhaseIdx = fluidState_.wettingPhase();
        const int nPhaseIdx = 1 - wPhaseIdx;
99

100
101
102
        mobility_[wPhaseIdx] =
            MaterialLaw::krw(materialParams, fluidState_.saturation(wPhaseIdx))
            / fluidState_.viscosity(wPhaseIdx);
103
104

        mobility_[nPhaseIdx] =
105
            MaterialLaw::krn(materialParams, fluidState_.saturation(wPhaseIdx))
106
107
            / fluidState_.viscosity(nPhaseIdx);

108
        // porosity calculation over inert volumefraction
Katharina Heck's avatar
Katharina Heck committed
109
        updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
110
        EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
111
        permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
112
        EnergyVolVars::updateEffectiveThermalConductivity();
Katharina Heck's avatar
Katharina Heck committed
113
    }
114
115

    /*!
116
     * \brief Sets complete fluid state.
117
118
     *
     * \param elemSol A vector containing all primary variables connected to the element
119
120
121
122
123
     * \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
     * \param fluidState A container with the current (physical) state of the fluid
     * \param solidState A container with the current (physical) state of the solid
124
125
     *
     * Set temperature, saturations, capillary pressures, viscosities, densities and enthalpies.
126
     */
127
    template<class ElemSol, class Problem, class Element, class Scv>
128
129
130
131
    void completeFluidState(const ElemSol& elemSol,
                            const Problem& problem,
                            const Element& element,
                            const Scv& scv,
132
133
                            FluidState& fluidState,
                            SolidState& solidState)
134
    {
135
        EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState);
136

137
        using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
138
        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
139
        const auto& priVars = elemSol[scv.localDofIndex()];
140

141
142
        const auto wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
        fluidState.setWettingPhase(wPhaseIdx);
143
144
145
        if (formulation == TwoPFormulation::p0s1)
        {
            fluidState.setPressure(phase0Idx, priVars[pressureIdx]);
146
            if (fluidState.wettingPhase() == phase1Idx)
147
148
149
            {
                fluidState.setSaturation(phase1Idx, priVars[saturationIdx]);
                fluidState.setSaturation(phase0Idx, 1 - priVars[saturationIdx]);
150
                pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
151
152
153
154
155
156
157
158
                fluidState.setPressure(phase1Idx, priVars[pressureIdx] - pc_);
            }
            else
            {
                const auto Sn = Traits::SaturationReconstruction::reconstructSn(problem.spatialParams(), element,
                                                                                scv, elemSol, priVars[saturationIdx]);
                fluidState.setSaturation(phase1Idx, Sn);
                fluidState.setSaturation(phase0Idx, 1 - Sn);
159
                pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
160
161
                fluidState.setPressure(phase1Idx, priVars[pressureIdx] + pc_);
            }
162
        }
163
164
165
        else if (formulation == TwoPFormulation::p1s0)
        {
            fluidState.setPressure(phase1Idx, priVars[pressureIdx]);
166
            if (wPhaseIdx == phase1Idx)
167
168
169
170
171
            {
                const auto Sn = Traits::SaturationReconstruction::reconstructSn(problem.spatialParams(), element,
                                                                                scv, elemSol, priVars[saturationIdx]);
                fluidState.setSaturation(phase0Idx, Sn);
                fluidState.setSaturation(phase1Idx, 1 - Sn);
172
                pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
173
174
175
176
                fluidState.setPressure(phase0Idx, priVars[pressureIdx] + pc_);
            }
            else
            {
177
178
                fluidState.setSaturation(phase0Idx, priVars[saturationIdx]);
                fluidState.setSaturation(phase1Idx, 1 - priVars[saturationIdx]);
179
                pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
180
181
                fluidState.setPressure(phase0Idx, priVars[pressureIdx] - pc_);
            }
182
183
184
185
186
        }

        typename FluidSystem::ParameterCache paramCache;
        paramCache.updateAll(fluidState);

187
        for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx) {
188
189
190
191
192
193
194
195
196
            // compute and set the viscosity
            Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
            fluidState.setViscosity(phaseIdx, mu);

            // compute and set the density
            Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
            fluidState.setDensity(phaseIdx, rho);

            // compute and set the enthalpy
197
            Scalar h = EnergyVolVars::enthalpy(fluidState, paramCache, phaseIdx);
198
199
200
201
202
203
204
205
206
207
            fluidState.setEnthalpy(phaseIdx, h);
        }
    }

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

208
209
210
211
212
213
    /*!
     * \brief Returns the phase state for the control volume.
     */
    const SolidState &solidState() const
    { return solidState_; }

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

    /*!
     * \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(int phaseIdx) const
    { return fluidState_.density(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(int phaseIdx) const
    { return fluidState_.pressure(phaseIdx); }

    /*!
     * \brief Returns the capillary pressure within the control volume
     * in \f$[kg/(m*s^2)=N/m^2=Pa]\f$.
     */
    Scalar capillaryPressure() const
246
    { return pc_; }
247
248
249
250
251
252
253
254
255
256
257
258

    /*!
     * \brief Returns temperature inside the sub-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); }

259
260
261
262
263
264
265
266
267
    /*!
     * \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(int phaseIdx) const
    { return fluidState_.viscosity(phaseIdx); }

268
269
270
271
272
273
274
275
276
277
278
279
280
    /*!
     * \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(int phaseIdx) const
    { return mobility_[phaseIdx]; }

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

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

289
290
291
    /*!
     * \brief Returns the wetting phase index
     */
292
    int wettingPhase() const
293
    {  return fluidState_.wettingPhase(); }
294

295
296
protected:
    FluidState fluidState_;
297
    SolidState solidState_;
298
299

private:
300
    Scalar pc_;
301
    Scalar porosity_;
302
    PermeabilityType permeability_;
303
    Scalar mobility_[ModelTraits::numFluidPhases()];
304
305
};

306
} // end namespace Dumux
307
308

#endif