spatialparams.hh 17.7 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:
3
/*****************************************************************************
4
 *   See the file COPYING for full copying permissions.                      *
5
6
7
 *                                                                           *
 *   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    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
9
10
11
12
 *   (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
 *   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/>.   *
 *****************************************************************************/
19
/*!
Katharina Heck's avatar
Katharina Heck committed
20
21
 * \file
 * \ingroup MPNCTests
22
 * \brief Spatial parameters for the kinetic test-case of the mpnc model.
23
 *
24
 * "Poor-mans" coupling of free-flow and porous medium.
25
 */
26
27
28
#ifndef DUMUX_EVAPORATION_ATMOSPHERE_SPATIALPARAMS_HH
#define DUMUX_EVAPORATION_ATMOSPHERE_SPATIALPARAMS_HH

Bernd Flemisch's avatar
Bernd Flemisch committed
29
#include <dumux/porousmediumflow/properties.hh>
30
#include <dumux/material/spatialparams/fvnonequilibrium.hh>
31
32
33

#include <dumux/material/fluidmatrixinteractions/fluidmatrixinteraction.hh>
#include <dumux/material/fluidmatrixinteractions/2p/brookscorey.hh>
34
#include <dumux/material/fluidmatrixinteractions/mp/mpadapter.hh>
35

36
#include <dumux/common/parameters.hh>
37

38
// material laws for interfacial area
39
40
41
42
#include <dumux/material/fluidmatrixinteractions/2p/interfacialarea/polynomial2ndorder.hh>
#include <dumux/material/fluidmatrixinteractions/2p/interfacialarea/pcmax.hh>
#include <dumux/material/fluidmatrixinteractions/2p/interfacialarea/exponentialcubic.hh>
#include <dumux/material/fluidmatrixinteractions/2p/interfacialarea/interfacialarea.hh>
43

44
45
namespace Dumux {

46
47
/**
 * \brief Definition of the spatial parameters for the evaporation atmosphere Problem (using a "poor man's coupling")
48
 */
49
template<class GridGeometry, class Scalar>
50
class EvaporationAtmosphereSpatialParams
51
52
: public FVNonEquilibriumSpatialParams<GridGeometry, Scalar,
                                       EvaporationAtmosphereSpatialParams<GridGeometry, Scalar>>
53
{
54
55
    using GridView = typename GridGeometry::GridView;
    using FVElementGeometry = typename GridGeometry::LocalView;
56
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
57
    using Element = typename GridView::template Codim<0>::Entity;
58
59
    using ThisType = EvaporationAtmosphereSpatialParams<GridGeometry, Scalar>;
    using ParentType = FVNonEquilibriumSpatialParams<GridGeometry, Scalar, ThisType>;
60

61
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
62

63
    static constexpr auto dimWorld = GridView::dimensionworld;
64
65
66
67
68
69
70
71
72
73
74
75

    using PcKrSwCurve = FluidMatrix::BrooksCoreyDefault<Scalar>;

    using NonwettingSolidInterfacialArea = FluidMatrix::InterfacialArea<Scalar,
                                                                        FluidMatrix::InterfacialAreaExponentialCubic,
                                                                        FluidMatrix::NonwettingSolidInterfacialAreaPcSw>;
    using WettingSolidInterfacialArea = FluidMatrix::InterfacialArea<Scalar,
                                                                     FluidMatrix::InterfacialAreaPolynomialSecondOrder,
                                                                     FluidMatrix::WettingSolidInterfacialAreaPcSw>;
    using WettingNonwettingInterfacialArea = FluidMatrix::InterfacialArea<Scalar,
                                                                          FluidMatrix::InterfacialAreaPcMax,
                                                                          FluidMatrix::WettingNonwettingInterfacialAreaPcSw>;
76
public:
77
    //! Export the type used for the permeability
78
79
    using PermeabilityType = Scalar;

80
81
    EvaporationAtmosphereSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
    : ParentType(gridGeometry)
82
    {
83
84
        heightPM_ = getParam<std::vector<Scalar>>("Grid.Positions1")[1];
        heightDomain_ = getParam<std::vector<Scalar>>("Grid.Positions1")[2];
85

86
87
        porosityPM_ = getParam<Scalar>("SpatialParams.PorousMedium.porosity");
        intrinsicPermeabilityPM_ = getParam<Scalar>("SpatialParams.PorousMedium.permeability");
88

89
90
        porosityFF_ = getParam<Scalar>("SpatialParams.FreeFlow.porosity");
        intrinsicPermeabilityFF_ = getParam<Scalar>("SpatialParams.FreeFlow.permeability");
91

92
93
94
        aWettingNonwettingA1_ = getParam<Scalar>("SpatialParams.soil.aWettingNonwettingA1");
        aWettingNonwettingA2_ = getParam<Scalar>("SpatialParams.soil.aWettingNonwettingA2");
        aWettingNonwettingA3_ = getParam<Scalar>("SpatialParams.soil.aWettingNonwettingA3");
95

96
97
98
        aNonwettingSolidA1_ = getParam<Scalar>("SpatialParams.soil.aNonwettingSolidA1");
        aNonwettingSolidA2_ = getParam<Scalar>("SpatialParams.soil.aNonwettingSolidA2");
        aNonwettingSolidA3_ = getParam<Scalar>("SpatialParams.soil.aNonwettingSolidA3");
99

100
101
102
103
        BCPd_ = getParam<Scalar>("SpatialParams.soil.BCPd");
        BClambda_ = getParam<Scalar>("SpatialParams.soil.BClambda");
        Swr_ = getParam<Scalar>("SpatialParams.soil.Swr");
        Snr_  = getParam<Scalar>("SpatialParams.soil.Snr");
104

105
106
        characteristicLengthFF_ = getParam<Scalar>("SpatialParams.FreeFlow.meanPoreSize");
        characteristicLengthPM_ = getParam<Scalar>("SpatialParams.PorousMedium.meanPoreSize");
107
108
109
110

        factorEnergyTransfer_ = getParam<Scalar>("SpatialParams.PorousMedium.factorEnergyTransfer");
        factorMassTransfer_ = getParam<Scalar>("SpatialParams.PorousMedium.factorMassTransfer");

111
112
113
114
        // PM parameters for Brooks-Corey
        typename PcKrSwCurve::BasicParams paramsPM(BCPd_, BClambda_);
        typename PcKrSwCurve::EffToAbsParams effToAbsParamsPM(Swr_, Snr_);
        pcKrSwCurvePM_ = std::make_unique<PcKrSwCurve>(paramsPM, effToAbsParamsPM);
115

116
117
118
119
        // FF parameters for Brooks-Corey
        typename PcKrSwCurve::BasicParams paramsFF(0/*dummy pe*/, 42/*dummy lambda*/);
        typename PcKrSwCurve::EffToAbsParams effToAbsParamsFF(0.0/*swr*/, 0.0/*snr*/);
        pcKrSwCurveFF_ = std::make_unique<PcKrSwCurve>(paramsFF, effToAbsParamsFF);
120

121
        // determine maximum capillary pressure for wetting-nonwetting surface
122
123
124
125
126
127
128
129
        /* Of course physically there is no such thing as a maximum capillary pressure.
         * The parametrization (VG/BC) goes to infinity and physically there is only one pressure
         * for single phase conditions.
         * Here, this is used for fitting the interfacial area surface: the capillary pressure,
         * where the interfacial area is zero.
         * Technically this value is obtained as the capillary pressure of saturation zero.
         * This value of course only exists for the case of a regularized pc-Sw relation.
         */
130
131
132
133
134
135
136
137
138
        const auto pcMax = pcKrSwCurvePM_->pc(/*sw = */0.0);

        // non-wetting-solid
        using NonwettingSolidInterfacialAreaParams = typename NonwettingSolidInterfacialArea::BasicParams;
        NonwettingSolidInterfacialAreaParams ansParams;
        ansParams.setA1(aNonwettingSolidA1_);
        ansParams.setA2(aNonwettingSolidA2_);
        ansParams.setA3(aNonwettingSolidA3_);
        aNs_ = std::make_unique<NonwettingSolidInterfacialArea>(ansParams, effToAbsParamsPM);
139

140
        // wetting-non wetting: surface which goes to zero on the edges, but is a polynomial
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
        using WettingNonwettingInterfacialAreaParams = typename WettingNonwettingInterfacialArea::BasicParams;
        WettingNonwettingInterfacialAreaParams anwParams;
        anwParams.setPcMax(pcMax);
        anwParams.setA1(aWettingNonwettingA1_);
        anwParams.setA2(aWettingNonwettingA2_);
        anwParams.setA3(aWettingNonwettingA3_);
        aNw_ = std::make_unique<WettingNonwettingInterfacialArea>(anwParams, effToAbsParamsPM);

        // zero-initialized dummys for free flow: no interface where there is only one phase
        aNwFreeFlow_ = std::make_unique<WettingNonwettingInterfacialArea>(WettingNonwettingInterfacialAreaParams(), effToAbsParamsFF);
        aNsFreeFlow_ = std::make_unique<NonwettingSolidInterfacialArea>(NonwettingSolidInterfacialAreaParams(), effToAbsParamsFF);

        // dummy
        using WettingSolidInterfacialAreaParams = typename WettingSolidInterfacialArea::BasicParams;
        WettingSolidInterfacialAreaParams awsParams;
        aWs_ = std::make_unique<WettingSolidInterfacialArea>(awsParams, effToAbsParamsPM); // use ctor without effToAbs
157
    }
158

159
160
    template<class ElementSolution>
    PermeabilityType permeability(const Element& element,
161
                                  const SubControlVolume& scv,
162
                                  const ElementSolution& elemSol) const
163
    {
164
165
        const  auto & globalPos = scv.dofPosition();
        if (inFF_(globalPos))
166
167
168
169
            return intrinsicPermeabilityFF_ ;
        else if (inPM_(globalPos))
            return intrinsicPermeabilityPM_ ;
        else
170
            DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
171
172
    }

Katharina Heck's avatar
Katharina Heck committed
173
174
175
    /*!
     * \brief Function for defining the porosity.
     *        That is possibly solution dependent.
176
     *
Katharina Heck's avatar
Katharina Heck committed
177
178
179
     * \param element The current element
     * \param scv The sub-control volume inside the element.
     * \param elemSol The solution at the dofs connected to the element.
180
     * \return The porosity
Katharina Heck's avatar
Katharina Heck committed
181
182
183
184
185
     */
    template<class ElementSolution>
    Scalar porosity(const Element& element,
                    const SubControlVolume& scv,
                    const ElementSolution& elemSol) const
186
    {
187
        const auto& globalPos =  scv.dofPosition();
188

189
        if (inFF_(globalPos))
190
191
192
193
            return porosityFF_ ;
        else if (inPM_(globalPos))
            return porosityPM_ ;
        else
194
            DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
195
196
    }

197
    /*!
198
     * \brief Returns the fluid-matrix interaction law at a given location
199
200
201
     * \param globalPos A global coordinate vector
     */
    const NonwettingSolidInterfacialArea& nonwettingSolidInterfacialAreaAtPos(const GlobalPosition &globalPos) const
202
    {
203
204
        if (inFF_(globalPos) )
            return *aNsFreeFlow_  ;
205
        else if (inPM_(globalPos))
206
            return *aNs_ ;
207
        else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
208
209
    }

210
    /*!
211
     * \brief Returns the fluid-matrix interaction law at a given location
212
     * \param globalPos A global coordinate vector
213
     */
214
    const WettingSolidInterfacialArea& wettingSolidInterfacialAreaAtPos(const GlobalPosition &globalPos) const
215
    {
216
        DUNE_THROW(Dune::InvalidStateException, "Should not be called in this test.");
217
    }
218

219
    /*!
220
     * \brief Returns the fluid-matrix interaction law at a given location
221
     * \param globalPos A global coordinate vector
222
     */
223
    const WettingNonwettingInterfacialArea& wettingNonwettingInterfacialAreaAtPos(const GlobalPosition &globalPos) const
224
225
    {
        if (inFF_(globalPos) )
226
            return *aNwFreeFlow_  ;
227
        else if (inPM_(globalPos))
228
            return *aNw_ ;
229
        else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
230

231
    }
232

233
    /*!
234
     * \brief Returns the fluid-matrix interaction law at a given location
235
     */
236
    template<class ElementSolution>
237
238
239
240
241
242
243
244
    auto fluidMatrixInteraction(const Element& element,
                                                         const SubControlVolume& scv,
                                                         const ElementSolution& elemSol) const
    {
        return fluidMatrixInteractionAtPos(scv.dofPosition());
    }

    /*!
245
     * \brief Returns the fluid-matrix interaction law at a given location
246
247
     */
    auto fluidMatrixInteractionAtPos(const GlobalPosition &globalPos) const
248
    {
249
        if (inFF_(globalPos))
250
            return makeFluidMatrixInteraction(FluidMatrix::MPAdapter(*pcKrSwCurveFF_), *aNsFreeFlow_, *aNwFreeFlow_, *aWs_);
251
        else if (inPM_(globalPos))
252
            return makeFluidMatrixInteraction(FluidMatrix::MPAdapter(*pcKrSwCurvePM_), *aNs_, *aNw_, *aWs_);
253
        else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
254
    }
255

256
257
258
259
    /*!
     * \brief Returns the characteristic length for the mass transfer.
     * \param globalPos The position in global coordinates.
     */
260
261
262
263
264
265
    const Scalar characteristicLengthAtPos(const  GlobalPosition & globalPos) const
    {
        if (inFF_(globalPos) )
            return characteristicLengthFF_ ;
        else if (inPM_(globalPos))
            return characteristicLengthPM_ ;
266
        else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
267
268
    }

269
270
271
272
    /*!
     * \brief Return the pre factor the the energy transfer
     * \param globalPos The position in global coordinates.
     */
273
    const Scalar factorEnergyTransferAtPos(const  GlobalPosition & globalPos) const
274
275
276
277
278
    {
        if (inFF_(globalPos) )
            return factorEnergyTransfer_ ;
        else if (inPM_(globalPos))
            return factorEnergyTransfer_ ;
279
        else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
280
281
    }

282
283
284
285
    /*!
     * \brief Return the pre factor the the mass transfer
     * \param globalPos The position in global coordinates.
     */
286
    const Scalar factorMassTransferAtPos(const  GlobalPosition & globalPos) const
287
288
289
290
291
    {
        if (inFF_(globalPos) )
            return factorMassTransfer_ ;
        else if (inPM_(globalPos))
            return factorMassTransfer_ ;
292
        else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
293
294
    }

295
296
297
298
    /*!
     * \brief Function for defining which phase is to be considered as the wetting phase.
     *
     * \param globalPos The global position
299
     * \return The wetting phase index
300
301
302
303
304
305
306
     */
    template<class FluidSystem>
    int wettingPhaseAtPos(const GlobalPosition& globalPos) const
    {
        return FluidSystem::phase0Idx;
    }

307
308
    /*!
     * \brief Returns whether the tested position is in the porous medium part of the domain
309
     *
310
311
312
313
314
315
316
317
     * This setting ensures that the boundary between the two domains has porous
     * medium properties. This is desirable, because I want to observe the
     * boundary of the porous domain.
     * However, the position has to be the coordinate of the vertex and not the
     * integration point of the boundary flux. If the position is the ip of the
     * Neumann flux this leads to a situation
     * where the vertex belongs to porous medium and there is nonetheless
     * injection over the boundary. That does not work.
318
319
     * -> be careful with neumannAtPos
     */
320
    bool inPM_(const GlobalPosition & globalPos) const
321
    { return ( (globalPos[dimWorld-1] > 0. - eps_) and (globalPos[dimWorld-1] < (heightPM_ + eps_) ) );   }
322
323

    /*!
324
     * \brief Returns whether the tested position is above PM / "free flow" in the domain.
325
     *
326
327
328
329
330
331
332
     * This setting ensures that the boundary between the two domains has porous
     * medium properties. This is desirable, because I want to observe the
     * boundary of the porous domain.
     * However, the position has to be the coordinate of the vertex and not the
     * integration point of the boundary flux. If the position is the ip of the
     * Neumann flux this leads to a situation where the vertex belongs to porous
     * medium and there is nonetheless injection over the boundary.
333
334
335
     * That does not work.
     * -> be careful with neumannAtPos
     */
336
    bool inFF_(const GlobalPosition & globalPos) const
337
    { return ( (globalPos[dimWorld-1] < heightDomain_ + eps_) and (globalPos[dimWorld-1] > (heightPM_ + eps_) ) );   }
338
339
340
341
342
343

    /*! \brief access function for the depth / height of the porous medium */
    const Scalar heightPM() const
    { return heightPM_; }

private:
344
    static constexpr Scalar eps_  = 1e-6;
345
346
347
348
349
350
351
352
353
    Scalar heightDomain_ ;

    // Porous Medium Domain
    Scalar intrinsicPermeabilityPM_ ;
    Scalar porosityPM_ ;
    Scalar heightPM_ ;
    Scalar factorEnergyTransfer_ ;
    Scalar factorMassTransfer_ ;
    Scalar characteristicLengthPM_ ;
354
    std::unique_ptr<PcKrSwCurve> pcKrSwCurvePM_;
355
356
357
358
359

    // Free Flow Domain
    Scalar porosityFF_ ;
    Scalar intrinsicPermeabilityFF_ ;
    Scalar characteristicLengthFF_ ;
360
    std::unique_ptr<PcKrSwCurve> pcKrSwCurveFF_;
361
362

    // interfacial area parameters
363
364
365
    Scalar aWettingNonwettingA1_ ;
    Scalar aWettingNonwettingA2_ ;
    Scalar aWettingNonwettingA3_ ;
366

367
368
369
    Scalar aNonwettingSolidA1_;
    Scalar aNonwettingSolidA2_;
    Scalar aNonwettingSolidA3_;
370
371
372
373
374
375

    // capillary pressures parameters
    Scalar BCPd_ ;
    Scalar BClambda_ ;
    Scalar Swr_ ;
    Scalar Snr_ ;
376
    std::vector<Scalar> gridVector_;
377
378
379
380
381
382
383
384
385

    std::unique_ptr<NonwettingSolidInterfacialArea> aNs_;
    std::unique_ptr<WettingNonwettingInterfacialArea> aNw_;

    std::unique_ptr<WettingNonwettingInterfacialArea> aNwFreeFlow_;
    std::unique_ptr<NonwettingSolidInterfacialArea> aNsFreeFlow_;

    std::unique_ptr<WettingSolidInterfacialArea> aWs_; // dummy, never used

386
387
};

388
} // end namespace Dumux
389

390
#endif