problem.hh 12.4 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
68
69
70
71
72
73
74
}

/*!
 * \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
75
 * upper part is the vadose zone.
76
77
78
79
80
81
82
83
84
85
86
87
 * 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.
 *
88
 * This problem uses the \ref ThreePThreeCModel or the \ref ThreePModel.
89
90
91
 *
 * This problem should typically be simulated for 30 days.
 * A good choice for the initial time step size is 60 s.
92
 * To adjust the simulation time it is necessary to edit the file test_naplfiltrationexercise.input
93
94
 *
 * To run the simulation execute the following line in shell:
95
 * <tt>./test_naplfiltrationexercise -ParameterFile test_naplfiltrationexercise.input</tt>
96
97
 *  */
template <class TypeTag >
Bernd Flemisch's avatar
Bernd Flemisch committed
98
class InfiltrationProblem : public ImplicitPorousMediaProblem<TypeTag>
99
{
Bernd Flemisch's avatar
Bernd Flemisch committed
100
    typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
101
102
103
104
105
106
107
108

    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,
109
110
111
112
113
114
        switch1Idx = Indices::swIdx,
        switch2Idx = Indices::snIdx,

        contiWEqIdx = Indices::contiWEqIdx,
        contiNEqIdx = Indices::contiNEqIdx,
        contiGEqIdx = Indices::contiGEqIdx,
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130

        // 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;
131
132
    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
    typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
133
134
135
136
137
138
139
140
141
142
143

    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)
144
        : ParentType(timeManager, gridView), eps_(1e-5)
145
146
    {
        temperature_ = 273.15 + 10.0; // -> 10 degrees Celsius
147

148
        episodeLength_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, EpisodeLength);
149
150
        this->timeManager().startNextEpisode(episodeLength_);

151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
        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.
     */
169
170
    std::string name() const
    { return "infiltration3p"; }
171
172
173
174
175
176

    /*!
     * \brief Returns the temperature within the domain.
     *
     * This problem assumes a temperature of 10 degrees Celsius.
     */
Holger Class's avatar
Holger Class committed
177
    Scalar temperature() const
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
    {
        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
200
     * \param globalPos The position of the finite volume in global coordinates
201
     */
202
203
    void boundaryTypesAtPos(BoundaryTypes &values,
                            const GlobalPosition &globalPos) const
204
    {
205
        if(globalPos[0] > this->bBoxMax()[0] - eps_)
206
            values.setAllDirichlet();
207
        else if(globalPos[0] < this->bBoxMin()[0] + eps_)
208
209
210
211
212
213
214
            values.setAllDirichlet();
        else
            values.setAllNeumann();
    }

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

    /*!
     * \brief Evaluate the boundary conditions for a neumann
     *        boundary segment.
     *
234
235
236
     * \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
237
238
239
240
     *
     * For this method, the \a values parameter stores the mass flux
     * in normal direction of each phase. Negative values mean influx.
     */
241
242
    void neumannAtPos(PrimaryVariables &values,
                      const GlobalPosition &globalPos) const
243
244
245
246
    {
        values = 0;

        // negative values for injection
247
        if (this->timeManager().time() < 2592000)
248
        {
249
250
251
            if (globalPos[0] < 175.0 + eps_
                && globalPos[0] > 150.0 - eps_
                && globalPos[1] > this->bBoxMax()[1] - eps_) // upper boundary
252
            {
253
                values[contiWEqIdx] = 0.0;
254
255
                // 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
256
                values[contiGEqIdx] = 0.0;
257
            }
258
259
260
261
262
263
264
265
266
267
268
269
270
        }
    }

    // \}

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

    /*!
     * \brief Evaluate the initial value for a control volume.
     *
271
272
273
274
     * \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)
275
     *
276
     * For this method, the \a values parameter stores primary variables.
277
     */
278
279
    void initialAtPos(PrimaryVariables &values,
                      const GlobalPosition &globalPos) const
280
    {
281
        initial_(values, globalPos);
282
283
    }

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

291
292
293
    void episodeEnd()
    {
        this->timeManager().startNextEpisode(episodeLength_);
294
    }
295
296
297
298
299

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

306
        if(globalPos[1] >= -1e-3*globalPos[0] + 5)
307
        {
308
309
            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)));
310

311
            values[pressureIdx] = 1.0e5;
312
            values[switch1Idx] = sw;
313
314
315
316
317
318
319
            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;
320
321
322
        }
    }

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

331
        Scalar sw, pcGW;
332
        for (unsigned int k = 1; k <= maxIterations; k++)
333
        {
334
            sw = 0.5*(upper + lower);
335
            pcGW = MaterialLaw::pcgw(pcParams, sw);
336
337
338
339
340
341
342
343
344
345
346
            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;
347
        }
348
        return sw;
349
350
351
    }

    Scalar temperature_;
352
    const Scalar eps_;
353
    Scalar episodeLength_;
354
355
356
357
};
} //end namespace

#endif