2p2cvolumevariables.hh 9.82 KB
Newer Older
1
// $Id$
2
3
4
5
6
7
8
9
/*****************************************************************************
 *   Copyright (C) 2008,2009 by Klaus Mosthaf,                               *
 *                              Andreas Lauser,                              *
 *                              Bernd Flemisch                               *
 *   Institute of Hydraulic Engineering                                      *
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
10
 *   This program is free software: you can redistribute it and/or modify    *
11
 *   it under the terms of the GNU General Public License as published by    *
12
13
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
14
 *                                                                           *
15
16
17
18
19
20
21
 *   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/>.   *
22
23
24
25
 *****************************************************************************/
/*!
 * \file
 *
Felix Bode's avatar
Felix Bode committed
26
 * \brief Contains the quantities which are constant within a
27
28
 *        finite volume in the two-phase, two-component model.
 */
29
30
#ifndef DUMUX_2P2C_VOLUME_VARIABLES_HH
#define DUMUX_2P2C_VOLUME_VARIABLES_HH
31

32
#include <dumux/boxmodels/common/boxmodel.hh>
33
34
35
36
37
38
#include <dumux/common/math.hh>

#include <dune/common/collectivecommunication.hh>
#include <vector>
#include <iostream>

39
#include "2p2cproperties.hh"
40
#include "2p2cfluidstate.hh"
41
42
43
44
45

namespace Dumux
{

/*!
Melanie Darcis's avatar
doc    
Melanie Darcis committed
46
 * \ingroup TwoPTwoCModel
47
48
49
50
 * \brief Contains the quantities which are are constant within a
 *        finite volume in the two-phase, two-component model.
 */
template <class TypeTag>
51
class TwoPTwoCVolumeVariables : public BoxVolumeVariables<TypeTag>
52
{
53
54
    typedef BoxVolumeVariables<TypeTag> ParentType;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) Implementation;
55

56
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
57
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
58
59
60
61
62
63
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLawParams)) MaterialLawParams;
64

65
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
66
    enum {
67
        dim = GridView::dimension,
68

69
70
        numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
        formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation)),
71

72
73
        lCompIdx = Indices::lCompIdx,
        gCompIdx = Indices::gCompIdx,
74

75
        lPhaseIdx = Indices::lPhaseIdx,
76
        gPhaseIdx = Indices::gPhaseIdx
77
78
    };

79
    typedef typename GridView::template Codim<0>::Entity Element;
80
    typedef TwoPTwoCFluidState<TypeTag> FluidState;
81

82
83
84
public:
    /*!
     * \brief Update all quantities for a given control volume.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
85
86
87
88
     *
     * \param priVars The primary variables
     * \param problem The problem
     * \param element The element
Karin Erbertseder's avatar
Karin Erbertseder committed
89
     * \param elemGeom The finite-volume geometry in the box scheme
Melanie Darcis's avatar
doc    
Melanie Darcis committed
90
91
     * \param scvIdx The local index of the SCV (sub-control volume)
     * \param isOldSol Evaluate function with solution of current or previous time step
92
     */
93
94
95
    void update(const PrimaryVariables &priVars,
                const Problem &problem,
                const Element &element,
96
                const FVElementGeometry &elemGeom,
97
98
                int scvIdx,
                bool isOldSol)
99
    {
100
101
102
103
104
105
        ParentType::update(priVars,
                           problem,
                           element,
                           elemGeom,
                           scvIdx,
                           isOldSol);
106

107
108
109
110
111
        asImp().updateTemperature_(priVars,
                                   element,
                                   elemGeom,
                                   scvIdx,
                                   problem);
112
113
114
115

        // capillary pressure parameters
        const MaterialLawParams &materialParams =
            problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
116
117
118
119
120
121
122
123
124

        int globalVertIdx = problem.model().dofMapper().map(element, scvIdx, dim);
        int phasePresence = problem.model().phasePresence(globalVertIdx, isOldSol);

        // calculate phase state
        fluidState_.update(priVars, materialParams, temperature(), phasePresence);
        Valgrind::CheckDefined(fluidState_);

        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
            // Mobilities
            const Scalar mu =
                FluidSystem::phaseViscosity(phaseIdx,
                                            fluidState().temperature(),
                                            fluidState().phasePressure(lPhaseIdx),
                                            fluidState());
            Scalar kr;
            if (phaseIdx == lPhaseIdx)
                kr = MaterialLaw::krw(materialParams, saturation(lPhaseIdx));
            else // ATTENTION: krn requires the liquid saturation
                // as parameter!
                kr = MaterialLaw::krn(materialParams, saturation(lPhaseIdx));
            mobility_[phaseIdx] = kr / mu;
            Valgrind::CheckDefined(mobility_[phaseIdx]);
            
            // binary diffusion coefficents
            diffCoeff_[phaseIdx] =
                FluidSystem::diffCoeff(phaseIdx,
                                       lCompIdx,
                                       gCompIdx,
                                       fluidState_.temperature(),
                                       fluidState_.phasePressure(phaseIdx),
                                       fluidState_);
            Valgrind::CheckDefined(diffCoeff_[phaseIdx]);
149
        }
150

151
152
153
154
155
        // porosity
        porosity_ = problem.spatialParameters().porosity(element,
                                                         elemGeom,
                                                         scvIdx);
        Valgrind::CheckDefined(porosity_);
156
   }
157
158
159
160
161
162
163
164
165
166

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

    /*!
     * \brief Returns the effective saturation of a given phase within
     *        the control volume.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
167
168
     *
     * \param phaseIdx The phase index
169
170
171
172
173
174
175
     */
    Scalar saturation(int phaseIdx) const
    { return fluidState_.saturation(phaseIdx); }

    /*!
     * \brief Returns the mass density of a given phase within the
     *        control volume.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
176
177
     *
     * \param phaseIdx The phase index
178
179
180
181
182
183
184
     */
    Scalar density(int phaseIdx) const
    { return fluidState_.density(phaseIdx); }

    /*!
     * \brief Returns the mass density of a given phase within the
     *        control volume.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
185
186
     *
     * \param phaseIdx The phase index
187
188
     */
    Scalar molarDensity(int phaseIdx) const
189
    { return fluidState_.density(phaseIdx) / fluidState_.averageMolarMass(phaseIdx); }
190
191
192
193

    /*!
     * \brief Returns the effective pressure of a given phase within
     *        the control volume.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
194
195
     *
     * \param phaseIdx The phase index
196
197
     */
    Scalar pressure(int phaseIdx) const
198
    { return fluidState_.phasePressure(phaseIdx); }
199
200
201
202
203
204
205
206
207

    /*!
     * \brief Returns temperature inside the sub-control volume.
     *
     * Note that we assume thermodynamic equilibrium, i.e. the
     * temperature of the rock matrix and of all fluid phases are
     * identical.
     */
    Scalar temperature() const
208
    { return temperature_; }
209
210
211
212

    /*!
     * \brief Returns the effective mobility of a given phase within
     *        the control volume.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
213
214
     *
     * \param phaseIdx The phase index
215
216
217
218
219
220
221
222
223
224
     */
    Scalar mobility(int phaseIdx) const
    {
        return mobility_[phaseIdx];
    }

    /*!
     * \brief Returns the effective capillary pressure within the control volume.
     */
    Scalar capillaryPressure() const
225
    { return fluidState_.capillaryPressure(); }
226
227
228
229
230
231
232
233
234
235
236
237

    /*!
     * \brief Returns the average porosity within the control volume.
     */
    Scalar porosity() const
    { return porosity_; }

    /*!
     * \brief Returns the binary diffusion coefficients for a phase
     */
    Scalar diffCoeff(int phaseIdx) const
    { return diffCoeff_[phaseIdx]; }
238

239
240

protected:
241

242
243
244
245
246
    void updateTemperature_(const PrimaryVariables &priVars,
                            const Element &element,
                            const FVElementGeometry &elemGeom,
                            int scvIdx,
                            const Problem &problem)
Melanie Darcis's avatar
doc    
Melanie Darcis committed
247
    {
248
        temperature_ = problem.temperature(element, elemGeom, scvIdx);
Melanie Darcis's avatar
doc    
Melanie Darcis committed
249
250
    }

251
    Scalar temperature_;     //!< Temperature within the control volume
252
253
    Scalar porosity_;        //!< Effective porosity within the control volume
    Scalar mobility_[numPhases];  //!< Effective mobility within the control volume
254
    Scalar diffCoeff_[numPhases]; //!< Binary diffusion coefficients for the phases
255
256
257
    FluidState fluidState_;

private:
258
    Implementation &asImp()
259
260
    { return *static_cast<Implementation*>(this); }

261
    const Implementation &asImp() const
262
263
264
265
266
267
    { return *static_cast<const Implementation*>(this); }
};

} // end namepace

#endif