volumevariables.hh 12.2 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
    using PermeabilityType = typename Traits::PermeabilityType;
    using ModelTraits = typename Traits::ModelTraits;
50
    using Idx = typename ModelTraits::Indices;
51
    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 = Idx::pressureIdx,
        saturationIdx = Idx::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
71
    //! Export the indices
    using Indices = typename ModelTraits::Indices;
72
    //! Export type of solid state
73
    using SolidState = typename Traits::SolidState;
74
    //! Export type of solid system
75
    using SolidSystem = typename Traits::SolidSystem;
76
77

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

94
        completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
95

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

99
100
        const int wPhaseIdx = fluidState_.wettingPhase();
        const int nPhaseIdx = 1 - wPhaseIdx;
101

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

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

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

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

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

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

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

189
        for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx) {
190
191
192
193
194
195
196
197
198
            // 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
199
            Scalar h = EnergyVolVars::enthalpy(fluidState, paramCache, phaseIdx);
200
201
202
203
204
205
206
207
208
209
            fluidState.setEnthalpy(phaseIdx, h);
        }
    }

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

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

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
    /*!
     * \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
248
    { return pc_; }
249
250
251
252
253
254
255
256
257
258
259
260

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

261
262
263
264
265
266
267
268
269
    /*!
     * \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); }

270
271
272
273
274
275
276
277
278
279
280
281
282
    /*!
     * \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
283
    { return solidState_.porosity(); }
284

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

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

297
298
protected:
    FluidState fluidState_;
299
    SolidState solidState_;
300
301

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

308
} // end namespace Dumux
309
310

#endif