parkervangen3p.hh 13.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
// -*- 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/>.   *
 *****************************************************************************/
/*!
20
21
 * \file
 *
22
 * \brief Implementation of van Genuchten's capillary pressure-saturation relation.
23
 *
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
 */
#ifndef PARKERVANGEN_3P_HH
#define PARKERVANGEN_3P_HH


#include "parkervangen3pparams.hh"

#include <algorithm>


namespace Dumux
{
/*!
 * \ingroup material
 *
 * \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
 */
template <class ScalarT, class ParamsT = ParkerVanGen3PParams<ScalarT> >
class ParkerVanGen3P
{

public:
    typedef ParamsT Params;
    typedef typename Params::Scalar Scalar;

    /*!
     * \brief The capillary pressure-saturation curve.
Alexander Kissinger's avatar
Alexander Kissinger committed
56
57
     * \param params Array of parameters
     * \param sw wetting phase saturation
58
59
60
61
62
63
     *
     */
    static Scalar pc(const Params &params, Scalar sw)
    {
        DUNE_THROW(Dune::NotImplemented, "Capillary pressures for three phases is not so simple! Use pcgn, pcnw, and pcgw");
    }
Alexander Kissinger's avatar
Alexander Kissinger committed
64
65
66
67
68
69
   /*!
     * \brief The capillary pressure-saturation curve copied from MUFTE/pml/constrel3p3cni.c 
     * \param params Array of parameters
     * \param sw wetting phase saturation or sum of wetting phase saturations
     *
     */
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
    static Scalar pcgw(const Params &params, Scalar sw)
    {
    /*
         sw = wetting phase saturation, or,
              sum of wetting phase saturations
         alpha : VanGenuchten-alpha
    this function is just copied from MUFTE/pml/constrel3p3cni.c
    that is why variable names do not yet fulfill Dumux rules, TODO Change */

    Scalar r,se,x,vgm;
    Scalar pc,pcPrime,seRegu;
    Scalar pcvgReg = 0.01;

    se   = (sw-params.swr())/(1.-params.sgr());

    /* Snr  = 0.0;   test version   */

    /* regularization */
    if (se<0.0) se=0.0;
    if (se>1.0) se=1.0;
    vgm = 1.-1./params.vgn();

        if (se>pcvgReg && se<1-pcvgReg)
        {
            r = std::pow(se,-1/vgm);
            x = r-1;
            vgm = 1-vgm;
            x = std::pow(x,vgm);
            r = x/params.vgAlpha();
            return(r);
        }
        else
        {
            /* value and derivative at regularization point */
            if (se<=pcvgReg) seRegu = pcvgReg; else seRegu = 1-pcvgReg;
            pc       = std::pow(std::pow(seRegu,-1/vgm)-1,1/params.vgn())/params.vgAlpha();
106
107
            pcPrime = std::pow(std::pow(seRegu,-1/vgm)-1,1/params.vgn()-1)*std::pow(seRegu,-1/vgm-1)
                      *(-1/vgm)/params.vgAlpha()/(1-params.sgr()-params.swr())/params.vgn();
108
109
110
111
112
113

            /* evaluate tangential */
            r        = (se-seRegu)*pcPrime+pc;
            return(r/params.betaGw());
        }
    }
Alexander Kissinger's avatar
Alexander Kissinger committed
114
115
116
117
118
  /*!
     * \brief The capillary pressure-saturation curve copied from MUFTE/pml/constrel3p3cni.c 
     * \param params Array of parameters
     * \param sw wetting phase saturation or sum of wetting phase saturations
     */
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
    static Scalar pcnw(const Params &params, Scalar sw)
    {
    /*
         sw = wetting phase saturation, or,
              sum of wetting phase saturations
         alpha : VanGenuchten-alpha
    this function is just copied from MUFTE/pml/constrel3p3cni.c
    that is why variable names do not yet fulfill Dumux rules, TODO Change */

    Scalar r,se,x,vgm;
    Scalar pc,pcPrime,seRegu;
    Scalar pcvgReg = 0.01;

    se   = (sw-params.swr())/(1.-params.snr());

    /* Snr  = 0.0;   test version   */

    /* regularization */
    if (se<0.0) se=0.0;
    if (se>1.0) se=1.0;
    vgm = 1.-1./params.vgn();

        if (se>pcvgReg && se<1-pcvgReg)
        {
            r = std::pow(se,-1/vgm);
            x = r-1;
            vgm = 1-vgm;
            x = std::pow(x,vgm);
            r = x/params.vgAlpha();
            return(r);
        }
        else
        {
            /* value and derivative at regularization point */
            if (se<=pcvgReg) seRegu = pcvgReg; else seRegu = 1-pcvgReg;
            pc       = std::pow(std::pow(seRegu,-1/vgm)-1,1/params.vgn())/params.vgAlpha();
155
156
            pcPrime = std::pow(std::pow(seRegu,-1/vgm)-1,1/params.vgn()-1)*std::pow(seRegu,-1/vgm-1)
                      *(-1/vgm)/params.vgAlpha()/(1-params.snr()-params.swr())/params.vgn();
157
158
159
160
161
162

            /* evaluate tangential */
            r        = (se-seRegu)*pcPrime+pc;
            return(r/params.betaNw());
        }
    }
Alexander Kissinger's avatar
Alexander Kissinger committed
163
164
165
    /*!
     * \brief The capillary pressure-saturation curve copied from MUFTE/pml/constrel3p3cni.c 
     * \param params Array of parameters
166
     * \param St sum of wetting (liquid) phase saturations
Alexander Kissinger's avatar
Alexander Kissinger committed
167
     */
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
    static Scalar pcgn(const Params &params, Scalar St)
    {
    /*
         St = sum of wetting (liquid) phase saturations
         alpha : VanGenuchten-alpha
    this function is just copied from MUFTE/pml/constrel3p3cni.c
    that is why variable names do not yet fulfill Dumux rules, TODO Change */

    Scalar r,se,x,vgm;
    Scalar pc,pcPrime,seRegu;
    Scalar pcvgReg = 0.01;

    se   = (St-params.swrx())/(1.-params.swrx());

    /* Snr  = 0.0;   test version   */

    /* regularization */
    if (se<0.0) se=0.0;
    if (se>1.0) se=1.0;
    vgm = 1.-1./params.vgn();

        if (se>pcvgReg && se<1-pcvgReg)
        {
            r = std::pow(se,-1/vgm);
            x = r-1;
            vgm = 1-vgm;
            x = std::pow(x,vgm);
            r = x/params.vgAlpha();
            return(r);
        }
        else
        {
            /* value and derivative at regularization point */
            if (se<=pcvgReg) seRegu = pcvgReg; else seRegu = 1-pcvgReg;
            pc       = std::pow(std::pow(seRegu,-1/vgm)-1,1/params.vgn())/params.vgAlpha();
203
204
            pcPrime = std::pow(std::pow(seRegu,-1/vgm)-1,1/params.vgn()-1)*std::pow(seRegu,-1/vgm-1)
                      *(-1/vgm)/params.vgAlpha()/(1-params.sgr()-params.swrx())/params.vgn();
205
206
207
208
209
210

            /* evaluate tangential */
            r        = (se-seRegu)*pcPrime+pc;
            return(r/params.betaGn());
        }
    }
Alexander Kissinger's avatar
Alexander Kissinger committed
211
212
213
214
215
 /*!
     * \brief The capillary pressure-saturation curve copied from MUFTE/pml/constrel3p3cni.c 
     * \param params Array of parameters
     * \param sn Non-wetting liquid saturation
     */
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
    static Scalar pcAlpha(const Params &params, Scalar sn)
    {
        /* continuous transition to zero */
        Scalar alpha,sne;

        sne=sn;
        /* regularization */
        if (sne<=0.001) sne=0.0;
        if (sne>=1.0) sne=1.0;

        if (sne>params.snr()) alpha = 1.0;
        else
        {
         if (params.snr()>=0.001) alpha = sne/params.snr();
         else          alpha = 0.0;
        }
        return(alpha);
    }

    /*!
     * \brief The saturation-capillary pressure curve.
Alexander Kissinger's avatar
Alexander Kissinger committed
237
238
     * \param params Array of parameters
     * \param pc Capillary pressure in \f$\mathrm{[Pa]}\f$
239
240
241
242
243
244
245
246
247
248
     *
     */
    static Scalar sw(const Params &params, Scalar pc)
    {
        DUNE_THROW(Dune::NotImplemented, "sw(pc) for three phases not implemented! Do it yourself!");
    }

    /*!
     * \brief Returns the partial derivative of the capillary
     *        pressure to the effective saturation.
Alexander Kissinger's avatar
Alexander Kissinger committed
249
250
     * \param params Array of parameters
     * \param sw Wetting liquid saturation
251
252
253
254
255
256
257
258
259
    */
    static Scalar dpc_dsw(const Params &params, Scalar sw)
    {
        DUNE_THROW(Dune::NotImplemented, "dpc/dsw for three phases not implemented! Do it yourself!");
    }

    /*!
     * \brief Returns the partial derivative of the effective
     *        saturation to the capillary pressure.
Alexander Kissinger's avatar
Alexander Kissinger committed
260
261
     * \param params Array of parameters
     * \param pc Capillary pressure in \f$\mathrm{[Pa]}\f$
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
     */
    static Scalar dsw_dpc(const Params &params, Scalar pc)
    {
        DUNE_THROW(Dune::NotImplemented, "dsw/dpc for three phases not implemented! Do it yourself!");
    }

    /*!
     * \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.)
     *
Alexander Kissinger's avatar
Alexander Kissinger committed
277
278
279
     * \param sn Non-wetting liquid saturation
     * \param sg Gas saturation
     * \param saturation wetting liquid saturation
280
281
282
283
284
285
286
287
288
289
290
291
292
293
     * \param params Array of parameters.
     */
    static Scalar krw(const Params &params,  Scalar saturation, Scalar sn, Scalar sg)
    {

        //transformation to effective saturation
        Scalar se = (saturation - params.swr()) / (1-params.swr());

        /* regularization */
        if(se > 1.0) return 1.;
        if(se < 0.0) return 0.;

        Scalar r = 1. - std::pow(1 - std::pow(se, 1/params.vgm()), params.vgm());
        return std::sqrt(se)*r*r;
294
    }
295
296
297
298
299
300
301
302
303
304
305
306
307

    /*!
     * \brief The relative permeability for the non-wetting phase
     *        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
     *
     *
Alexander Kissinger's avatar
Alexander Kissinger committed
308
309
310
     * \param sw Wetting liquid saturation
     * \param sg Gas saturation
     * \param saturation Non-wetting liquid saturation
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
     * \param params Array of parameters.
     */
    static Scalar krn(const Params &params, Scalar sw, Scalar saturation, Scalar sg)
    {

        Scalar swe = std::min((sw - params.swr()) / (1 - params.swr()), 1.);
        Scalar ste = std::min((sw +  saturation - params.swr()) / (1 - params.swr()), 1.);

        // regularization
        if(swe <= 0.0) swe = 0.;
        if(ste <= 0.0) ste = 0.;
        if(ste - swe <= 0.0) return 0.;

        Scalar krn_;
        krn_ = std::pow(1 - std::pow(swe, 1/params.vgm()), params.vgm());
        krn_ -= std::pow(1 - std::pow(ste, 1/params.vgm()), params.vgm());
        krn_ *= krn_;

        if (params.krRegardsSnr())
        {
            // regard Snr in the permeability of the n-phase, see Helmig1997
            Scalar resIncluded = std::max(std::min((saturation - params.snr()/ (1-params.swr())), 1.), 0.);
            krn_ *= std::sqrt(resIncluded );
        }
        else
            krn_ *= std::sqrt(saturation / (1 - params.swr()));   // Hint: (ste - swe) = sn / (1-Srw)


        return krn_;
340
    }
341
342
343
344
345
346
347
348
349
350
351


    /*!
     * \brief The relative permeability for the non-wetting phase
     *        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.)
     *
Alexander Kissinger's avatar
Alexander Kissinger committed
352
353
354
     * \param sw Wetting liquid saturation
     * \param sn Non-wetting liquid saturation
     * \param saturation Gas saturation
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
     * \param params Array of parameters.
     */
    static Scalar krg(const Params &params, Scalar sw, Scalar sn, Scalar saturation)
    {

        // se = (sw+sn - Sgr)/(1-Sgr)
        Scalar se = std::min(((1-saturation) - params.sgr()) / (1 - params.sgr()), 1.);


        /* regularization */
        if(se > 1.0) return 0.0;
        if(se < 0.0) return 1.0;
        Scalar scalFact = 1.;
        if (saturation<=0.1)
        {
          scalFact = (saturation - params.sgr())/(0.1 - params.sgr());
          if (scalFact < 0.) scalFact = 0.;
        }

        Scalar result = scalFact * std::pow(1 - se, 1.0/3.) * std::pow(1 - std::pow(se, 1/params.vgm()), 2*params.vgm());

        return result;
377
    }
378
379
380

    /*!
     * \brief The relative permeability for a phase.
Alexander Kissinger's avatar
Alexander Kissinger committed
381
382
383
     * \param sw Wetting liquid saturation
     * \param sg Gas saturation
     * \param sn Non-wetting liquid saturation
384
     * \param params Array of parameters.
385
     * \param phaseIdx indicator, The saturation of all phases.
386
     */
387
    static Scalar kr(const Params &params, const int phaseIdx, const Scalar sw, const Scalar sn, const Scalar sg)
388
    {
389
        switch (phaseIdx)
390
391
392
393
394
395
396
397
398
399
400
401
        {
        case 0:
            return krw(params, sw, sn, sg);
            break;
        case 1:
            return krn(params, sw, sn, sg);
            break;
        case 2:
            return krg(params, sw, sn, sg);
            break;
        }
        return 0;
402
    }
403

Alexander Kissinger's avatar
Alexander Kissinger committed
404
   /*!
405
    * \brief the basis for calculating adsorbed NAPL in storage term
Alexander Kissinger's avatar
Alexander Kissinger committed
406
    * \param params Array of parameters
407
408
409
410
    */
   static Scalar bulkDensTimesAdsorpCoeff (const Params &params)
   {
      return params.rhoBulk() * params.KdNAPL();
411
   }
412
413
414
415
};
}

#endif