My Project
ThreePhaseParkerVanGenuchten.hpp
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef OPM_THREE_PHASE_PARKER_VAN_GENUCHTEN_HPP
28 #define OPM_THREE_PHASE_PARKER_VAN_GENUCHTEN_HPP
29 
31 
33 
34 #include <algorithm>
35 
36 namespace Opm {
48 template <class TraitsT,
49  class ParamsT = ThreePhaseParkerVanGenuchtenParams<TraitsT> >
51 {
52 public:
53  static_assert(TraitsT::numPhases == 3,
54  "The number of phases considered by this capillary pressure "
55  "law is always three!");
56 
57  typedef TraitsT Traits;
58  typedef ParamsT Params;
59  typedef typename Traits::Scalar Scalar;
60 
61  static const int numPhases = 3;
62  static const int wettingPhaseIdx = Traits::wettingPhaseIdx;
63  static const int nonWettingPhaseIdx = Traits::nonWettingPhaseIdx;
64  static const int gasPhaseIdx = Traits::gasPhaseIdx;
65 
68  static const bool implementsTwoPhaseApi = false;
69 
72  static const bool implementsTwoPhaseSatApi = false;
73 
76  static const bool isSaturationDependent = true;
77 
80  static const bool isPressureDependent = false;
81 
84  static const bool isTemperatureDependent = false;
85 
88  static const bool isCompositionDependent = false;
89 
101  template <class ContainerT, class FluidState>
102  static void capillaryPressures(ContainerT& values,
103  const Params& params,
104  const FluidState& fluidState)
105  {
106  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
107 
108  values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, fluidState);
109  values[nonWettingPhaseIdx] = 0;
110  values[wettingPhaseIdx] = - pcnw<FluidState, Evaluation>(params, fluidState);
111  }
112 
122  template <class FluidState, class Evaluation = typename FluidState::Scalar>
123  static Evaluation pcgn(const Params& params, const FluidState& fluidState)
124  {
125  Scalar PC_VG_REG = 0.01;
126 
127  // sum of liquid saturations
128  const auto& St =
129  decay<Evaluation>(fluidState.saturation(wettingPhaseIdx))
130  + decay<Evaluation>(fluidState.saturation(nonWettingPhaseIdx));
131 
132  Evaluation Se = (St - params.Swrx())/(1. - params.Swrx());
133 
134  // regularization
135  if (Se < 0.0)
136  Se=0.0;
137  if (Se > 1.0)
138  Se=1.0;
139 
140  if (Se>PC_VG_REG && Se<1-PC_VG_REG)
141  {
142  const Evaluation& x = pow(Se,-1/params.vgM()) - 1;
143  return pow(x, 1.0 - params.vgM())/params.vgAlpha();
144  }
145 
146  // value and derivative at regularization point
147  Scalar Se_regu;
148  if (Se<=PC_VG_REG)
149  Se_regu = PC_VG_REG;
150  else
151  Se_regu = 1-PC_VG_REG;
152  const Evaluation& x = std::pow(Se_regu,-1/params.vgM())-1;
153  const Evaluation& pc = pow(x, 1.0/params.vgN())/params.vgAlpha();
154  const Evaluation& pc_prime =
155  pow(x, 1/params.vgN()-1)
156  * std::pow(Se_regu,-1/params.vgM()-1)
157  / (-params.vgM())
158  / params.vgAlpha()
159  / (1 - params.Sgr() - params.Swrx())
160  / params.vgN();
161 
162  // evaluate tangential
163  return ((Se-Se_regu)*pc_prime + pc)/params.betaGN();
164  }
165 
175  template <class FluidState, class Evaluation = typename FluidState::Scalar>
176  static Evaluation pcnw(const Params& params, const FluidState& fluidState)
177  {
178  const Evaluation& Sw =
179  decay<Evaluation>(fluidState.saturation(wettingPhaseIdx));
180  Evaluation Se = (Sw-params.Swr())/(1.-params.Snr());
181 
182  Scalar PC_VG_REG = 0.01;
183 
184  // regularization
185  if (Se<0.0)
186  Se=0.0;
187  if (Se>1.0)
188  Se=1.0;
189 
190  if (Se>PC_VG_REG && Se<1-PC_VG_REG) {
191  Evaluation x = pow(Se,-1/params.vgM()) - 1.0;
192  x = pow(x, 1 - params.vgM());
193  return x/params.vgAlpha();
194  }
195 
196  // value and derivative at regularization point
197  Scalar Se_regu;
198  if (Se<=PC_VG_REG)
199  Se_regu = PC_VG_REG;
200  else
201  Se_regu = 1.0 - PC_VG_REG;
202 
203  const Evaluation& x = std::pow(Se_regu,-1/params.vgM())-1;
204  const Evaluation& pc = pow(x, 1/params.vgN())/params.vgAlpha();
205  const Evaluation& pc_prime =
206  pow(x,1/params.vgN()-1)
207  * std::pow(Se_regu, -1.0/params.vgM() - 1)
208  / (-params.vgM())
209  / params.vgAlpha()
210  / (1-params.Snr()-params.Swr())
211  / params.vgN();
212 
213  // evaluate tangential
214  return ((Se-Se_regu)*pc_prime + pc)/params.betaNW();
215  }
216 
221  template <class ContainerT, class FluidState>
222  static void saturations(ContainerT& /*values*/,
223  const Params& /*params*/,
224  const FluidState& /*fluidState*/)
225  { throw std::logic_error("Not implemented: inverse capillary pressures"); }
226 
230  template <class FluidState, class Evaluation = typename FluidState::Scalar>
231  static Evaluation Sg(const Params& /*params*/, const FluidState& /*fluidState*/)
232  { throw std::logic_error("Not implemented: Sg()"); }
233 
237  template <class FluidState, class Evaluation = typename FluidState::Scalar>
238  static Evaluation Sn(const Params& /*params*/, const FluidState& /*fluidState*/)
239  { throw std::logic_error("Not implemented: Sn()"); }
240 
244  template <class FluidState, class Evaluation = typename FluidState::Scalar>
245  static Evaluation Sw(const Params& /*params*/, const FluidState& /*fluidState*/)
246  { throw std::logic_error("Not implemented: Sw()"); }
247 
251  template <class ContainerT, class FluidState>
252  static void relativePermeabilities(ContainerT& values,
253  const Params& params,
254  const FluidState& fluidState)
255  {
256  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
257 
258  values[wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
259  values[nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
260  values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
261  }
262 
271  template <class FluidState, class Evaluation = typename FluidState::Scalar>
272  static Evaluation krw(const Params& params, const FluidState& fluidState)
273  {
274  const Evaluation& Sw =
275  decay<Evaluation>(fluidState.saturation(wettingPhaseIdx));
276  // transformation to effective saturation
277  const Evaluation& Se = (Sw - params.Swr()) / (1-params.Swr());
278 
279  // regularization
280  if(Se > 1.0) return 1.;
281  if(Se < 0.0) return 0.;
282 
283  const Evaluation& r = 1. - pow(1 - pow(Se, 1/params.vgM()), params.vgM());
284  return sqrt(Se)*r*r;
285  }
286 
299  template <class FluidState, class Evaluation = typename FluidState::Scalar>
300  static Evaluation krn(const Params& params, const FluidState& fluidState)
301  {
302  const Evaluation& Sn =
303  decay<Evaluation>(fluidState.saturation(nonWettingPhaseIdx));
304  const Evaluation& Sw =
305  decay<Evaluation>(fluidState.saturation(wettingPhaseIdx));
306  Evaluation Swe = min((Sw - params.Swr()) / (1 - params.Swr()), 1.);
307  Evaluation Ste = min((Sw + Sn - params.Swr()) / (1 - params.Swr()), 1.);
308 
309  // regularization
310  if(Swe <= 0.0) Swe = 0.;
311  if(Ste <= 0.0) Ste = 0.;
312  if(Ste - Swe <= 0.0) return 0.;
313 
314  Evaluation krn_;
315  krn_ = pow(1 - pow(Swe, 1/params.vgM()), params.vgM());
316  krn_ -= pow(1 - pow(Ste, 1/params.vgM()), params.vgM());
317  krn_ *= krn_;
318 
319  if (params.krRegardsSnr())
320  {
321  // regard Snr in the permeability of the non-wetting
322  // phase, see Helmig1997
323  const Evaluation& resIncluded =
324  max(min(Sw - params.Snr() / (1-params.Swr()), 1.0), 0.0);
325  krn_ *= sqrt(resIncluded );
326  }
327  else
328  krn_ *= sqrt(Sn / (1 - params.Swr()));
329 
330  return krn_;
331  }
332 
333 
344  template <class FluidState, class Evaluation = typename FluidState::Scalar>
345  static Evaluation krg(const Params& params, const FluidState& fluidState)
346  {
347  const Evaluation& Sg =
348  decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
349  const Evaluation& Se = min(((1-Sg) - params.Sgr()) / (1 - params.Sgr()), 1.);
350 
351  // regularization
352  if(Se > 1.0)
353  return 0.0;
354  if(Se < 0.0)
355  return 1.0;
356 
357  Evaluation scaleFactor = 1.;
358  if (Sg<=0.1) {
359  scaleFactor = (Sg - params.Sgr())/(0.1 - params.Sgr());
360  if (scaleFactor < 0.)
361  return 0.0;
362  }
363 
364  return scaleFactor
365  * pow(1 - Se, 1.0/3.)
366  * pow(1 - pow(Se, 1/params.vgM()), 2*params.vgM());
367  }
368 };
369 } // namespace Opm
370 
371 #endif
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Specification of the material params for the three-phase van Genuchten capillary pressure model.
Implementation of three-phase capillary pressure and relative permeability relations proposed by Park...
Definition: ThreePhaseParkerVanGenuchten.hpp:51
static Evaluation Sn(const Params &, const FluidState &)
The saturation of the non-wetting (i.e., oil) phase.
Definition: ThreePhaseParkerVanGenuchten.hpp:238
static void capillaryPressures(ContainerT &values, const Params &params, const FluidState &fluidState)
Implements the three phase capillary pressure law proposed by Parker and van Genuchten.
Definition: ThreePhaseParkerVanGenuchten.hpp:102
static const bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition: ThreePhaseParkerVanGenuchten.hpp:88
static void relativePermeabilities(ContainerT &values, const Params &params, const FluidState &fluidState)
The relative permeability of all phases.
Definition: ThreePhaseParkerVanGenuchten.hpp:252
static Evaluation Sg(const Params &, const FluidState &)
The saturation of the gas phase.
Definition: ThreePhaseParkerVanGenuchten.hpp:231
static const bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition: ThreePhaseParkerVanGenuchten.hpp:80
static Evaluation Sw(const Params &, const FluidState &)
The saturation of the wetting (i.e., water) phase.
Definition: ThreePhaseParkerVanGenuchten.hpp:245
static const bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition: ThreePhaseParkerVanGenuchten.hpp:72
static void saturations(ContainerT &, const Params &, const FluidState &)
The saturation-capillary pressure curve.
Definition: ThreePhaseParkerVanGenuchten.hpp:222
static const bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition: ThreePhaseParkerVanGenuchten.hpp:76
static Evaluation pcnw(const Params &params, const FluidState &fluidState)
Capillary pressure between the non-wetting liquid (i.e., oil) and the wetting liquid (i....
Definition: ThreePhaseParkerVanGenuchten.hpp:176
static Evaluation krg(const Params &params, const FluidState &fluidState)
The relative permeability for the non-wetting phase of the medium implied by van Genuchten's paramete...
Definition: ThreePhaseParkerVanGenuchten.hpp:345
static Evaluation pcgn(const Params &params, const FluidState &fluidState)
Capillary pressure between the gas and the non-wetting liquid (i.e., oil) phase.
Definition: ThreePhaseParkerVanGenuchten.hpp:123
static Evaluation krw(const Params &params, const FluidState &fluidState)
The relative permeability for the wetting phase of the medium implied by van Genuchten's parameteriza...
Definition: ThreePhaseParkerVanGenuchten.hpp:272
static Evaluation krn(const Params &params, const FluidState &fluidState)
The relative permeability for the non-wetting phase due to the model of Parker et al.
Definition: ThreePhaseParkerVanGenuchten.hpp:300
static const bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition: ThreePhaseParkerVanGenuchten.hpp:84
static const bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition: ThreePhaseParkerVanGenuchten.hpp:68