waterairproblem.hh 11.2 KB
Newer Older
1
2
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
Bernd Flemisch's avatar
Bernd Flemisch committed
3
4
5
/*****************************************************************************
 *   Copyright (C) 2009 by Klaus Mosthaf                                     *
 *   Copyright (C) 2009 by Andreas Lauser                                    *
Andreas Lauser's avatar
Andreas Lauser committed
6
 *   Institute for Modelling Hydraulic and Environmental Systems             *
Bernd Flemisch's avatar
Bernd Flemisch committed
7
8
9
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
10
 *   This program is free software: you can redistribute it and/or modify    *
Bernd Flemisch's avatar
Bernd Flemisch committed
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.                                     *
Bernd Flemisch's avatar
Bernd Flemisch committed
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/>.   *
Bernd Flemisch's avatar
Bernd Flemisch committed
22
 *****************************************************************************/
Melanie Darcis's avatar
doc    
Melanie Darcis committed
23
24
25
/*!
 * \file
 *
26
27
 * \brief Non-isothermal gas injection problem where a gas (e.g. air)
 *        is injected into a fully water saturated medium.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
28
 */
29
30
#ifndef DUMUX_WATER_AIR_PROBLEM_HH
#define DUMUX_WATER_AIR_PROBLEM_HH
Bernd Flemisch's avatar
Bernd Flemisch committed
31
32
33
34
35

#include <dune/grid/io/file/dgfparser/dgfug.hh>
#include <dune/grid/io/file/dgfparser/dgfs.hh>
#include <dune/grid/io/file/dgfparser/dgfyasp.hh>

36
#include <dumux/material/fluidsystems/h2on2fluidsystem.hh>
Bernd Flemisch's avatar
Bernd Flemisch committed
37

38
#include <dumux/boxmodels/2p2cni/2p2cnimodel.hh>
39
#include <dumux/boxmodels/common/porousmediaboxproblem.hh>
Bernd Flemisch's avatar
Bernd Flemisch committed
40

41
#include "waterairspatialparams.hh"
Bernd Flemisch's avatar
Bernd Flemisch committed
42
43
44
45
46
47
48
49
50
51

#define ISOTHERMAL 0

namespace Dumux
{
template <class TypeTag>
class WaterAirProblem;

namespace Properties
{
52
NEW_TYPE_TAG(WaterAirProblem, INHERITS_FROM(BoxTwoPTwoCNI, WaterAirSpatialParams));
Bernd Flemisch's avatar
Bernd Flemisch committed
53
54
55
56

// Set the grid type
SET_PROP(WaterAirProblem, Grid)
{
57
    typedef Dune::YaspGrid<2> type;
Bernd Flemisch's avatar
Bernd Flemisch committed
58
59
60
61
62
};

// Set the problem property
SET_PROP(WaterAirProblem, Problem)
{
63
    typedef Dumux::WaterAirProblem<TypeTag> type;
Bernd Flemisch's avatar
Bernd Flemisch committed
64
65
};

Andreas Lauser's avatar
Andreas Lauser committed
66
// Set the wetting phase
67
SET_TYPE_PROP(WaterAirProblem, FluidSystem, Dumux::FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), false>);
Bernd Flemisch's avatar
Bernd Flemisch committed
68
69
70
71

// Enable gravity
SET_BOOL_PROP(WaterAirProblem, EnableGravity, true);

72
73
74
// Use forward differences instead of central differences
SET_INT_PROP(WaterAirProblem, NumericDifferenceMethod, +1);

Bernd Flemisch's avatar
Bernd Flemisch committed
75
// Write newton convergence
76
SET_BOOL_PROP(WaterAirProblem, NewtonWriteConvergence, false);
Bernd Flemisch's avatar
Bernd Flemisch committed
77
78
79
80
}


/*!
81
82
 * \ingroup TwoPTwoCNIModel
 * \ingroup BoxTestProblems
83
84
85
86
 * \brief Non-isothermal gas injection problem where a gas (e.g. air)
 *        is injected into a fully water saturated medium. During
 *        buoyancy driven upward migration the gas passes a high
 *        temperature area.
Bernd Flemisch's avatar
Bernd Flemisch committed
87
 *
88
89
90
 * The domain is sized 40 m times 40 m in a depth of 1000 m. The rectangular area
 * with the increased temperature (380 K) starts at (20 m, 5 m) and ends at
 * (30 m, 35 m).
Bernd Flemisch's avatar
Bernd Flemisch committed
91
92
93
94
95
96
97
98
99
100
101
102
103
104
 *
 * For the mass conservation equation neumann boundary conditions are used on
 * the top and on the bottom of the domain, while dirichlet conditions
 * apply on the left and the right boundary.
 * For the energy conservation equation dirichlet boundary conditions are applied
 * on all boundaries.
 *
 * Gas is injected at the bottom boundary from 15 m to 25 m at a rate of
 * 0.001 kg/(s m), the remaining neumann boundaries are no-flow
 * boundaries.
 *
 * At the dirichlet boundaries a hydrostatic pressure, a gas saturation of zero and
 * a geothermal temperature gradient of 0.03 K/m are applied.
 *
105
 * This problem uses the \ref TwoPTwoCNIModel.
Bernd Flemisch's avatar
Bernd Flemisch committed
106
107
108
109
110
 *
 * This problem should typically be simulated for 300000 s.
 * A good choice for the initial time step size is 1000 s.
 *
 * To run the simulation execute the following line in shell:
111
 * <tt>./test_2p2cni -parameterFile test_2p2cni.input</tt>
Bernd Flemisch's avatar
Bernd Flemisch committed
112
 *  */
Andreas Lauser's avatar
Andreas Lauser committed
113
template <class TypeTag >
114
class WaterAirProblem : public PorousMediaBoxProblem<TypeTag>
Bernd Flemisch's avatar
Bernd Flemisch committed
115
{
116
117
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
118
    typedef typename GridView::Grid Grid;
Bernd Flemisch's avatar
Bernd Flemisch committed
119

120
    typedef PorousMediaBoxProblem<TypeTag> ParentType;
Bernd Flemisch's avatar
Bernd Flemisch committed
121
122

    // copy some indices for convenience
123
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
Bernd Flemisch's avatar
Bernd Flemisch committed
124
125
126
    enum {

        pressureIdx = Indices::pressureIdx,
127
        switchIdx = Indices::switchIdx,
Bernd Flemisch's avatar
Bernd Flemisch committed
128
129
130
131
132
133
#if !ISOTHERMAL
        temperatureIdx = Indices::temperatureIdx,
        energyEqIdx = Indices::energyEqIdx,
#endif

        // Phase State
134
        wPhaseOnly = Indices::wPhaseOnly,
Bernd Flemisch's avatar
Bernd Flemisch committed
135
136

        // Grid and world dimension
137
        dim = GridView::dimension,
138
139
140
141
        dimWorld = GridView::dimensionworld,

        conti0EqIdx = Indices::conti0EqIdx,
        contiNEqIdx = conti0EqIdx + Indices::nCompIdx
Bernd Flemisch's avatar
Bernd Flemisch committed
142
143
144
    };


145
146
147
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
    typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
148
149
150
151

    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<dim>::Entity Vertex;
    typedef typename GridView::Intersection Intersection;
Bernd Flemisch's avatar
Bernd Flemisch committed
152

153
154
    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
Bernd Flemisch's avatar
Bernd Flemisch committed
155

156
    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
Bernd Flemisch's avatar
Bernd Flemisch committed
157
158

public:
Melanie Darcis's avatar
doc    
Melanie Darcis committed
159
160
161
162
163
164
    /*!
     * \brief The constructor
     *
     * \param timeManager The time manager
     * \param gridView The grid view
     */
165
166
    WaterAirProblem(TimeManager &timeManager, const GridView &gridView)
        : ParentType(timeManager, gridView)
Bernd Flemisch's avatar
Bernd Flemisch committed
167
    {
168
169
170
        maxDepth_ = 1000.0; // [m]
        eps_ = 1e-6;

Bernd Flemisch's avatar
Bernd Flemisch committed
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
        FluidSystem::init();
    }

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

    /*!
     * \brief The problem name.
     *
     * This is used as a prefix for files generated by the simulation.
     */
    const char *name() const
    { return "waterair"; }

#if ISOTHERMAL
    /*!
     * \brief Returns the temperature within the domain.
     *
Melanie Darcis's avatar
doc    
Melanie Darcis committed
191
     * \param element The element
192
     * \param fvGeometry The finite-volume geometry in the box scheme
Melanie Darcis's avatar
doc    
Melanie Darcis committed
193
194
     * \param scvIdx The local vertex index (SCV index)
     *
Bernd Flemisch's avatar
Bernd Flemisch committed
195
196
     * This problem assumes a temperature of 10 degrees Celsius.
     */
197
    Scalar temperature() const
Bernd Flemisch's avatar
Bernd Flemisch committed
198
199
200
201
202
    {
        return 273.15 + 10; // -> 10°C
    };
#endif

203
    void sourceAtPos(PrimaryVariables &values,
204
                     const GlobalPosition &globalPos) const
205
206
207
208
    {
        values = 0;
    }

Bernd Flemisch's avatar
Bernd Flemisch committed
209
210
211
212
213
214
215
216
217
218
    // \}

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

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
219
220
     *
     * \param values The boundary types for the conservation equations
221
     * \param vertex The vertex for which the boundary type is set
Bernd Flemisch's avatar
Bernd Flemisch committed
222
     */
223
    void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
Bernd Flemisch's avatar
Bernd Flemisch committed
224
    {
225
        const GlobalPosition globalPos = vertex.geometry().center();
Bernd Flemisch's avatar
Bernd Flemisch committed
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240

        if(globalPos[0] > 40 - eps_ || globalPos[0] < eps_)
            values.setAllDirichlet();
        else
            values.setAllNeumann();

#if !ISOTHERMAL
        values.setDirichlet(temperatureIdx, energyEqIdx);
#endif
    }

    /*!
     * \brief Evaluate the boundary conditions for a dirichlet
     *        boundary segment.
     *
Melanie Darcis's avatar
doc    
Melanie Darcis committed
241
     * \param values The dirichlet values for the primary variables
242
     * \param vertex The vertex for which the boundary type is set
Melanie Darcis's avatar
doc    
Melanie Darcis committed
243
     *
Bernd Flemisch's avatar
Bernd Flemisch committed
244
245
     * For this method, the \a values parameter stores primary variables.
     */
246
    void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
Bernd Flemisch's avatar
Bernd Flemisch committed
247
    {
248
        const GlobalPosition globalPos = vertex.geometry().center();
249

Bernd Flemisch's avatar
Bernd Flemisch committed
250
251
252
253
254
255
256
        initial_(values, globalPos);
    }

    /*!
     * \brief Evaluate the boundary conditions for a neumann
     *        boundary segment.
     *
Melanie Darcis's avatar
doc    
Melanie Darcis committed
257
258
     * \param values The neumann values for the conservation equations
     * \param element The finite element
259
     * \param fvGeometry The finite-volume geometry in the box scheme
Melanie Darcis's avatar
doc    
Melanie Darcis committed
260
261
262
263
     * \param is The intersection between element and boundary
     * \param scvIdx The local vertex index
     * \param boundaryFaceIdx The index of the boundary face
     *
Bernd Flemisch's avatar
Bernd Flemisch committed
264
     * For this method, the \a values parameter stores the mass flux
Melanie Darcis's avatar
doc    
Melanie Darcis committed
265
     * in normal direction of each phase. Negative values mean influx.
Bernd Flemisch's avatar
Bernd Flemisch committed
266
     */
267
268
    void neumann(PrimaryVariables &values,
                 const Element &element,
269
                 const FVElementGeometry &fvGeometry,
270
                 const Intersection &is,
271
272
                 const int scvIdx,
                 const int boundaryFaceIdx) const
Bernd Flemisch's avatar
Bernd Flemisch committed
273
274
275
276
277
278
279
280
    {
        const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
        values = 0;

        // negative values for injection
        if (globalPos[0] > 15 && globalPos[0] < 25 &&
            globalPos[1] < eps_)
        {
281
            values[contiNEqIdx] = -1e-3;
Bernd Flemisch's avatar
Bernd Flemisch committed
282
283
284
285
286
287
288
289
290
291
292
293
294
        }
    }

    // \}

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

    /*!
     * \brief Evaluate the initial value for a control volume.
     *
Melanie Darcis's avatar
doc    
Melanie Darcis committed
295
296
     * \param values The initial values for the primary variables
     * \param element The finite element
297
     * \param fvGeometry The finite-volume geometry in the box scheme
Melanie Darcis's avatar
doc    
Melanie Darcis committed
298
299
     * \param scvIdx The local vertex index
     *
Bernd Flemisch's avatar
Bernd Flemisch committed
300
301
302
     * For this method, the \a values parameter stores primary
     * variables.
     */
303
304
    void initial(PrimaryVariables &values,
                 const Element &element,
305
                 const FVElementGeometry &fvGeometry,
306
                 int scvIdx) const
Bernd Flemisch's avatar
Bernd Flemisch committed
307
    {
308
        const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
Bernd Flemisch's avatar
Bernd Flemisch committed
309
310
311
312

        initial_(values, globalPos);

#if !ISOTHERMAL
313
314
        if (globalPos[0] > 20 && globalPos[0] < 30 && globalPos[1] < 30)
            values[temperatureIdx] = 380;
Bernd Flemisch's avatar
Bernd Flemisch committed
315
316
317
318
319
#endif
    }

    /*!
     * \brief Return the initial phase state inside a control volume.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
320
321
322
323
     *
     * \param vert The vertex
     * \param globalIdx The index of the global vertex
     * \param globalPos The global position
Bernd Flemisch's avatar
Bernd Flemisch committed
324
     */
325
326
    int initialPhasePresence(const Vertex &vert,
                             int &globalIdx,
Bernd Flemisch's avatar
Bernd Flemisch committed
327
328
                             const GlobalPosition &globalPos) const
    {
329
        return wPhaseOnly;
Bernd Flemisch's avatar
Bernd Flemisch committed
330
331
332
333
334
    }

private:
    // internal method for the initial condition (reused for the
    // dirichlet conditions!)
335
    void initial_(PrimaryVariables &values,
Bernd Flemisch's avatar
Bernd Flemisch committed
336
337
338
                  const GlobalPosition &globalPos) const
    {
        Scalar densityW = 1000.0;
339
        values[pressureIdx] = 1e5 + (maxDepth_ - globalPos[1])*densityW*9.81;
Bernd Flemisch's avatar
Bernd Flemisch committed
340
341
        values[switchIdx] = 0.0;
#if !ISOTHERMAL
342
        values[temperatureIdx] = 283.0 + (maxDepth_ - globalPos[1])*0.03;
Bernd Flemisch's avatar
Bernd Flemisch committed
343
344
345
#endif
    }

346
347
    Scalar maxDepth_;
    Scalar eps_;
Bernd Flemisch's avatar
Bernd Flemisch committed
348
349
350
351
};
} //end namespace

#endif