My Project
RegularizedBrooksCorey.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 REGULARIZED_BROOKS_COREY_HPP
28 #define REGULARIZED_BROOKS_COREY_HPP
29 
30 #include "BrooksCorey.hpp"
32 
34 
35 namespace Opm {
64 template <class TraitsT, class ParamsT = RegularizedBrooksCoreyParams<TraitsT> >
65 class RegularizedBrooksCorey : public TraitsT
66 {
67  typedef ::Opm::BrooksCorey<TraitsT, ParamsT> BrooksCorey;
68 
69 public:
70  typedef TraitsT Traits;
71  typedef ParamsT Params;
72  typedef typename Traits::Scalar Scalar;
73 
75  static const int numPhases = Traits::numPhases;
76  static_assert(numPhases == 2,
77  "The regularized Brooks-Corey capillary pressure law only "
78  "applies to the case of two fluid phases");
79 
82  static const bool implementsTwoPhaseApi = true;
83 
86  static const bool implementsTwoPhaseSatApi = true;
87 
90  static const bool isSaturationDependent = true;
91 
94  static const bool isPressureDependent = false;
95 
98  static const bool isTemperatureDependent = false;
99 
102  static const bool isCompositionDependent = false;
103 
104  static_assert(Traits::numPhases == 2,
105  "The number of fluid phases must be two if you want to use "
106  "this material law!");
107 
118  template <class Container, class FluidState>
119  static void capillaryPressures(Container& values, const Params& params, const FluidState& fs)
120  {
121  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
122 
123  values[Traits::wettingPhaseIdx] = 0.0; // reference phase
124  values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
125  }
126 
131  template <class Container, class FluidState>
132  static void saturations(Container& values, const Params& params, const FluidState& fs)
133  {
134  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
135 
136  values[Traits::wettingPhaseIdx] = Sw<FluidState, Evaluation>(params, fs);
137  values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
138  }
139 
150  template <class Container, class FluidState>
151  static void relativePermeabilities(Container& values, const Params& params, const FluidState& fs)
152  {
153  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
154 
155  values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
156  values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
157  }
158 
173  template <class FluidState, class Evaluation = typename FluidState::Scalar>
174  static Evaluation pcnw(const Params& params, const FluidState& fs)
175  {
176  const auto& Sw = decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
177  return twoPhaseSatPcnw(params, Sw);
178  }
179 
180  template <class Evaluation>
181  static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
182  {
183  const Scalar Sthres = params.pcnwLowSw();
184 
185  if (Sw <= Sthres) {
186  Scalar m = params.pcnwSlopeLow();
187  Scalar pcnw_SwLow = params.pcnwLow();
188  return pcnw_SwLow + m*(Sw - Sthres);
189  }
190  else if (Sw >= 1.0) {
191  Scalar m = params.pcnwSlopeHigh();
192  Scalar pcnw_SwHigh = params.pcnwHigh();
193  return pcnw_SwHigh + m*(Sw - 1.0);
194  }
195 
196  // if the effective saturation is in an 'reasonable'
197  // range, we use the real Brooks-Corey law...
198  return BrooksCorey::twoPhaseSatPcnw(params, Sw);
199  }
200 
207  template <class FluidState, class Evaluation = typename FluidState::Scalar>
208  static Evaluation Sw(const Params& params, const FluidState& fs)
209  {
210  const Evaluation& pC =
211  decay<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
212  - decay<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
213  return twoPhaseSatSw(params, pC);
214  }
215 
216  template <class Evaluation>
217  static Evaluation twoPhaseSatSw(const Params& params, const Evaluation& pcnw)
218  {
219  const Scalar Sthres = params.pcnwLowSw();
220 
221  // calculate the saturation which corrosponds to the
222  // saturation in the non-regularized version of the
223  // Brooks-Corey law. If the input capillary pressure is
224  // smaller than the entry pressure, make sure that we will
225  // regularize.
226  Evaluation Sw = 1.5;
227  if (pcnw >= params.entryPressure())
228  Sw = BrooksCorey::twoPhaseSatSw(params, pcnw);
229 
230  // make sure that the capilary pressure observes a
231  // derivative != 0 for 'illegal' saturations. This is
232  // required for example by newton solvers (if the
233  // derivative calculated numerically) in order to get the
234  // saturation moving to the right direction if it
235  // temporarily is in an 'illegal' range.
236  if (Sw <= Sthres) {
237  // invert the low saturation regularization of pcnw()
238  Scalar m = params.pcnwSlopeLow();
239  Scalar pcnw_SwLow = params.pcnwLow();
240  return Sthres + (pcnw - pcnw_SwLow)/m;
241  }
242  else if (Sw > 1.0) {
243  Scalar m = params.pcnwSlopeHigh();
244  Scalar pcnw_SwHigh = params.pcnwHigh();
245  return 1.0 + (pcnw - pcnw_SwHigh)/m;;
246  }
247 
248  return BrooksCorey::twoPhaseSatSw(params, pcnw);
249  }
250 
255  template <class FluidState, class Evaluation = typename FluidState::Scalar>
256  static Evaluation Sn(const Params& params, const FluidState& fs)
257  { return 1 - Sw<FluidState, Evaluation>(params, fs); }
258 
259  template <class Evaluation>
260  static Evaluation twoPhaseSatSn(const Params& params, const Evaluation& pcnw)
261  { return 1 - twoPhaseSatSw(params, pcnw); }
262 
277  template <class FluidState, class Evaluation = typename FluidState::Scalar>
278  static Evaluation krw(const Params& params, const FluidState& fs)
279  {
280  const auto& Sw = decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
281  return twoPhaseSatKrw(params, Sw);
282  }
283 
284  template <class Evaluation>
285  static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
286  {
287  if (Sw <= 0.0)
288  return 0.0;
289  else if (Sw >= 1.0)
290  return 1.0;
291 
292  return BrooksCorey::twoPhaseSatKrw(params, Sw);
293  }
294 
295  template <class Evaluation>
296  static Evaluation twoPhaseSatKrwInv(const Params& params, const Evaluation& krw)
297  {
298  if (krw <= 0.0)
299  return 0.0;
300  else if (krw >= 1.0)
301  return 1.0;
302 
303  return BrooksCorey::twoPhaseSatKrwInv(params, krw);
304  }
305 
320  template <class FluidState, class Evaluation = typename FluidState::Scalar>
321  static Evaluation krn(const Params& params, const FluidState& fs)
322  {
323  const Evaluation& Sw =
324  1.0 - decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
325  return twoPhaseSatKrn(params, Sw);
326  }
327 
328  template <class Evaluation>
329  static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
330  {
331  if (Sw >= 1.0)
332  return 0.0;
333  else if (Sw <= 0.0)
334  return 1.0;
335 
336  return BrooksCorey::twoPhaseSatKrn(params, Sw);
337  }
338 
339  template <class Evaluation>
340  static Evaluation twoPhaseSatKrnInv(const Params& params, const Evaluation& krn)
341  {
342  if (krn <= 0.0)
343  return 1.0;
344  else if (krn >= 1.0)
345  return 0.0;
346 
347  return BrooksCorey::twoPhaseSatKrnInv(params, krn);
348  }
349 };
350 } // namespace Opm
351 
352 #endif
Implementation of the Brooks-Corey capillary pressure <-> saturation relation.
Parameters that are necessary for the regularization of the Brooks-Corey capillary pressure model.
Class implementing cubic splines.
Implementation of the Brooks-Corey capillary pressure <-> saturation relation.
Definition: BrooksCorey.hpp:54
Implementation of the regularized Brooks-Corey capillary pressure / relative permeability <-> saturat...
Definition: RegularizedBrooksCorey.hpp:66
static Evaluation krn(const Params &params, const FluidState &fs)
Regularized version of the relative permeability of the non-wetting phase of the Brooks-Corey curves.
Definition: RegularizedBrooksCorey.hpp:321
static const bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition: RegularizedBrooksCorey.hpp:82
static const bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition: RegularizedBrooksCorey.hpp:86
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
The capillary pressure-saturation curves depending on absolute saturations.
Definition: RegularizedBrooksCorey.hpp:119
static const bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition: RegularizedBrooksCorey.hpp:94
static Evaluation Sn(const Params &params, const FluidState &fs)
Calculate the non-wetting phase saturations depending on the phase pressures.
Definition: RegularizedBrooksCorey.hpp:256
static const int numPhases
The number of fluid phases.
Definition: RegularizedBrooksCorey.hpp:75
static Evaluation pcnw(const Params &params, const FluidState &fs)
A regularized Brooks-Corey capillary pressure-saturation curve.
Definition: RegularizedBrooksCorey.hpp:174
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
The relative permeability-saturation curves depending on absolute saturations.
Definition: RegularizedBrooksCorey.hpp:151
static const bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition: RegularizedBrooksCorey.hpp:102
static const bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition: RegularizedBrooksCorey.hpp:98
static Evaluation krw(const Params &params, const FluidState &fs)
Regularized version of the relative permeability of the wetting phase of the Brooks-Corey curves.
Definition: RegularizedBrooksCorey.hpp:278
static void saturations(Container &values, const Params &params, const FluidState &fs)
Calculate the saturations of the phases starting from their pressure differences.
Definition: RegularizedBrooksCorey.hpp:132
static const bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition: RegularizedBrooksCorey.hpp:90
static Evaluation Sw(const Params &params, const FluidState &fs)
A regularized Brooks-Corey saturation-capillary pressure curve.
Definition: RegularizedBrooksCorey.hpp:208