heterogeneousproblem.hh 17.7 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
 *   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 Definition of a problem, where CO2 is injected in a reservoir.
23
24
25
26
 */
#ifndef DUMUX_HETEROGENEOUS_PROBLEM_HH
#define DUMUX_HETEROGENEOUS_PROBLEM_HH

Timo Koch's avatar
Timo Koch committed
27
28
29
30
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
#include <dumux/discretization/box/properties.hh>

#include <dumux/porousmediumflow/problem.hh>
Timo Koch's avatar
Timo Koch committed
31
#include <dumux/porousmediumflow/co2/model.hh>
32

Timo Koch's avatar
Timo Koch committed
33
34
#include <dumux/material/fluidsystems/brineco2.hh>
#include <dumux/discretization/box/scvftoscvboundarytypes.hh>
35

36
37
38
#include "heterogeneousspatialparameters.hh"
#include "heterogeneousco2tables.hh"

Timo Koch's avatar
Timo Koch committed
39
40
41
42
43
// per default use isothermal model
#ifndef ISOTHERMAL
#define ISOTHERMAL 1
#endif

44
45
46
47
48
49
50
51
namespace Dumux
{

template <class TypeTag>
class HeterogeneousProblem;

namespace Properties
{
Timo Koch's avatar
Timo Koch committed
52
53
54
NEW_TYPE_TAG(HeterogeneousTypeTag, INHERITS_FROM(TwoPTwoCCO2, HeterogeneousSpatialParams));
NEW_TYPE_TAG(HeterogeneousBoxTypeTag, INHERITS_FROM(BoxModel, HeterogeneousTypeTag));
NEW_TYPE_TAG(HeterogeneousCCTpfaTypeTag, INHERITS_FROM(CCTpfaModel, HeterogeneousTypeTag));
55
56

// Set the grid type
Timo Koch's avatar
Timo Koch committed
57
SET_TYPE_PROP(HeterogeneousTypeTag, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>);
58
59

// Set the problem property
Timo Koch's avatar
Timo Koch committed
60
SET_TYPE_PROP(HeterogeneousTypeTag, Problem, HeterogeneousProblem<TypeTag>);
61
62

// Set fluid configuration
Timo Koch's avatar
Timo Koch committed
63
64
65
66
67
SET_TYPE_PROP(HeterogeneousTypeTag, FluidSystem, FluidSystems::BrineCO2<typename GET_PROP_TYPE(TypeTag, Scalar),
                                                                        HeterogeneousCO2Tables::CO2Tables>);

// Use Moles
SET_BOOL_PROP(HeterogeneousTypeTag, UseMoles, false);
68

Timo Koch's avatar
Timo Koch committed
69
70
71
72
73
74
75
#if !ISOTHERMAL
NEW_TYPE_TAG(HeterogeneousNITypeTag, INHERITS_FROM(TwoPTwoCCO2NI, HeterogeneousSpatialParams));
NEW_TYPE_TAG(HeterogeneousNIBoxTypeTag, INHERITS_FROM(BoxModel, HeterogeneousNITypeTag));
NEW_TYPE_TAG(HeterogeneousNICCTpfaTypeTag, INHERITS_FROM(CCTpfaModel, HeterogeneousNITypeTag));

// Set the grid type
SET_TYPE_PROP(HeterogeneousNITypeTag, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>);
76

Timo Koch's avatar
Timo Koch committed
77
78
// Set the problem property
SET_TYPE_PROP(HeterogeneousNITypeTag, Problem, HeterogeneousProblem<TypeTag>);
79

Timo Koch's avatar
Timo Koch committed
80
81
82
// Set fluid configuration
SET_TYPE_PROP(HeterogeneousNITypeTag, FluidSystem, FluidSystems::BrineCO2<typename GET_PROP_TYPE(TypeTag, Scalar),
                                                                        HeterogeneousCO2Tables::CO2Tables>);
83

84
// Use Moles
Timo Koch's avatar
Timo Koch committed
85
86
SET_BOOL_PROP(HeterogeneousNITypeTag, UseMoles, false);
#endif
87
88
89
90
91
}


/*!
 * \ingroup CO2Model
92
 * \ingroup ImplicitTestProblems
93
 * \brief Definition of a problem, where CO2 is injected in a reservoir.
94
95
96
97
98
99
100
101
102
103
104
105
 *
 * The domain is sized 200m times 100m and consists of four layers, a
 * permeable reservoir layer at the bottom, a barrier rock layer with reduced permeability, another reservoir layer
 * and at the top a barrier rock layer with a very low permeablility.
 *
 * CO2 is injected at the permeable bottom layer
 * from the left side. The domain is initially filled with brine.
 *
 * The grid is unstructered and permeability and porosity for the elements are read in from the grid file. The grid file
 * also contains so-called boundary ids which can be used assigned during the grid creation in order to differentiate
 * between different parts of the boundary.
 * These boundary ids can be imported into the problem where the boundary conditions can then be assigned accordingly.
106
 *
107
108
109
 * The model is able to use either mole or mass fractions. The property useMoles can be set to either true or false in the
 * problem file. Make sure that the according units are used in the problem setup. The default setting for useMoles is false.
 *
110
 * To run the simulation execute the following line in shell (works with the box and cell centered spatial discretization method):
111
 * <tt>./test_ccco2 </tt> or <tt>./test_boxco2 </tt>
112
 */
113
template <class TypeTag >
Timo Koch's avatar
Timo Koch committed
114
class HeterogeneousProblem : public PorousMediumFlowProblem<TypeTag>
115
{
Timo Koch's avatar
Timo Koch committed
116
117
118
119
120
121
122
123
    using ParentType = PorousMediumFlowProblem<TypeTag>;
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
124
125
126
127
128
129
130
131
132
133

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

    // copy some indices for convenience
    enum {
        lPhaseIdx = Indices::wPhaseIdx,
134
135
136
        gPhaseIdx = Indices::nPhaseIdx
    };
    enum {
137
        wCompIdx = FluidSystem::wCompIdx,
138
139
140
        nCompIdx = FluidSystem::nCompIdx
    };
    enum {
141
        BrineIdx = FluidSystem::BrineIdx,
142
143
144
        CO2Idx = FluidSystem::CO2Idx
    };
    enum {
145
146
147
148
        conti0EqIdx = Indices::conti0EqIdx,
        contiCO2EqIdx = conti0EqIdx + CO2Idx
    };

Timo Koch's avatar
Timo Koch committed
149
150
151
152
153
154
#if !ISOTHERMAL
    enum {
        temperatureIdx = Indices::temperatureIdx,
        energyEqIdx = Indices::energyEqIdx,
    };
#endif
155

Timo Koch's avatar
Timo Koch committed
156
157
158
159
160
161
162
163
164
165
    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
    using Element = typename GridView::template Codim<0>::Entity;
    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
    using CO2 = Dumux::CO2<Scalar, HeterogeneousCO2Tables::CO2Tables>;
166

167
    //! property that defines whether mole or mass fractions are used
168
    static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
169

Timo Koch's avatar
Timo Koch committed
170
171
172
    // the discretization method we are using
    static constexpr auto discMethod = GET_PROP_VALUE(TypeTag, DiscretizationMethod);

173
174
175
176
177
178
179
public:
    /*!
     * \brief The constructor
     *
     * \param timeManager The time manager
     * \param gridView The grid view
     */
Timo Koch's avatar
Timo Koch committed
180
181
182
183
184
185
    HeterogeneousProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
    : ParentType(fvGridGeometry)
    , injectionTop_(1)
    , injectionBottom_(2)
    , dirichletBoundary_(3)
    , noFlowBoundary_(4)
186
    {
Timo Koch's avatar
Timo Koch committed
187
188
189
190
191
192
193
194
195
196
197
198
199
200
        nTemperature_       = getParam<int>("FluidSystem.NTemperature");
        nPressure_          = getParam<int>("FluidSystem.NPressure");
        pressureLow_        = getParam<Scalar>("FluidSystem.PressureLow");
        pressureHigh_       = getParam<Scalar>("FluidSystem.PressureHigh");
        temperatureLow_     = getParam<Scalar>("FluidSystem.TemperatureLow");
        temperatureHigh_    = getParam<Scalar>("FluidSystem.TemperatureHigh");
        depthBOR_           = getParam<Scalar>("Problem.DepthBOR");
        name_               = getParam<std::string>("Problem.Name");
        injectionRate_      = getParam<Scalar>("Problem.InjectionRate");

#if !ISOTHERMAL
        injectionPressure_ = getParam<Scalar>("Problem.InjectionPressure");
        injectionTemperature_ = getParam<Scalar>("Problem.InjectionTemperature");
#endif
201
202

        // set the spatial parameters by reading the DGF grid file
Timo Koch's avatar
Timo Koch committed
203
        this->spatialParams().getParamsFromGrid();
204
205
206
207
208
209
210
211

        // initialize the tables of the fluid system
        FluidSystem::init(/*Tmin=*/temperatureLow_,
                          /*Tmax=*/temperatureHigh_,
                          /*nT=*/nTemperature_,
                          /*pmin=*/pressureLow_,
                          /*pmax=*/pressureHigh_,
                          /*np=*/nPressure_);
212

213
214
215
        //stating in the console whether mole or mass fractions are used
        if(useMoles)
            std::cout<<"problem uses mole fractions"<<std::endl;
216
        else
217
            std::cout<<"problem uses mass fractions"<<std::endl;
218

Timo Koch's avatar
Timo Koch committed
219
220
        // precompute the boundary types for the box method from the cell-centered boundary types
        scvfToScvBoundaryTypes_.computeBoundaryTypes(*this);
221
222
    }

223
224
225
226
227
    /*!
     * \brief Append all quantities of interest which can be derived
     *        from the solution of the current time step to the VTK
     *        writer.
     */
Timo Koch's avatar
Timo Koch committed
228
229
    template<class VTKWriter>
    void addFieldsToWriter(VTKWriter& vtk)
230
    {
Timo Koch's avatar
Timo Koch committed
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
        const auto numElements = this->fvGridGeometry().gridView().size(0);
        const auto numDofs = this->fvGridGeometry().numDofs();

        vtkKxx_.resize(numElements);
        vtkPorosity_.resize(numElements);
        vtkBoxVolume_.resize(numDofs, 0.0);

        vtk.addField(vtkKxx_, "Kxx");
        vtk.addField(vtkPorosity_, "cellwisePorosity");
        vtk.addField(vtkBoxVolume_, "boxVolume");

#if !ISOTHERMAL
        vtk.addVolumeVariable([](const VolumeVariables& v){ return v.enthalpy(Indices::wPhaseIdx); }, "enthalpyW");
        vtk.addVolumeVariable([](const VolumeVariables& v){ return v.enthalpy(Indices::nPhaseIdx); }, "enthalpyN");
#else
        vtkTemperature_.resize(numDofs, 0.0);
        vtk.addField(vtkTemperature_, "temperature");
#endif

        const auto& gridView = this->fvGridGeometry().gridView();
        for (const auto& element : elements(gridView))
        {
            const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
            auto fvGeometry = localView(this->fvGridGeometry());
            fvGeometry.bindElement(element);

            for (const auto& scv : scvs(fvGeometry))
            {
                const auto dofIdxGlobal = scv.dofIndex();
                vtkBoxVolume_[dofIdxGlobal] += scv.volume();
#if ISOTHERMAL
                vtkTemperature_[dofIdxGlobal] = initialTemperatureField_(scv.dofPosition());
#endif
            }

            vtkKxx_[eIdx] = this->spatialParams().permeability(eIdx);
            vtkPorosity_[eIdx] = this->spatialParams().porosity(eIdx);
        }
269
270
    }

271
272
273
274
275
276
277
278
279
280
    /*!
     * \name Problem parameters
     */
    // \{

    /*!
     * \brief The problem name.
     *
     * This is used as a prefix for files generated by the simulation.
     */
Timo Koch's avatar
Timo Koch committed
281
    const std::string& name() const
282
283
284
285
286
    { return name_; }

    /*!
     * \brief Returns the temperature within the domain.
     *
287
288
     * \param globalPos The position
     *
289
     * This problem assumes a geothermal gradient with
290
     * a surface temperature of 10 degrees Celsius.
291
     */
292
    Scalar temperatureAtPos(const GlobalPosition &globalPos) const
Timo Koch's avatar
Timo Koch committed
293
    { return initialTemperatureField_(globalPos); }
294
295
296
297
298
299
300
301
302
303
304
305

    // \}

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

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     *
Timo Koch's avatar
Timo Koch committed
306
307
     * \param element The finite element
     * \param scv The sub control volume
308
     */
Timo Koch's avatar
Timo Koch committed
309
310
311
    BoundaryTypes boundaryTypes(const Element &element,
                                const SubControlVolume &scv) const
    { return scvfToScvBoundaryTypes_.boundaryTypes(scv); }
312
313
314
315
316

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     *
Timo Koch's avatar
Timo Koch committed
317
318
     * \param element The finite element
     * \param scvf The sub control volume face
319
     */
Timo Koch's avatar
Timo Koch committed
320
321
    BoundaryTypes boundaryTypes(const Element &element,
                                const SubControlVolumeFace &scvf) const
322
    {
Timo Koch's avatar
Timo Koch committed
323
324
325
        BoundaryTypes bcTypes;
        const auto boundaryId = scvf.boundaryFlag();

326
        if (boundaryId < 1 || boundaryId > 4)
Timo Koch's avatar
Timo Koch committed
327
328
            DUNE_THROW(Dune::InvalidStateException, "Invalid boundary ID: " << boundaryId);

329
        if (boundaryId == dirichletBoundary_)
Timo Koch's avatar
Timo Koch committed
330
            bcTypes.setAllDirichlet();
331
        else
Timo Koch's avatar
Timo Koch committed
332
333
334
            bcTypes.setAllNeumann();

        return bcTypes;
335
336
337
    }

    /*!
338
339
     * \brief Evaluates the boundary conditions for a Dirichlet
     *        boundary segment
340
     *
Timo Koch's avatar
Timo Koch committed
341
     * \param returns the Dirichlet values for the conservation equations in
342
     *               \f$ [ \textnormal{unit of primary variable} ] \f$
343
     * \param globalPos The global position
344
     */
Timo Koch's avatar
Timo Koch committed
345
346
    PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
    { return initial_(globalPos); }
347

348
    /*!
Timo Koch's avatar
Timo Koch committed
349
     * \brief Evaluate the boundary conditions for a neumann
350
351
     *        boundary segment.
     *
Timo Koch's avatar
Timo Koch committed
352
353
354
     * This is the method for the case where the Neumann condition is
     * potentially solution dependent and requires some quantities that
     * are specific to the fully-implicit method.
355
     *
Timo Koch's avatar
Timo Koch committed
356
357
358
359
360
361
     * \param values The neumann values for the conservation equations in units of
     *                 \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$
     * \param element The finite element
     * \param fvGeometry The finite-volume geometry
     * \param elemVolVars All volume variables for the element
     * \param scvf The sub control volume face
362
     *
Timo Koch's avatar
Timo Koch committed
363
364
365
     * For this method, the \a values parameter stores the flux
     * in normal direction of each phase. Negative values mean influx.
     * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
366
     */
Timo Koch's avatar
Timo Koch committed
367
368
369
370
    NeumannFluxes neumann(const Element& element,
                          const FVElementGeometry& fvGeometry,
                          const ElementVolumeVariables& elemVolvars,
                          const SubControlVolumeFace& scvf) const
371
    {
Timo Koch's avatar
Timo Koch committed
372
        const auto boundaryId = scvf.boundaryFlag();
373

Timo Koch's avatar
Timo Koch committed
374
375
        NeumannFluxes fluxes(0.0);
         // kg/(m^2*s) or mole/(m^2*s) depending on useMoles
376
377
        if (boundaryId == injectionBottom_)
        {
Timo Koch's avatar
Timo Koch committed
378
379
380
381
382
383
            fluxes[contiCO2EqIdx] = useMoles ? -injectionRate_/FluidSystem::molarMass(nCompIdx) : -injectionRate_;
#if !ISOTHERMAL
            // energy fluxes are always mass specific
            fluxes[energyEqIdx] = -injectionRate_/*kg/(m^2 s)*/*CO2::gasEnthalpy(
                                    injectionTemperature_, injectionPressure_)/*J/kg*/; // W/(m^2)
#endif
384
        }
Timo Koch's avatar
Timo Koch committed
385
386

        return fluxes;
387
388
389
390
391
392
393
394
395
396
    }

    // \}

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

    /*!
Timo Koch's avatar
Timo Koch committed
397
     * \brief Evaluates the initial values at a position
398
     *
Timo Koch's avatar
Timo Koch committed
399
400
     * \returns the initial values for the conservation equations in
     *           \f$ [ \textnormal{unit of primary variables} ] \f$
401
     * \param globalPos The global position
402
     */
Timo Koch's avatar
Timo Koch committed
403
    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
404
    {
Timo Koch's avatar
Timo Koch committed
405
        return initial_(globalPos);
406
407
408
409
410
    }

    // \}

private:
411
412
413
414
415
416
417
418
419
    /*!
     * \brief Evaluates the initial values for a control volume
     *
     * The internal method for the initial condition
     *
     * \param values Stores the initial values for the conservation equations in
     *               \f$ [ \textnormal{unit of primary variables} ] \f$
     * \param globalPos The global position
     */
Timo Koch's avatar
Timo Koch committed
420
    PrimaryVariables initial_(const GlobalPosition &globalPos) const
421
    {
Timo Koch's avatar
Timo Koch committed
422
423
424
425
426
427
428
429
        PrimaryVariables values(0.0);
        values.setState(Indices::wPhaseOnly);

        const Scalar temp = initialTemperatureField_(globalPos);
        const Scalar densityW = FluidSystem::Brine::liquidDensity(temp, 1e7);

        const Scalar moleFracLiquidCO2 = 0.00;
        const Scalar moleFracLiquidBrine = 1.0 - moleFracLiquidCO2;
430

Timo Koch's avatar
Timo Koch committed
431
432
        const Scalar meanM = FluidSystem::molarMass(BrineIdx)*moleFracLiquidBrine
                             + FluidSystem::molarMass(CO2Idx)*moleFracLiquidCO2;
433

434
435
436
        if(useMoles) // mole-fraction formulation
            values[Indices::switchIdx] = moleFracLiquidCO2;
        else // mass-fraction formulation
Timo Koch's avatar
Timo Koch committed
437
438
439
440
441
442
443
444
            values[Indices::switchIdx] = moleFracLiquidCO2*FluidSystem::molarMass(CO2Idx)/meanM;

        values[Indices::pressureIdx] = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(depthBOR_ - globalPos[dimWorld-1]);

#if !ISOTHERMAL
        values[temperatureIdx] = temp;
#endif
        return values;
445
446
    }

Timo Koch's avatar
Timo Koch committed
447
    Scalar initialTemperatureField_(const GlobalPosition globalPos) const
448
    {
Timo Koch's avatar
Timo Koch committed
449
        return 283.0 + (depthBOR_ - globalPos[dimWorld-1])*0.03;
450
    }
451
452
453
454

    Scalar depthBOR_;
    Scalar injectionRate_;

Timo Koch's avatar
Timo Koch committed
455
456
457
458
#if !ISOTHERMAL
    Scalar injectionPressure_, injectionTemperature_;
#endif

459
460
461
462
463
464
465
466
467
468
469
470
471
    int nTemperature_;
    int nPressure_;

    std::string name_ ;

    Scalar pressureLow_, pressureHigh_;
    Scalar temperatureLow_, temperatureHigh_;

    int injectionTop_;
    int injectionBottom_;
    int dirichletBoundary_;
    int noFlowBoundary_;

Timo Koch's avatar
Timo Koch committed
472
473
474
    // vtk output
    std::vector<Scalar> vtkKxx_, vtkPorosity_, vtkBoxVolume_, vtkTemperature_;
    ScvfToScvBoundaryTypes<BoundaryTypes, discMethod> scvfToScvBoundaryTypes_;
475
};
Timo Koch's avatar
Timo Koch committed
476
477

} //end namespace Dumux
478
479

#endif