2p2cproblem.hh 8.49 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// -*- 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
 *
22
 * \brief Tutorial problem for a fully coupled two phase-two component box model.
23
 */
24
25
#ifndef DUMUX_EXERCISE_FLUIDSYSTEM_B_PROBLEM_HH
#define DUMUX_EXERCISE_FLUIDSYSTEM_B_PROBLEM_HH
26

27
28
29
// The grid manager
#include <dune/grid/yaspgrid.hh>

30
31
32
// The numerical model
#include <dumux/porousmediumflow/2p2c/model.hh>

33
// The box discretization
34
#include <dumux/discretization/box.hh>
35
36
37
38
39
40
41
42
43
44
45
46

// The base porous media box problem
#include <dumux/porousmediumflow/problem.hh>

// Spatially dependent parameters
#include "spatialparams.hh"

// The fluid system that is created in this exercise
#include "fluidsystems/h2omycompressiblecomponent.hh"

namespace Dumux{
// Forward declaration of the problem class
47
template <class TypeTag> class ExerciseFluidsystemProblemTwoPTwoC;
48
49
50

namespace Properties {
// Create a new type tag for the problem
51
52
// Create new type tags
namespace TTag {
53
struct ExerciseFluidsystemTwoPTwoC { using InheritsFrom = std::tuple<TwoPTwoC, BoxModel>; };
54
} // end namespace TTag
55
56

// Set the "Problem" property
57
template<class TypeTag>
Felix Weinhardt's avatar
Felix Weinhardt committed
58
struct Problem<TypeTag, TTag::ExerciseFluidsystemTwoPTwoC> { using type = ExerciseFluidsystemProblemTwoPTwoC<TypeTag>; };
59
60

// Set the spatial parameters
61
62
63
64
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::ExerciseFluidsystemTwoPTwoC>
{
private:
65
    using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
66
67
68
69
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
    using type = ExerciseFluidsystemSpatialParams<FVGridGeometry, Scalar>;
};
70
71

// Set grid and the grid creator to be used
72
template<class TypeTag>
Felix Weinhardt's avatar
Felix Weinhardt committed
73
struct Grid<TypeTag, TTag::ExerciseFluidsystemTwoPTwoC> { using type = Dune::YaspGrid<2>; };
74
75

 // The fluid system property
76
template<class TypeTag>
Felix Weinhardt's avatar
Felix Weinhardt committed
77
struct FluidSystem<TypeTag, TTag::ExerciseFluidsystemTwoPTwoC>
78
79
{
private:
80
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
81
82
83
84
85
86
87
88
89
public:
    using type = FluidSystems::H2OMyCompressibleComponent<Scalar>;
};

}

/*!
 * \ingroup TwoPBoxModel
 *
90
 * \brief  Tutorial problem for a fully coupled two phase-two component box model.
91
92
 */
template <class TypeTag>
93
class ExerciseFluidsystemProblemTwoPTwoC : public PorousMediumFlowProblem<TypeTag>
94
95
{
    using ParentType = PorousMediumFlowProblem<TypeTag>;
96
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
97
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
98
99
100
101
102
103
104
105
106

    // Grid dimension
    enum { dim = GridView::dimension,
           dimWorld = GridView::dimensionworld
    };
    using Element = typename GridView::template Codim<0>::Entity;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;

    // Dumux specific types
107
108
109
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
110
111
    using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
112
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
113
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
114
115

public:
116
    ExerciseFluidsystemProblemTwoPTwoC(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
117
118
119
120
121
122
123
124
    : ParentType(fvGridGeometry)
        , eps_(3e-6)
    {

        // initialize the fluid system
        FluidSystem::init();

        // set the depth of the bottom of the reservoir
125
        depthBOR_ = this->gridGeometry().bBoxMax()[dimWorld-1];
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
    }

    /*!
     * \name Problem parameters
     */
    // \{

    /*!
     * \brief Returns the temperature \f$ K \f$
     */
    Scalar temperature() const
    { return 283.15; }

     /*!
     * \name Boundary conditions
     */
    // \{

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     *
     * \param globalPos The position for which the bc type should be evaluated
     */
    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
    {
         BoundaryTypes bcTypes;

154
        if (globalPos[0] < eps_ || globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
           bcTypes.setAllDirichlet();
        else
            bcTypes.setAllNeumann();

        return bcTypes;
    }

    /*!
     * \brief Evaluates the boundary conditions for a Dirichlet
     *        boundary segment
     *
     * \param globalPos The global position
     */
    PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
    {
170
171
        // use initial as Dirichlet conditions
        return initialAtPos(globalPos);
172
173
174
175
176
177
178
179
    }

    /*!
     * \brief Evaluate the boundary conditions for a neumann
     *        boundary segment.
     *
     * \param globalPos The position of the integration point of the boundary segment.
     */
180
    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
181
182
    {
        // initialize values to zero, i.e. no-flow Neumann boundary conditions
183
        NumEqVector values(0.0);
184

185
        Scalar up = this->gridGeometry().bBoxMax()[dimWorld-1];
186
187
188
        // extraction of oil (30 g/m/s) on a segment of the upper boundary
        if (globalPos[dimWorld-1] > up - eps_ && globalPos[0] > 20 && globalPos[0] < 40)
        {
189
190
191
            // we solve for the mole balance, so we have to divide by the molar mass
            values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0;
            values[Indices::conti0EqIdx + FluidSystem::NAPLIdx] = -3e-2/FluidSystem::MyCompressibleComponent::molarMass();
192
193
194
        }
        else
        {
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
            // no-flow on the remaining Neumann-boundaries.
            values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0;
            values[Indices::conti0EqIdx + FluidSystem::NAPLIdx] = 0;
        }

        return values;
    }

    // \}

    /*!
     * \name Volume terms
     */
    // \{

    /*!
     * \brief Evaluate the initial value for a control volume.
     * \param globalPos The position for which the initial condition should be evaluated
     *
     * For this method, the \a values parameter stores primary
     * variables.
     */
    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const

    {
        PrimaryVariables values(0.0);
221
222
223
224

        // tell the primary variables the phase state, i.e. which phase/phases
        // is/are present, because this changes the meaning of the primary variable
        // value at the index Indices::switchIdx
225
226
        values.setState(Indices::firstPhaseOnly);

227
        // use hydrostatic pressure distribution with 2 bar at the top and zero saturation
228
229
230
231
232
233
234
235
236
237
238
239
        values[Indices::pressureIdx] = 200.0e3 + 9.81*1000*(depthBOR_ - globalPos[dimWorld-1]); // 200 kPa = 2 bar
        values[Indices::switchIdx] = 0.0;

        return values;
    }

    // \}

    /*!
     * \brief Returns the source term
     * \param globalPos The global position
     */
240
    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
241
    {
242
        // we do not define any sources
243
        NumEqVector values(0.0);
244
245
246
247
248
        return values;
    }


private:
249
250
    Scalar eps_; //! small epsilon value
    Scalar depthBOR_; //! depth at the bottom of the reservoir
251
};
252
253

} // end namespace Dumux
254
255

#endif