2pproblem.hh 8.96 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
// -*- 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 Tutorial problem for a fully coupled twophase box model.
 */
24
25
#ifndef DUMUX_EXERCISE_FLUIDSYSTEM_A_PROBLEM_HH
#define DUMUX_EXERCISE_FLUIDSYSTEM_A_PROBLEM_HH
26

27
// The porous media base problem
28
#include <dumux/porousmediumflow/problem.hh>
29
#include <dumux/common/properties.hh>
30
31
32
33
34

// The water component
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/components/h2o.hh>

35
36
37
// The interface to create plots during simulation
#include <dumux/io/gnuplotinterface.hh>

38
namespace Dumux {
39
40
41
42
43
44

/*!
 * \ingroup TwoPBoxModel
 * \brief  Tutorial problem for a fully coupled twophase box model.
 */
template <class TypeTag>
45
class ExerciseFluidsystemProblemTwoP : public PorousMediumFlowProblem<TypeTag>
46
47
{
    using ParentType = PorousMediumFlowProblem<TypeTag>;
48
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
49
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
50
51
52
53
54
55
56
57
58

    // 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
59
60
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
61
    using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
62
63
    using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
64
65
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using FluidState = GetPropType<TypeTag, Properties::FluidState>;
66
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
67
68
69
70
71
72
73
74
75

    enum {
        waterPressureIdx = Indices::pressureIdx,
        naplSaturationIdx = Indices::saturationIdx,
        contiWEqIdx = Indices::conti0EqIdx + FluidSystem::comp0Idx, // water transport equation index
        contiNEqIdx = Indices::conti0EqIdx + FluidSystem::comp1Idx // napl transport equation index
    };

public:
76
    ExerciseFluidsystemProblemTwoP(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
77
    : ParentType(fvGridGeometry)
78
    , eps_(3e-6)
79
80
    {

81
82
83
84
85
86
87
88
89
        // initialize the tables for the water properties
        Components::TabulatedComponent<Components::H2O<Scalar>>::init(/*tempMin=*/273.15,
                                                                      /*tempMax=*/623.15,
                                                                      /*numTemp=*/100,
                                                                      /*pMin=*/0.0,
                                                                      /*pMax=*/20e6,
                                                                      /*numP=*/200);

        // set the depth of the bottom of the reservoir
90
        depthBOR_ = this->gridGeometry().bBoxMax()[dimWorld-1];
91
92

        // plot density over pressure of the phase consisting of your component
93
94
        if(getParam<bool>("Output.PlotDensity"))
            plotDensity_();
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
    }

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

123
        if (globalPos[0] < eps_ || globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
124
125
126
127
128
129
130
131
132
133
134
135
136
           bcTypes.setAllDirichlet();
        else
            bcTypes.setAllNeumann();

        return bcTypes;
    }

    /*!
     * \brief Evaluates the boundary conditions for a Dirichlet
     *        boundary segment
     *
     * \param globalPos The global position
     */
137
    PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
138
    {
139
140
        // use the initial values as Dirichlet values
        return initialAtPos(globalPos);
141
142
143
144
145
146
147
148
149
150
151
    }

    /*!
     * \brief Evaluate the boundary conditions for a neumann
     *        boundary segment.
     *
     * \param globalPos The position of the integration point of the boundary segment.
     *
     * For this method, the \a values parameter stores the mass flux
     * in normal direction of each phase. Negative values mean influx.
     */
152
    NumEqVector neumannAtPos(const GlobalPosition& globalPos) const
153
154
    {
        // initialize values to zero, i.e. no-flow Neumann boundary conditions
155
        NumEqVector values(0.0);
156

157
        Scalar up = this->gridGeometry().bBoxMax()[dimWorld-1];
158
159
160
161

        // influx of oil (30 g/m/s) over a segment of the top boundary
        if (globalPos[dimWorld-1] > up - eps_ && globalPos[0] > 20 && globalPos[0] < 40)
        {
162
163
            values[contiWEqIdx] = 0;
            values[contiNEqIdx] = -3e-2;
164
165
166
        }
        else
        {
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
            // no-flow on the remaining Neumann-boundaries.
            values[contiWEqIdx] = 0;
            values[contiNEqIdx] = 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.
     */
190
    PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
191
192
193
194

    {
        PrimaryVariables values(0.0);

195
196
197
        // use hydrostatic pressure distribution with 2 bar at the top and zero saturation
        values[waterPressureIdx] = 200.0e3 + 9.81*1000*(depthBOR_ - globalPos[dimWorld-1]);
        values[naplSaturationIdx] = 0.0;
198
199
200
201
202
203
204
205
206
207

        return values;
    }
    // \}

    /*!
     * \brief Returns the source term
     *
     * \param globalPos The global position
     */
208
    NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
209
    {
210
        // we define no source term here
211
        NumEqVector values(0.0);
212
213
214
215
        return values;
    }

private:
216
217
    void plotDensity_()
    {
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
        FluidState fluidState;
        fluidState.setTemperature(temperature());
        int numberOfPoints = 100;
        Scalar xMin = 1e4;
        Scalar xMax = 1e7;
        Scalar spacing = std::pow((xMax/xMin), 1.0/(numberOfPoints-1));
        std::vector<double> x(numberOfPoints);
        std::vector<double> y(numberOfPoints);
        for (int i=0; i<numberOfPoints; ++i)
            x[i] = xMin*std::pow(spacing,i);
        for (int i=0; i<numberOfPoints; ++i)
        {
            fluidState.setPressure(FluidSystem::phase1Idx, x[i]);
            y[i] = FluidSystem::density(fluidState, FluidSystem::phase1Idx);
        }
233
        gnuplot_.resetPlot();
234
        gnuplot_.setXRange(xMin, xMax);
235
        gnuplot_.setOption("set logscale x 10");
236
        gnuplot_.setOption("set key left top");
237
238
239
        gnuplot_.setYRange(1440, 1480);
        gnuplot_.setXlabel("pressure [Pa]");
        gnuplot_.setYlabel("density [kg/m^3]");
240
241
        gnuplot_.addDataSetToPlot(x, y, "YourComponentPhase_density.dat", "w l t 'Phase consisting of your component'");
        gnuplot_.plot("YourComponentPhase_density");
242
243
    }

244
245
    Scalar eps_; //! small epsilon value
    Scalar depthBOR_; //! depth at the bottom of the reservoir
246
    Dumux::GnuplotInterface<double> gnuplot_; //! collects data for plotting
247
};
248
249

} // end namespace Dumux
250
251

#endif