problem.hh 12.5 KB
Newer Older
1
2
3
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
4
 *   See the file COPYING for full copying permissions.                      *
5
6
7
8
9
10
11
12
 *                                                                           *
 *   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          *
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
14
15
16
17
18
19
20
21
22
23
24
 *   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 Isothermal NAPL infiltration problem: LNAPL contaminates
 *        the unsaturated and the saturated groundwater zone.
 */
25
26
#ifndef DUMUX_NAPLINFILTRATIONPROBLEM_3P_HH
#define DUMUX_NAPLINFILTRATIONPROBLEM_3P_HH
27

28
29
#include <dumux/porousmediumflow/3p/implicit/propertydefaults.hh>
#include <dumux/porousmediumflow/implicit/problem.hh>
30

31
32
#include <lecture/mm/naplinfiltration/h2oairnaplfluidsystem.hh>
#include <lecture/mm/naplinfiltration/spatialparams.hh>
33
34
35
36
37
38
39
40

namespace Dumux
{
template <class TypeTag>
class InfiltrationProblem;

namespace Properties
{
41
NEW_TYPE_TAG(InfiltrationProblem, INHERITS_FROM(BoxThreeP, InfiltrationSpatialParams));
42
43

// Set the grid type
44
SET_TYPE_PROP(InfiltrationProblem, Grid, Dune::YaspGrid<2>);
45
46
47
48
49
50
51

// Set the problem property
SET_TYPE_PROP(InfiltrationProblem, Problem, Dumux::InfiltrationProblem<TypeTag>);

// Set the fluid system
SET_TYPE_PROP(InfiltrationProblem,
              FluidSystem,
Katharina Heck's avatar
Katharina Heck committed
52
              FluidSystems::H2OAirNAPL<TypeTag, typename GET_PROP_TYPE(TypeTag, Scalar)>);
53
54

// Enable gravity?
55
SET_BOOL_PROP(InfiltrationProblem, ProblemEnableGravity, true);
56
57
58
59
60

// Write newton convergence?
SET_BOOL_PROP(InfiltrationProblem, NewtonWriteConvergence, false);

// Maximum tolerated relative error in the Newton method
61
SET_SCALAR_PROP(InfiltrationProblem, NewtonMaxRelativeShift, 1e-4);
62
63

// -1 backward differences, 0: central differences, +1: forward differences
64
SET_INT_PROP(InfiltrationProblem, ImplicitNumericDifferenceMethod, 1);
65
66
67

// Set the linear solver
SET_TYPE_PROP(InfiltrationProblem, LinearSolver, ILU0BiCGSTABBackend<TypeTag>);
68
69
70
71
72
73
74
75
76
77
}

/*!
 * \ingroup ThreePThreeCBoxModel
 * \ingroup BoxTestProblems
 * \brief Isothermal NAPL infiltration problem: LNAPL contaminates
 *        the unsaturated and the saturated groundwater zone.
 *
 * The 2D domain of this test problem is 500 m long and 10 m deep, where
 * the lower part represents a slightly inclined groundwater table, and the
78
 * upper part is the vadose zone.
79
80
81
82
83
84
85
86
87
88
89
90
 * A LNAPL (Non-Aqueous Phase Liquid which is lighter than water) infiltrates
 * (modelled with a Neumann boundary condition) into the vadose zone. Upon
 * reaching the water table, it spreads (since lighter than water) and migrates
 * on top of the water table in the direction of the slope.
 * On its way through the vadose zone, it leaves a trace of residually trapped
 * immobile NAPL, which can in the following dissolve and evaporate slowly,
 * and eventually be transported by advection and diffusion.
 *
 * Left and right boundaries are constant head boundaries (Dirichlet),
 * Top and bottom are Neumann boundaries, all no-flow except for the small
 * infiltration zone in the upper left part.
 *
91
 * This problem uses the \ref ThreePThreeCModel or the \ref ThreePModel.
92
93
94
 *
 * This problem should typically be simulated for 30 days.
 * A good choice for the initial time step size is 60 s.
95
 * To adjust the simulation time it is necessary to edit the file test_naplfiltrationexercise.input
96
97
 *
 * To run the simulation execute the following line in shell:
98
 * <tt>./test_naplfiltrationexercise -ParameterFile test_naplfiltrationexercise.input</tt>
99
100
 *  */
template <class TypeTag >
Bernd Flemisch's avatar
Bernd Flemisch committed
101
class InfiltrationProblem : public ImplicitPorousMediaProblem<TypeTag>
102
{
Bernd Flemisch's avatar
Bernd Flemisch committed
103
    typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
104
105
106
107
108
109
110
111

    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;

    // copy some indices for convenience
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
    enum {
        pressureIdx = Indices::pressureIdx,
112
113
114
115
116
117
        switch1Idx = Indices::swIdx,
        switch2Idx = Indices::snIdx,

        contiWEqIdx = Indices::contiWEqIdx,
        contiNEqIdx = Indices::contiNEqIdx,
        contiGEqIdx = Indices::contiGEqIdx,
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133

        // Grid and world dimension
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
    };


    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
    typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;

    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<dim>::Entity Vertex;

    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
134
135
    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
    typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
136
137
138
139
140
141
142
143
144
145
146

    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;

public:
    /*!
     * \brief The constructor
     *
     * \param timeManager The time manager
     * \param gridView The grid view
     */
    InfiltrationProblem(TimeManager &timeManager, const GridView &gridView)
147
        : ParentType(timeManager, gridView), eps_(1e-5)
148
149
    {
        temperature_ = 273.15 + 10.0; // -> 10 degrees Celsius
150

151
        episodeLength_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, EpisodeLength);
152
153
        this->timeManager().startNextEpisode(episodeLength_);

154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
        FluidSystem::init(/*tempMin=*/temperature_ - 1,
                          /*tempMax=*/temperature_ + 1,
                          /*nTemp=*/3,
                          /*pressMin=*/0.8*1e5,
                          /*pressMax=*/3*1e5,
                          /*nPress=*/200);
    }

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

    /*!
     * \brief The problem name.
     *
     * This is used as a prefix for files generated by the simulation.
     */
172
173
    std::string name() const
    { return "infiltration3p"; }
174
175
176
177
178
179

    /*!
     * \brief Returns the temperature within the domain.
     *
     * This problem assumes a temperature of 10 degrees Celsius.
     */
Holger Class's avatar
Holger Class committed
180
    Scalar temperature() const
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
    {
        return temperature_;
    };

    void sourceAtPos(PrimaryVariables &values,
                     const GlobalPosition &globalPos) const
    {
        values = 0;
    }

    // \}

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

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     *
     * \param values The boundary types for the conservation equations
203
     * \param globalPos The position of the finite volume in global coordinates
204
     */
205
206
    void boundaryTypesAtPos(BoundaryTypes &values,
                            const GlobalPosition &globalPos) const
207
    {
208
        if(globalPos[0] > this->bBoxMax()[0] - eps_)
209
            values.setAllDirichlet();
210
        else if(globalPos[0] < this->bBoxMin()[0] + eps_)
211
212
213
214
215
216
217
            values.setAllDirichlet();
        else
            values.setAllNeumann();
    }

    /*!
     * \brief Evaluate the boundary conditions for a dirichlet
218
     *        control volume.
219
220
     *
     * \param values The dirichlet values for the primary variables
221
222
223
     * \param globalPos The position of the center of the finite volume
     *            for which the dirichlet condition ought to be
     *            set in global coordinates
224
225
226
     *
     * For this method, the \a values parameter stores primary variables.
     */
227
228
    void dirichletAtPos(PrimaryVariables &values,
                        const GlobalPosition &globalPos) const
229
    {
230
        initial_(values, globalPos);
231
232
233
234
235
236
    }

    /*!
     * \brief Evaluate the boundary conditions for a neumann
     *        boundary segment.
     *
237
238
239
     * \param values The neumann values for the conservation equations in units of
     *                 \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$
     * \param globalPos The position of the boundary face's integration point in global coordinates
240
241
242
243
     *
     * For this method, the \a values parameter stores the mass flux
     * in normal direction of each phase. Negative values mean influx.
     */
244
245
    void neumannAtPos(PrimaryVariables &values,
                      const GlobalPosition &globalPos) const
246
247
248
249
    {
        values = 0;

        // negative values for injection
250
        if (this->timeManager().time() < 2592000)
251
        {
252
253
254
            if (globalPos[0] < 175.0 + eps_
                && globalPos[0] > 150.0 - eps_
                && globalPos[1] > this->bBoxMax()[1] - eps_) // upper boundary
255
            {
256
                values[contiWEqIdx] = 0.0;
257
258
                // mole flux conversion via M(mesit.) = 0,120 kg/mol --> 1.2e-4  kg/(s*m)
                values[contiNEqIdx] = -0.001*0.12; // 3p needs a mass fraction
259
                values[contiGEqIdx] = 0.0;
260
            }
261
262
263
264
265
266
267
268
269
270
271
272
273
        }
    }

    // \}

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

    /*!
     * \brief Evaluate the initial value for a control volume.
     *
274
275
276
277
     * \param values The dirichlet values for the primary variables
     * \param globalPos The position of the center of the finite volume
     *            for which the initial values ought to be
     *            set (in global coordinates)
278
     *
279
     * For this method, the \a values parameter stores primary variables.
280
     */
281
282
    void initialAtPos(PrimaryVariables &values,
                      const GlobalPosition &globalPos) const
283
    {
284
        initial_(values, globalPos);
285
286
    }

287
288
    bool shouldWriteOutput() const
    {
289
    return this->timeManager().timeStepIndex() == 0 ||
290
           this->timeManager().episodeWillBeFinished() ||
291
           this->timeManager().willBeFinished();
292
293
    }

294
295
296
    void episodeEnd()
    {
        this->timeManager().startNextEpisode(episodeLength_);
297
    }
298
299
300
301
302

private:
    // internal method for the initial condition (reused for the
    // dirichlet conditions!)
    void initial_(PrimaryVariables &values,
303
                  const GlobalPosition &globalPos) const
304
    {
305
306
307
        const auto& materialLawParams = this->spatialParams().materialLawParams();
        const Scalar swr = materialLawParams.swr();
        const Scalar sgr = materialLawParams.sgr();
308

309
        if(globalPos[1] >= -1e-3*globalPos[0] + 5)
310
        {
311
312
            const Scalar pc = std::max(0.0, 9.81*1000.0*(globalPos[1] + 5e-4*globalPos[0] - 5));
            const Scalar sw = std::min(1.0-sgr, std::max(swr, invertPcGW_(pc, materialLawParams)));
313

314
            values[pressureIdx] = 1.0e5;
315
            values[switch1Idx] = sw;
316
317
318
319
320
321
322
            values[switch2Idx] = 0.0;
        }
        else
        {
            values[pressureIdx] = 1.0e5 + 9.81*1000.0*(5 - 5.0e-4*globalPos[0] - globalPos[1]);
            values[switch1Idx] = 1.0 - sgr;
            values[switch2Idx] = 0.0;
323
324
325
        }
    }

326
327
    static Scalar invertPcGW_(const Scalar pcIn,
                              const MaterialLawParams &pcParams)
328
    {
329
330
331
332
333
        Scalar lower(0.0);
        Scalar upper(1.0);
        const unsigned int maxIterations = 25;
        const Scalar bisLimit = 1.0;

334
        Scalar sw, pcGW;
335
        for (unsigned int k = 1; k <= maxIterations; k++)
336
        {
337
            sw = 0.5*(upper + lower);
338
            pcGW = MaterialLaw::pcgw(pcParams, sw);
339
340
341
342
343
344
345
346
347
348
349
            const Scalar delta = std::abs(pcGW - pcIn);
            if (delta < bisLimit)
                return sw;

            if (k == maxIterations)
                return sw;

            if (pcGW > pcIn)
                lower = sw;
            else
                upper = sw;
350
        }
351
        return sw;
352
353
354
    }

    Scalar temperature_;
355
    const Scalar eps_;
356
    Scalar episodeLength_;
357
358
359
360
};
} //end namespace

#endif