parkervangenuchtenzerohysteresis.hh 12.8 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
27
// -*- 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 Implementation of van Genuchten's capillary pressure-saturation relation.
 *
 */
#ifndef DUMUX_PARKERVANGENZEROHYST_3P_HH
#define DUMUX_PARKERVANGENZEROHYST_3P_HH

28
#include <cmath>
29
30
#include <algorithm>

31
#include <dune/common/fvector.hh>
32
#include <dumux/material/fluidmatrixinteractions/fluidmatrixinteraction.hh>
33

34
namespace Dumux::FluidMatrix {
35
36
37
38
39
40
41
42
/*!
 * \brief Implementation of van Genuchten's capillary pressure <->
 *        saturation relation. This class bundles the "raw" curves
 *        as static members and doesn't concern itself converting
 *        absolute to effective saturations and vince versa.
 *
 * \sa VanGenuchten, VanGenuchtenThreephase
 */
43
44
template <class ScalarT>
class ParkerVanGenZeroHyst3P : Adapter<ParkerVanGenZeroHyst3P<ScalarT>, ThreePhasePcKrSw>
45
46
{
public:
47
    using Scalar = ScalarT;
48
49

    /*!
50
51
     * \brief The parameter type
     * \tparam Scalar The scalar type
52
     */
53
    struct Params
54
    {
55
56
57
58
        Params()
        {
            betaGw_ = betaNw_ = betaGn_ = 1.;
        }
59

60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
        Params(Scalar vgAlpha,
               Scalar vgn,
               Dune::FieldVector<Scalar, 4> residualSaturation,
               Scalar betaNw = 1.,
               Scalar betaGn = 1.,
               Scalar betaGw = 1.,
               bool regardSnr = false,
               Scalar TrappedSatN = 1.)
        {
            setVgAlpha(vgAlpha);
            setVgn(vgn);
            setSwr(residualSaturation[0]);
            setSnr(residualSaturation[1]);
            setSgr(residualSaturation[2]);
            setSwrx(residualSaturation[3]);
            setBetaNw(betaNw);
            setBetaGn(betaGn);
            setBetaGw(betaGw);
            setTrappedSatN(TrappedSatN);
        };

        Scalar TrappedSatN() const
        {
            return TrappedSatN_;
        }
85

86
87
88
89
        void setTrappedSatN(Scalar v)
        {
            TrappedSatN_ = v;
        }
90

91
92
93
94
95
96
97
98
        /*!
         * \brief Return the \f$\alpha\f$ shape parameter of van Genuchten's
         *        curve.
         */
        Scalar vgAlpha() const
        {
            return vgAlpha_;
        }
99

100
101
102
103
104
105
106
107
        /*!
         * \brief Set the \f$\alpha\f$ shape parameter of van Genuchten's
         *        curve.
         */
        void setVgAlpha(Scalar v)
        {
            vgAlpha_ = v;
        }
108

109
110
111
112
113
114
115
116
        /*!
         * \brief Return the \f$m\f$ shape parameter of van Genuchten's
         *        curve.
         */
        Scalar vgm() const
        {
            return vgm_;
        }
117

118
119
120
121
122
123
124
125
126
127
        /*!
         * \brief Set the \f$m\f$ shape parameter of van Genuchten's
         *        curve.
         *
         * The \f$n\f$ shape parameter is set to \f$n = \frac{1}{1 - m}\f$
         */
        void setVgm(Scalar m)
        {
            vgm_ = m; vgn_ = 1/(1 - vgm_);
        }
128

129
130
131
132
133
134
135
136
        /*!
         * \brief Return the \f$n\f$ shape parameter of van Genuchten's
         *        curve.
         */
        Scalar vgn() const
        {
            return vgn_;
        }
137

138
139
140
141
142
143
144
145
146
147
        /*!
         * \brief Set the \f$n\f$ shape parameter of van Genuchten's
         *        curve.
         *
         * The \f$n\f$ shape parameter is set to \f$m = 1 - \frac{1}{n}\f$
         */
        void setVgn(Scalar n)
        {
            vgn_ = n; vgm_ = 1 - 1/vgn_;
        }
148

149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
        /*!
         * \brief Return the residual saturation.
         */
        Scalar satResidual(int phaseIdx) const
        {
            switch (phaseIdx)
            {
            case 0:
                return swr_;
                break;
            case 1:
                return snr_;
                break;
            case 2:
                return sgr_;
                break;
            }

            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
        }
169

170
171
172
173
174
175
176
177
178
        /*!
         * \brief Set all residual saturations.
         */
        void setResiduals(Dune::FieldVector<Scalar, 3> residualSaturation)
        {
            setSwr(residualSaturation[0]);
            setSnr(residualSaturation[1]);
            setSgr(residualSaturation[2]);
        }
179
180


181
182
183
184
185
186
187
        /*!
         * \brief Return the residual wetting saturation.
         */
        Scalar swr() const
        {
            return swr_;
        }
188

189
190
191
192
193
194
195
        /*!
         * \brief Set the residual wetting saturation.
         */
        void setSwr(Scalar input)
        {
            swr_ = input;
        }
196

197
198
199
200
201
202
203
        /*!
         * \brief Return the residual nonwetting saturation.
         */
        Scalar snr() const
        {
            return snr_;
        }
204

205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
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
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
        /*!
         * \brief Set the residual nonwetting saturation.
         */
        void setSnr(Scalar input)
        {
            snr_ = input;
        }

        /*!
         * \brief Return the residual gas saturation.
         */
        Scalar sgr() const
        {
            return sgr_;
        }

        /*!
         * \brief Set the residual gas saturation.
         */
        void setSgr(Scalar input)
        {
            sgr_ = input;
        }

        Scalar swrx() const
        {
            return swrx_;
        }

        /*!
         * \brief Set the residual gas saturation.
         */
        void setSwrx(Scalar input)
        {
            swrx_ = input;
        }

        /*!
         * \brief defines the scaling parameters of capillary pressure between the phases (=1 for Gas-Water)
         */
        void setBetaNw(Scalar input)
        {
            betaNw_ = input;
        }

        void setBetaGn(Scalar input)
        {
            betaGn_ = input;
        }

        void setBetaGw(Scalar input)
        {
            betaGw_ = input;
        }

        /*!
         * \brief Return the values for the beta scaling parameters of capillary pressure between the phases
         */
        Scalar betaNw() const
        {
            return betaNw_;
        }

        Scalar betaGn() const
        {
            return betaGn_;
        }

        Scalar betaGw() const
        {
            return betaGw_;
        }

        bool krRegardsSnr() const
        {
            return krRegardsSnr_;
        }

        /*!
         * \brief defines if residual n-phase saturation should be regarded in its relative permeability.
         */
        void setKrRegardsSnr(bool input)
        {
            krRegardsSnr_ = input;
        }

    private:
        Scalar vgAlpha_;
        Scalar vgm_;
        Scalar vgn_;
        Scalar swr_;
        Scalar snr_;
        Scalar sgr_;
        Scalar swrx_;     /* (sw+sn)_r */

        Scalar betaNw_;
        Scalar betaGn_;
        Scalar betaGw_;

        Scalar TrappedSatN_;

        bool krRegardsSnr_ ;
    };

    ParkerVanGenZeroHyst3P(const Params& params)
    : params_(params)
    {}

    Scalar pcgw(const Scalar sw, const Scalar sn) const
    { return 0;}

    Scalar pcnw(const Scalar sw, const Scalar sn) const
    { return 0; }

    Scalar pcgn(const Scalar sw, const Scalar sn) const
    { return 0; }

    Scalar pcAlpha(const Scalar sw, const Scalar sn) const
    { return 1; }

    void updateTappedSn(const Scalar trappedSn)
    { params_.setTrappedSatN(trappedSn); }
327
328
329
330
331
332
333
334
335
336

    /*!
     * \brief The relative permeability for the wetting phase of
     *        the medium implied by van Genuchten's
     *        parameterization.
     *
     * The permeability of water in a 3p system equals the standard 2p description.
     * (see p61. in "Comparison of the Three-Phase Oil Relative Permeability Models"
     * MOJDEH  DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.)
     *
337
     * \param sw saturation saturation of the water phase.
338
339
     * \param sn saturation of the NAPL phase.
     */
340
    Scalar krw(const Scalar sw, const Scalar sn) const
341
342
    {
        //transformation to effective saturation
343
        const Scalar se = (sw - params_.swr()) / (1 - params_.TrappedSatN() - params_.swr());
344

Kai Wendel's avatar
Kai Wendel committed
345
346
347
348
349
        // regularization
        if (se > 1.0)
            return 1.;
        if (se < 0.0)
            return 0.;
350

351
        const Scalar r = 1. - std::pow(1 - std::pow(se, 1/params_.vgm()), params_.vgm());
Kai Wendel's avatar
Kai Wendel committed
352

353
354
355
356
        return std::sqrt(se)*r*r;
    };

    /*!
357
     * \brief The relative permeability for the nonwetting phase
358
359
360
361
362
363
364
365
366
367
     *        after the Model of Parker et al. (1987).
     *
     * See model 7 in "Comparison of the Three-Phase Oil Relative Permeability Models"
     * MOJDEH  DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.
     * or more comprehensive in
     * "Estimation of primary drainage three-phase relative permeability for organic
     * liquid transport in the vadose zone", Leonardo I. Oliveira, Avery H. Demond,
     * Journal of Contaminant Hydrology 66 (2003), 261-285
     *
     *
368
369
     * \param sw saturation saturation of the water phase.
     * \param sn saturation of the NAPL phase.
370
     */
371
    Scalar krn(const Scalar sw, const Scalar sn) const
372
    {
373
374
        Scalar swe = std::min((sw - params_.swr()) / (1 - params_.TrappedSatN()- params_.swr()), 1.);
        Scalar ste = std::min((sw +  sn - params_.swr()) / (1 - params_.swr()), 1.);
375
376

        // regularization
Kai Wendel's avatar
Kai Wendel committed
377
378
379
380
381
382
        if (swe <= 0.0)
            swe = 0.;
        if (ste <= 0.0)
            ste = 0.;
        if (ste - swe <= 0.0)
            return 0.;
383
384

        Scalar krn_;
385
386
        krn_ = std::pow(1 - std::pow(swe, 1/params_.vgm()), params_.vgm());
        krn_ -= std::pow(1 - std::pow(ste, 1/params_.vgm()), params_.vgm());
387
388
        krn_ *= krn_;

389
        if (params_.krRegardsSnr())
390
391
        {
            // regard Snr in the permeability of the n-phase, see Helmig1997
392
            const Scalar resIncluded = std::max(std::min((sn - params_.snr()/ (1-params_.swr())), 1.), 0.);
393
394
            krn_ *= std::sqrt(resIncluded );
        }
Kai Wendel's avatar
Kai Wendel committed
395

396
        else
397
            krn_ *= std::sqrt(sn / (1 - params_.swr()));   // Hint: (ste - swe) = sn / (1-Srw)
398
399
400
401
402

        return krn_;
    };

    /*!
403
     * \brief The relative permeability for the nonwetting phase
404
405
406
407
408
409
410
     *        of the medium implied by van Genuchten's
     *        parameterization.
     *
     * The permeability of gas in a 3p system equals the standard 2p description.
     * (see p61. in "Comparison of the Three-Phase Oil Relative Permeability Models"
     * MOJDEH  DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.)
     *
411
     * \param sw saturation saturation of the water phase.
412
413
     * \param sn saturation of the NAPL phase.
     */
414
    Scalar krg(const Scalar sw, const Scalar sn) const
415
    {
416
417
        const Scalar st = sw + sn;
        const Scalar se = std::min(((1-st) - params_.sgr()) / (1 - params_.sgr()), 1.);
418

Kai Wendel's avatar
Kai Wendel committed
419
420
421
422
423
        // regularization
        if (se > 1.0)
            return 0.0;
        if (se < 0.0)
            return 1.0;
424
425

        Scalar scalFact = 1.;
Kai Wendel's avatar
Kai Wendel committed
426

427
        if (st<=0.1)
428
        {
429
          scalFact = (st - params_.sgr())/(0.1 - params_.sgr());
Kai Wendel's avatar
Kai Wendel committed
430
431
          if (scalFact < 0.)
                scalFact = 0.;
432
433
        }

434
        return scalFact * std::pow(1 - se, 1.0/3.) * std::pow(1 - std::pow(se, 1/params_.vgm()), 2*params_.vgm());
435
436
437
438
    };

    /*!
     * \brief The relative permeability for a phase.
439
     * \param sw saturation saturation of the water phase.
440
441
442
     * \param sn saturation of the NAPL phase.
     * \param phase indicator, The saturation of all phases.
     */
443
    Scalar kr(const int phaseIdx, const Scalar sw, const Scalar sn) const
444
    {
445
        switch (phaseIdx)
446
447
        {
        case 0:
448
            return krw(sw, sn);
449
450
            break;
        case 1:
451
            return krn(sw, sn);
452
453
            break;
        case 2:
454
            return krg(sw, sn);
455
456
            break;
        }
Kai Wendel's avatar
Kai Wendel committed
457

458
459
460
        return 0;
    };

461
462
private:
    Params params_;
463
464
};

465
} // end namespace Dumux::FluidMatrix
Kai Wendel's avatar
Kai Wendel committed
466

467
#endif