dissolutionproblem.hh 16.9 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
24
25
26
// -*- 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 Problem where water is injected in a for flushing precipitated salt clogging a gas reservoir.
 */
#ifndef DUMUX_DISSOLUTION_PROBLEM_HH
#define DUMUX_DISSOLUTION_PROBLEM_HH

27
#include <dumux/porousmediumflow/2pncmin/implicit/model.hh>
28
#include <dumux/porousmediumflow/implicit/problem.hh>
29
#include <dumux/material/fluidsystems/brineair.hh>
30
31
32
33
34

#include "dissolutionspatialparams.hh"

namespace Dumux
{
35

36
37
38
39
40
41
42
43
44
45
template <class TypeTag>
class DissolutionProblem;

namespace Properties
{
NEW_TYPE_TAG(DissolutionProblem, INHERITS_FROM(TwoPNCMin, DissolutionSpatialparams));
NEW_TYPE_TAG(DissolutionBoxProblem, INHERITS_FROM(BoxModel, DissolutionProblem));
NEW_TYPE_TAG(DissolutionCCProblem, INHERITS_FROM(CCModel, DissolutionProblem));

// Set the grid type
46
SET_TYPE_PROP(DissolutionProblem, Grid, Dune::YaspGrid<2>);
47
48

// Set the problem property
49
SET_TYPE_PROP(DissolutionProblem, Problem, DissolutionProblem<TypeTag>);
50
51
52
53
54

// Set fluid configuration
SET_PROP(DissolutionProblem, FluidSystem)
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
55
    typedef FluidSystems::BrineAir<Scalar, H2O<Scalar>, true/*useComplexrelations=*/> type;
56
57
58
};

// Set the spatial parameters
59
SET_TYPE_PROP(DissolutionProblem, SpatialParams, DissolutionSpatialparams<TypeTag>);
60
61
62
63
64
65

// Enable gravity
SET_BOOL_PROP(DissolutionProblem, ProblemEnableGravity, true);

//Set properties here to override the default property settings in the model.
SET_INT_PROP(DissolutionProblem, ReplaceCompEqIdx, 1);
66
SET_INT_PROP(DissolutionProblem, Formulation, TwoPNCFormulation::pnsw);
67
68
69
70
71
72
73
74
75
76
77
78
79
80
}

/*!
 * \ingroup TwoPNCMinModel
 * \ingroup ImplicitTestProblems
 * \brief Problem where water is injected to flush precipitated salt in a gas reservoir clogged due to precipitated salt.
 *
 * The domain is sized 10m times 20m and contains a vertical low-permeable layer of precipitated salt near an extraction well.
 *
 * To flush this precipitated salt, water is injected through the gas extraction well in order to dissolve the precipitated salt increasing the permeability and thereby achieving high gas extraction rates later. Here, the system is assumed to be isothermal.
 * Neumann no-flow boundary condition is applied at the top and bottom boundary and Dirichlet boundary condition is used on the right and left sides.
 * The injected water phase migrates downwards due to increase in density as the precipitated salt dissolves.
 *
 * The model uses mole fractions of dissolved components and volume fractions of precipitated salt as primary variables. Make sure that the according units are used in the problem setup.
81
 *
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
 * This problem uses the \ref TwoPNCMinModel.
 *
 * To run the simulation execute the following line in shell:
 * <tt>./test_box2pncmin</tt>
 */
template <class TypeTag>
class DissolutionProblem : public ImplicitPorousMediaProblem<TypeTag>
{
    typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;

    enum {

99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
        pressureIdx = Indices::pressureIdx,
        switchIdx = Indices::switchIdx, //Saturation
        xlNaClIdx = FluidSystem::NaClIdx,
        precipNaClIdx = FluidSystem::numComponents,

        //Indices of the components
        wCompIdx = FluidSystem::H2OIdx,
        nCompIdx = FluidSystem::AirIdx,
        NaClIdx = FluidSystem::NaClIdx,

        //Indices of the phases
        wPhaseIdx = FluidSystem::wPhaseIdx,
        nPhaseIdx = FluidSystem::nPhaseIdx,
        sPhaseIdx = FluidSystem::sPhaseIdx,

        //Index of the primary component of G and L phase
        conti0EqIdx = Indices::conti0EqIdx,
        contiTotalMassIdx = conti0EqIdx + FluidSystem::AirIdx,
        precipNaClEqIdx = Indices::conti0EqIdx + FluidSystem::numComponents,
        contiWEqIdx = conti0EqIdx + FluidSystem::H2OIdx,

        // Phase State
        wPhaseOnly = Indices::wPhaseOnly,
        nPhaseOnly = Indices::nPhaseOnly,
        bothPhases = Indices::bothPhases,

        // Grid and world dimension
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld,
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
    };


    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 GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;

    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<dim>::Entity Vertex;
    typedef typename GridView::Intersection Intersection;
    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;

public:
143
144
    DissolutionProblem(TimeManager &timeManager, const GridView &gridView)
        : ParentType(timeManager, gridView)
145
    {
146

147
148
149
150
        outerSalinity_          = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, OuterSalinity);
        temperature_            = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, Temperature);
        reservoirPressure_      = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, ReservoirPressure);
        initLiqSaturation_      = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, LiquidSaturation);
151
152
        initPrecipitatedSalt1_  = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InitPrecipitatedSalt1);
        initPrecipitatedSalt2_  = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InitPrecipitatedSalt2);
153
154
155
156
157
158
159

        outerLiqSaturation_     = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, OuterLiqSaturation);
        innerLiqSaturation_     = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InnerLiqSaturation);
        innerSalinity_          = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InnerSalinity);
        innerPressure_          = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InnerPressure);
        outerPressure_          = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, OuterPressure);
        reservoirSaturation_    = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, reservoirSaturation);
160

161
162
163
164
165
166
        nTemperature_           = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, FluidSystem, NTemperature);
        nPressure_              = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, FluidSystem, NPressure);
        pressureLow_            = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FluidSystem, PressureLow);
        pressureHigh_           = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FluidSystem, PressureHigh);
        temperatureLow_         = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FluidSystem, TemperatureLow);
        temperatureHigh_        = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FluidSystem, TemperatureHigh);
167
        name_                   = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
168
        freqMassOutput_         = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqMassOutput);
169
170
        storageLastTimestep_    = 0.0;
        lastMassOutputTime_     = 0.0;
171
172
173

        outfile.open("evaporation.out");
        outfile << "time; evaporationRate" << std::endl;
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189

        FluidSystem::init(/*Tmin=*/temperatureLow_,
                          /*Tmax=*/temperatureHigh_,
                          /*nT=*/nTemperature_,
                          /*pmin=*/pressureLow_,
                          /*pmax=*/pressureHigh_,
                          /*np=*/nPressure_);
    }

    ~DissolutionProblem()
    {
        outfile.close();
    }

    bool shouldWriteOutput() const
    {
190
191
192
        return this->timeManager().timeStepIndex() % 1 == 0 ||
               this->timeManager().episodeWillBeOver() ||
               this->timeManager().willBeFinished();
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
    }


    /*!
     * \name Problem parameters
     */


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

    /*!
     * \brief Returns the temperature within the domain.
     *
     * This problem assumes a temperature of 10 degrees Celsius.
     */
    Scalar temperature() const
215
    { return temperature_; }
216
217
218
219
220
221
222
223
224
225

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

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     */
226
    void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition &globalPos) const
227
    {
228
        const Scalar rmax = this->bBoxMax()[0]; // outerRadius_;
Thomas Fetzer's avatar
Thomas Fetzer committed
229
        const Scalar rmin = this->bBoxMin()[0];
230

231
        // default to Neumann
232
233
        bcTypes.setAllNeumann();

234
        // Constant pressure  at reservoir boundary (Dirichlet condition)
235
236
237
        if(globalPos[0] > rmax - eps_)
            bcTypes.setAllDirichlet();

238
        // Constant pressure at well (Dirichlet condition)
239
240
241
242
243
244
245
246
247
248
        if(globalPos[0] < rmin + eps_)
            bcTypes.setAllDirichlet();
    }

    /*!
     * \brief Evaluate the boundary conditions for a dirichlet
     *        boundary segment.
     *
     * For this method, the \a values parameter stores primary variables.
     */
249
    void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
250
    {
251
252
        const Scalar rmax = this->bBoxMax()[0];
        const Scalar rmin = this->bBoxMin()[0];
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280

        if(globalPos[0] > rmax - eps_)
        {
            values[pressureIdx]   = outerPressure_ ; // Outer boundary pressure bar
            values[switchIdx]     = outerLiqSaturation_; // Saturation outer boundary
            values[xlNaClIdx]     = massTomoleFrac_(outerSalinity_);// mole fraction salt
            values[precipNaClIdx] = 0.0;// precipitated salt
        }

        if(globalPos[0] < rmin + eps_)
        {

            values[pressureIdx]   = innerPressure_ ; // Inner boundary pressure bar
            values[switchIdx]     = innerLiqSaturation_; // Saturation inner boundary
            values[xlNaClIdx]     = massTomoleFrac_(innerSalinity_);// mole fraction salt
            values[precipNaClIdx] = 0.0;// precipitated salt
        }
    }

    /*!
     * \brief Evaluate the boundary conditions for a neumann
     *        boundary segment.
     *
     * For this method, the \a values parameter stores the mass flux
     * in normal direction of each component. Negative values mean
     * influx.
     */
    void neumann(PrimaryVariables &values,
281
282
283
284
285
                 const Element &element,
                 const FVElementGeometry &fvGeometry,
                 const Intersection &is,
                 int scvIdx,
                 int boundaryFaceIdx) const
286
    {
287
        values = 0.0;
288
289
290
291
292
293
294
295
    }

    /*!
     * \brief Evaluate the initial value for a control volume.
     *
     * For this method, the \a values parameter stores primary
     * variables.
     */
Thomas Fetzer's avatar
Thomas Fetzer committed
296
    void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
297
298
299
300
    {
        values[pressureIdx] = reservoirPressure_;
        values[switchIdx]   = initLiqSaturation_;                 // Sl primary variable
        values[xlNaClIdx]   = massTomoleFrac_(outerSalinity_);     // mole fraction
301
        if(globalPos[0] > 5.0 - eps_ && globalPos[0] < 20.0 - eps_)
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
            values[precipNaClIdx] = initPrecipitatedSalt2_; // [kg/m^3]
        else
            values[precipNaClIdx] = initPrecipitatedSalt1_; // [kg/m^3]
    }

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

    /*!
     * \brief Evaluate the source term for all phases within a given
     *        sub-control-volume.
     *
     * For this method, the \a values parameter stores the rate mass
     * of a component is generated or annihilate per volume
     * unit. Positive values mean that mass is created, negative ones
     * mean that it vanishes.
     */
    void solDependentSource(PrimaryVariables &source,
                            const Element &element,
                            const FVElementGeometry &fvGeometry,
324
325
                            int scvIdx,
                            const ElementVolumeVariables &elemVolVars) const
326
327
    {
        source = 0;
328
        const auto& volVars = elemVolVars[scvIdx];
329
330
        Scalar moleFracNaCl_lPhase = volVars.moleFraction(wPhaseIdx, NaClIdx);
        Scalar moleFracNaCl_gPhase = volVars.moleFraction(nPhaseIdx, NaClIdx);
331
332
        Scalar massFracNaCl_Max_lPhase = this->spatialParams().SolubilityLimit();
        Scalar moleFracNaCl_Max_lPhase = massTomoleFrac_(massFracNaCl_Max_lPhase);
333
        Scalar moleFracNaCl_Max_gPhase = moleFracNaCl_Max_lPhase / volVars.pressure(nPhaseIdx);
334
        Scalar saltPorosity = this->spatialParams().porosityMin(element, fvGeometry, scvIdx);
335
336
337

        // liquid phase
        Scalar precipSalt = volVars.porosity() * volVars.molarDensity(wPhaseIdx)
338
339
                                               * volVars.saturation(wPhaseIdx)
                                               * std::abs(moleFracNaCl_lPhase - moleFracNaCl_Max_lPhase);
340

341
342
343
344
345
346
        if (moleFracNaCl_lPhase < moleFracNaCl_Max_lPhase)
            precipSalt *= -1;

        // gas phase
        if (moleFracNaCl_gPhase > moleFracNaCl_Max_gPhase)
            precipSalt += volVars.porosity() * volVars.molarDensity(nPhaseIdx)
347
348
                                             * volVars.saturation(nPhaseIdx)
                                             * std::abs(moleFracNaCl_gPhase - moleFracNaCl_Max_gPhase);
349
350
351
352
353

        // make sure we don't disolve more salt than previously precipitated
        if (precipSalt*this->timeManager().timeStepSize() + volVars.precipitateVolumeFraction(sPhaseIdx)* volVars.molarDensity(sPhaseIdx)< 0)
            precipSalt = - volVars.precipitateVolumeFraction(sPhaseIdx)* volVars.molarDensity(sPhaseIdx)/this->timeManager().timeStepSize();

354
        if (volVars.precipitateVolumeFraction(sPhaseIdx) >= volVars.initialPorosity() - saltPorosity  && precipSalt > 0)
355
356
357
358
359
360
361
362
363
364
            precipSalt = 0;

        source[conti0EqIdx + NaClIdx] += -precipSalt;
        source[precipNaClEqIdx] += precipSalt;

        Valgrind::CheckDefined(source);
    }
    /*!
     * \brief Return the initial phase state inside a control volume.
     */
365
366
367
    int initialPhasePresence(const Element& element,
                             const FVElementGeometry& fvGeometry,
                             int scvIdx) const
368
369
370
371
372
373
374
    {
        return bothPhases;
    }

private:

    /*!
375
376
377
378
     * \brief Returns the molality of NaCl (mol NaCl / kg water) for a given mole fraction
     *
     * \param XlNaCl the XlNaCl [kg NaCl / kg solution]
     */
379
380
    static Scalar massTomoleFrac_(Scalar XlNaCl)
    {
381
382
       const Scalar Mw = 18.015e-3; /* molecular weight of water [kg/mol] */
       const Scalar Ms = 58.44e-3; /* molecular weight of NaCl  [kg/mol] */
383

384
385
       const Scalar X_NaCl = XlNaCl;
       /* XlNaCl: conversion from mass fraction to mol fraction */
386
       auto xlNaCl = -Mw * X_NaCl / ((Ms - Mw) * X_NaCl - Ms);
387
       return xlNaCl;
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
    }

    int nTemperature_;
    int nPressure_;
    int freqMassOutput_;
    PrimaryVariables storageLastTimestep_;
    Scalar lastMassOutputTime_;
    std::string name_;

    Scalar pressureLow_, pressureHigh_;
    Scalar temperatureLow_, temperatureHigh_;
    Scalar outerSalinity_;
    Scalar reservoirPressure_;
    Scalar innerPressure_;
    Scalar outerPressure_;
    Scalar temperature_;
    Scalar initLiqSaturation_;
    Scalar outerLiqSaturation_;
    Scalar innerLiqSaturation_;
    Scalar initPrecipitatedSalt1_;
    Scalar initPrecipitatedSalt2_;
    Scalar innerSalinity_;
    static constexpr Scalar eps_ = 1e-6;
    Scalar reservoirSaturation_;
    std::ofstream outfile;

};
} //end namespace

#endif