My Project
EclHysteresisTwoPhaseLaw.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_ECL_HYSTERESIS_TWO_PHASE_LAW_HPP
28 #define OPM_ECL_HYSTERESIS_TWO_PHASE_LAW_HPP
29 
31 
32 #include <stdexcept>
33 
34 namespace Opm {
35 
41 template <class EffectiveLawT,
42  class ParamsT = EclHysteresisTwoPhaseLawParams<EffectiveLawT> >
43 class EclHysteresisTwoPhaseLaw : public EffectiveLawT::Traits
44 {
45 public:
46  using EffectiveLaw = EffectiveLawT;
47  using EffectiveLawParams = typename EffectiveLaw::Params;
48 
49  using Traits = typename EffectiveLaw::Traits;
50  using Params = ParamsT;
51  using Scalar = typename EffectiveLaw::Scalar;
52 
53  enum { wettingPhaseIdx = Traits::wettingPhaseIdx };
54  enum { nonWettingPhaseIdx = Traits::nonWettingPhaseIdx };
55 
57  static constexpr int numPhases = EffectiveLaw::numPhases;
58  static_assert(numPhases == 2,
59  "The endpoint scaling applies to the nested twophase laws, not to "
60  "the threephase one!");
61 
64  static constexpr bool implementsTwoPhaseApi = true;
65 
66  static_assert(EffectiveLaw::implementsTwoPhaseApi,
67  "The material laws put into EclEpsTwoPhaseLaw must implement the "
68  "two-phase material law API!");
69 
72  static constexpr bool implementsTwoPhaseSatApi = true;
73 
74  static_assert(EffectiveLaw::implementsTwoPhaseSatApi,
75  "The material laws put into EclEpsTwoPhaseLaw must implement the "
76  "two-phase material law saturation API!");
77 
80  static constexpr bool isSaturationDependent = true;
81 
84  static constexpr bool isPressureDependent = false;
85 
88  static constexpr bool isTemperatureDependent = false;
89 
92  static constexpr bool isCompositionDependent = false;
93 
104  template <class Container, class FluidState>
105  static void capillaryPressures(Container& /* values */,
106  const Params& /* params */,
107  const FluidState& /* fs */)
108  {
109  throw std::invalid_argument("The capillaryPressures(fs) method is not yet implemented");
110  }
111 
122  template <class Container, class FluidState>
123  static void relativePermeabilities(Container& /* values */,
124  const Params& /* params */,
125  const FluidState& /* fs */)
126  {
127  throw std::invalid_argument("The pcnw(fs) method is not yet implemented");
128  }
129 
141  template <class FluidState, class Evaluation = typename FluidState::Scalar>
142  static Evaluation pcnw(const Params& /* params */,
143  const FluidState& /* fs */)
144  {
145  throw std::invalid_argument("The pcnw(fs) method is not yet implemented");
146  }
147 
148  template <class Evaluation>
149  static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
150  {
151  // if no pc hysteresis is enabled, use the drainage curve
152  if (!params.config().enableHysteresis() || params.config().pcHysteresisModel() < 0)
153  return EffectiveLaw::twoPhaseSatPcnw(params.drainageParams(), Sw);
154 
155  // Initial imbibition process
156  if (params.initialImb()) {
157  if (Sw >= params.pcSwMic()) {
158  return EffectiveLaw::twoPhaseSatPcnw(params.imbibitionParams(), Sw);
159  }
160  else { // Reversal
161  const Evaluation& F = (1.0/(params.pcSwMic()-Sw+params.curvatureCapPrs())-1.0/params.curvatureCapPrs())
162  / (1.0/(params.pcSwMic()-params.Swcrd()+params.curvatureCapPrs())-1.0/params.curvatureCapPrs());
163 
164  const Evaluation& Pcd = EffectiveLaw::twoPhaseSatPcnw(params.drainageParams(), Sw);
165  const Evaluation& Pci = EffectiveLaw::twoPhaseSatPcnw(params.imbibitionParams(), Sw);
166  const Evaluation& pc_Killough = Pci+F*(Pcd-Pci);
167 
168  return pc_Killough;
169  }
170  }
171 
172  // Initial drainage process
173  if (Sw <= params.pcSwMdc())
174  return EffectiveLaw::twoPhaseSatPcnw(params.drainageParams(), Sw);
175 
176  // Reversal
177  Scalar Swma = 1.0-params.Sncrt();
178  if (Sw >= Swma) {
179  const Evaluation& Pci = EffectiveLaw::twoPhaseSatPcnw(params.imbibitionParams(), Sw);
180  return Pci;
181  }
182  else {
183  Scalar pciwght = params.pcWght(); // Align pci and pcd at Swir
184  //const Evaluation& SwEff = params.Swcri()+(Sw-params.pcSwMdc())/(Swma-params.pcSwMdc())*(Swma-params.Swcri());
185  const Evaluation& SwEff = Sw; // This is Killough 1976, Gives significantly better fit compared to benchmark then the above "scaling"
186  const Evaluation& Pci = pciwght*EffectiveLaw::twoPhaseSatPcnw(params.imbibitionParams(), SwEff);
187 
188  const Evaluation& Pcd = EffectiveLaw::twoPhaseSatPcnw(params.drainageParams(), Sw);
189 
190  if (Pci == Pcd)
191  return Pcd;
192 
193  const Evaluation& F = (1.0/(Sw-params.pcSwMdc()+params.curvatureCapPrs())-1.0/params.curvatureCapPrs())
194  / (1.0/(Swma-params.pcSwMdc()+params.curvatureCapPrs())-1.0/params.curvatureCapPrs());
195 
196  const Evaluation& pc_Killough = Pcd+F*(Pci-Pcd);
197 
198  return pc_Killough;
199  }
200 
201  return 0.0;
202  }
203 
207  template <class Container, class FluidState>
208  static void saturations(Container& /* values */,
209  const Params& /* params */,
210  const FluidState& /* fs */)
211  {
212  throw std::invalid_argument("The saturations(fs) method is not yet implemented");
213  }
214 
219  template <class FluidState, class Evaluation = typename FluidState::Scalar>
220  static Evaluation Sw(const Params& /* params */,
221  const FluidState& /* fs */)
222  {
223  throw std::invalid_argument("The Sw(fs) method is not yet implemented");
224  }
225 
226  template <class Evaluation>
227  static Evaluation twoPhaseSatSw(const Params& /* params */,
228  const Evaluation& /* pc */)
229  {
230  throw std::invalid_argument("The twoPhaseSatSw(pc) method is not yet implemented");
231  }
232 
237  template <class FluidState, class Evaluation = typename FluidState::Scalar>
238  static Evaluation Sn(const Params& /* params */,
239  const FluidState& /* fs */)
240  {
241  throw std::invalid_argument("The Sn(pc) method is not yet implemented");
242  }
243 
244  template <class Evaluation>
245  static Evaluation twoPhaseSatSn(const Params& /* params */,
246  const Evaluation& /* pc */)
247  {
248  throw std::invalid_argument("The twoPhaseSatSn(pc) method is not yet implemented");
249  }
250 
260  template <class FluidState, class Evaluation = typename FluidState::Scalar>
261  static Evaluation krw(const Params& /* params */,
262  const FluidState& /* fs */)
263  {
264  throw std::invalid_argument("The krw(fs) method is not yet implemented");
265  }
266 
267  template <class Evaluation>
268  static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
269  {
270  // if no relperm hysteresis is enabled, use the drainage curve
271  if (!params.config().enableHysteresis() || params.config().krHysteresisModel() < 0)
272  return EffectiveLaw::twoPhaseSatKrw(params.drainageParams(), Sw);
273 
274  if (params.config().krHysteresisModel() == 0 || params.config().krHysteresisModel() == 2)
275  // use drainage curve for wetting phase
276  return EffectiveLaw::twoPhaseSatKrw(params.drainageParams(), Sw);
277 
278  // use imbibition curve for wetting phase
279  assert(params.config().krHysteresisModel() == 1 || params.config().krHysteresisModel() == 3);
280  return EffectiveLaw::twoPhaseSatKrw(params.imbibitionParams(), Sw);
281  }
282 
286  template <class FluidState, class Evaluation = typename FluidState::Scalar>
287  static Evaluation krn(const Params& /* params */,
288  const FluidState& /* fs */)
289  {
290  throw std::invalid_argument("The krn(fs) method is not yet implemented");
291  }
292 
293  template <class Evaluation>
294  static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
295  {
296  // if no relperm hysteresis is enabled, use the drainage curve
297  if (!params.config().enableHysteresis() || params.config().krHysteresisModel() < 0)
298  return EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Sw);
299 
300  // if it is enabled, use either the drainage or the imbibition curve. if the
301  // imbibition curve is used, the saturation must be shifted.
302  if (Sw <= params.krnSwMdc())
303  return EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Sw);
304 
305  if (params.config().krHysteresisModel() <= 1) { //Carlson
306  return EffectiveLaw::twoPhaseSatKrn(params.imbibitionParams(),
307  Sw + params.deltaSwImbKrn());
308  }
309 
310  // Killough
311  assert(params.config().krHysteresisModel() == 2 || params.config().krHysteresisModel() == 3);
312  Evaluation Snorm = params.Sncri()+(1.0-Sw-params.Sncrt())*(params.Snmaxd()-params.Sncri())/(params.Snhy()-params.Sncrt());
313  return params.krnWght()*EffectiveLaw::twoPhaseSatKrn(params.imbibitionParams(),1.0-Snorm);
314  }
315 };
316 
317 } // namespace Opm
318 
319 #endif
A default implementation of the parameters for the material law which implements the ECL relative per...
This material law implements the hysteresis model of the ECL file format.
Definition: EclHysteresisTwoPhaseLaw.hpp:44
static constexpr bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition: EclHysteresisTwoPhaseLaw.hpp:80
static Evaluation krn(const Params &, const FluidState &)
The relative permeability of the non-wetting phase.
Definition: EclHysteresisTwoPhaseLaw.hpp:287
static void capillaryPressures(Container &, const Params &, const FluidState &)
The capillary pressure-saturation curves depending on absolute saturations.
Definition: EclHysteresisTwoPhaseLaw.hpp:105
static constexpr bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition: EclHysteresisTwoPhaseLaw.hpp:64
static constexpr bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition: EclHysteresisTwoPhaseLaw.hpp:92
static Evaluation Sw(const Params &, const FluidState &)
Calculate wetting liquid phase saturation given that the rest of the fluid state has been initialized...
Definition: EclHysteresisTwoPhaseLaw.hpp:220
static Evaluation Sn(const Params &, const FluidState &)
Calculate non-wetting liquid phase saturation given that the rest of the fluid state has been initial...
Definition: EclHysteresisTwoPhaseLaw.hpp:238
static void saturations(Container &, const Params &, const FluidState &)
The saturation-capillary pressure curves.
Definition: EclHysteresisTwoPhaseLaw.hpp:208
static Evaluation krw(const Params &, const FluidState &)
The relative permeability for the wetting phase.
Definition: EclHysteresisTwoPhaseLaw.hpp:261
static constexpr int numPhases
The number of fluid phases.
Definition: EclHysteresisTwoPhaseLaw.hpp:57
static constexpr bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition: EclHysteresisTwoPhaseLaw.hpp:84
static constexpr bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition: EclHysteresisTwoPhaseLaw.hpp:72
static constexpr bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition: EclHysteresisTwoPhaseLaw.hpp:88
static Evaluation pcnw(const Params &, const FluidState &)
The capillary pressure-saturation curve.
Definition: EclHysteresisTwoPhaseLaw.hpp:142
static void relativePermeabilities(Container &, const Params &, const FluidState &)
The relative permeability-saturation curves depending on absolute saturations.
Definition: EclHysteresisTwoPhaseLaw.hpp:123