My Project
EclHysteresisTwoPhaseLawParams.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_PARAMS_HPP
28 #define OPM_ECL_HYSTERESIS_TWO_PHASE_LAW_PARAMS_HPP
29 
30 #include "EclHysteresisConfig.hpp"
31 #include "EclEpsScalingPoints.hpp"
32 
33 #if HAVE_ECL_INPUT
34 #include <opm/input/eclipse/Deck/Deck.hpp>
35 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
36 #endif
37 
38 #include <string>
39 #include <memory>
40 #include <cassert>
41 #include <algorithm>
42 
44 
45 namespace Opm {
52 template <class EffLawT>
54 {
55  typedef typename EffLawT::Params EffLawParams;
56  typedef typename EffLawParams::Traits::Scalar Scalar;
57 
58 public:
59  typedef typename EffLawParams::Traits Traits;
60 
62  {
63  // These are initialized to two (even though they represent saturations)
64  // to signify that the values are larger than physically possible, and force
65  // using the drainage curve before the first saturation update.
66  pcSwMdc_ = 2.0;
67  krnSwMdc_ = 2.0;
68  // krwSwMdc_ = 2.0;
69 
70  deltaSwImbKrn_ = 0.0;
71  // deltaSwImbKrw_ = 0.0;
72  }
73 
78  void finalize()
79  {
80  if (config().enableHysteresis()) {
81  //C_ = 1.0/(Sncri_ - Sncrd_) + 1.0/(Snmaxd_ - Sncrd_);
82 
83  updateDynamicParams_();
84  }
85 
87  }
88 
92  void setConfig(std::shared_ptr<EclHysteresisConfig> value)
93  { config_ = value; }
94 
98  const EclHysteresisConfig& config() const
99  { return *config_; }
100 
104  void setDrainageParams(const EffLawParams& value)
105  {
106  drainageParams_ = value;
107  }
108 
112  const EffLawParams& drainageParams() const
113  { return drainageParams_; }
114 
115  EffLawParams& drainageParams()
116  { return drainageParams_; }
117 
121  void setImbibitionParams(const EffLawParams& value,
122  const EclEpsScalingPointsInfo<Scalar>& /* info */,
123  EclTwoPhaseSystemType /* twoPhaseSystem */)
124  {
125  imbibitionParams_ = value;
126 
127 /*
128  if (twoPhaseSystem == EclGasOilSystem) {
129  Sncri_ = info.Sgcr;
130  }
131  else {
132  assert(twoPhaseSystem == EclOilWaterSystem);
133  Sncri_ = info.Sowcr;
134  }
135 */
136  }
137 
141  const EffLawParams& imbibitionParams() const
142  { return imbibitionParams_; }
143 
144  EffLawParams& imbibitionParams()
145  { return imbibitionParams_; }
146 
151  void setPcSwMdc(Scalar value)
152  { pcSwMdc_ = value; }
153 
158  Scalar pcSwMdc() const
159  { return pcSwMdc_; }
160 
166  void setKrwSwMdc(Scalar /* value */)
167  {}
168  // { krwSwMdc_ = value; };
169 
175  Scalar krwSwMdc() const
176  { return 2.0; }
177  // { return krwSwMdc_; };
178 
184  void setKrnSwMdc(Scalar value)
185  { krnSwMdc_ = value; }
186 
192  Scalar krnSwMdc() const
193  { return krnSwMdc_; }
194 
202  void setDeltaSwImbKrw(Scalar /* value */)
203  {}
204  // { deltaSwImbKrw_ = value; }
205 
213  Scalar deltaSwImbKrw() const
214  { return 0.0; }
215 // { return deltaSwImbKrw_; }
216 
224  void setDeltaSwImbKrn(Scalar value)
225  { deltaSwImbKrn_ = value; }
226 
234  Scalar deltaSwImbKrn() const
235  { return deltaSwImbKrn_; }
236 
243  void update(Scalar pcSw, Scalar /* krwSw */, Scalar krnSw)
244  {
245  bool updateParams = false;
246  if (pcSw < pcSwMdc_) {
247  pcSwMdc_ = pcSw;
248  updateParams = true;
249  }
250 
251 /*
252  // This is quite hacky: Eclipse says that it only uses relperm hysteresis for the
253  // wetting phase (indicated by '0' for the second item of the EHYSTER keyword),
254  // even though this makes about zero sense: one would expect that hysteresis can
255  // be limited to the oil phase, but the oil phase is the wetting phase of the
256  // gas-oil twophase system whilst it is non-wetting for water-oil.
257  if (krwSw < krwSwMdc_)
258  {
259  krwSwMdc_ = krwSw;
260  updateParams = true;
261  }
262 */
263 
264  if (krnSw < krnSwMdc_) {
265  krnSwMdc_ = krnSw;
266  updateParams = true;
267  }
268 
269  if (updateParams)
270  updateDynamicParams_();
271  }
272 
273 private:
274  void updateDynamicParams_()
275  {
276  // HACK: Eclipse seems to disable the wetting-phase relperm even though this is
277  // quite pointless from the physical POV. (see comment above)
278 /*
279  // calculate the saturation deltas for the relative permeabilities
280  Scalar krwMdcDrainage = EffLawT::twoPhaseSatKrw(drainageParams(), krwSwMdc_);
281  Scalar SwKrwMdcImbibition = EffLawT::twoPhaseSatKrwInv(imbibitionParams(), krwMdcDrainage);
282  deltaSwImbKrw_ = SwKrwMdcImbibition - krwSwMdc_;
283 */
284 
285  Scalar krnMdcDrainage = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_);
286  Scalar SwKrnMdcImbibition = EffLawT::twoPhaseSatKrnInv(imbibitionParams(), krnMdcDrainage);
287  deltaSwImbKrn_ = SwKrnMdcImbibition - krnSwMdc_;
288 
289  // Scalar pcMdcDrainage = EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_);
290  // Scalar SwPcMdcImbibition = EffLawT::twoPhaseSatPcnwInv(imbibitionParams(), pcMdcDrainage);
291  // deltaSwImbPc_ = SwPcMdcImbibition - pcSwMdc_;
292 
293 // assert(std::abs(EffLawT::twoPhaseSatPcnw(imbibitionParams(), pcSwMdc_ + deltaSwImbPc_)
294 // - EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_)) < 1e-8);
295  assert(std::abs(EffLawT::twoPhaseSatKrn(imbibitionParams(), krnSwMdc_ + deltaSwImbKrn_)
296  - EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_)) < 1e-8);
297 // assert(std::abs(EffLawT::twoPhaseSatKrw(imbibitionParams(), krwSwMdc_ + deltaSwImbKrw_)
298 // - EffLawT::twoPhaseSatKrw(drainageParams(), krwSwMdc_)) < 1e-8);
299 
300 #if 0
301  Scalar Snhy = 1.0 - SwMdc_;
302 
303  Sncrt_ = Sncrd_ + (Snhy - Sncrd_)/(1 + C_*(Snhy - Sncrd_));
304 #endif
305  }
306 
307  std::shared_ptr<EclHysteresisConfig> config_;
308  EffLawParams imbibitionParams_;
309  EffLawParams drainageParams_;
310 
311  // largest wettinging phase saturation which is on the main-drainage curve. These are
312  // three different values because the sourounding code can choose to use different
313  // definitions for the saturations for different quantities
314  //Scalar krwSwMdc_;
315  Scalar krnSwMdc_;
316  Scalar pcSwMdc_;
317 
318  // offsets added to wetting phase saturation uf using the imbibition curves need to
319  // be used to calculate the wetting phase relperm, the non-wetting phase relperm and
320  // the capillary pressure
321  //Scalar deltaSwImbKrw_;
322  Scalar deltaSwImbKrn_;
323  //Scalar deltaSwImbPc_;
324 
325  // trapped non-wetting phase saturation
326  //Scalar Sncrt_;
327 
328  // the following uses the conventions of the Eclipse technical description:
329  //
330  // Sncrd_: critical non-wetting phase saturation for the drainage curve
331  // Sncri_: critical non-wetting phase saturation for the imbibition curve
332  // Snmaxd_: non-wetting phase saturation where the non-wetting relperm reaches its
333  // maximum on the drainage curve
334  // C_: factor required to calculate the trapped non-wetting phase saturation using
335  // the Killough approach
336  //Scalar Sncrd_;
337  //Scalar Sncri_;
338  //Scalar Snmaxd_;
339  //Scalar C_;
340 };
341 
342 } // namespace Opm
343 
344 #endif
EclTwoPhaseSystemType
Specified which fluids are involved in a given twophase material law for endpoint scaling.
Definition: EclEpsConfig.hpp:49
Specifies the configuration used by the ECL kr/pC hysteresis code.
Default implementation for asserting finalization of parameter objects.
Specifies the configuration used by the ECL kr/pC hysteresis code.
Definition: EclHysteresisConfig.hpp:48
A default implementation of the parameters for the material law which implements the ECL relative per...
Definition: EclHysteresisTwoPhaseLawParams.hpp:54
Scalar pcSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition: EclHysteresisTwoPhaseLawParams.hpp:158
const EffLawParams & drainageParams() const
Returns the parameters used for the drainage curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:112
void setDeltaSwImbKrn(Scalar value)
Sets the saturation value which must be added if krn is calculated using the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:224
void setConfig(std::shared_ptr< EclHysteresisConfig > value)
Set the endpoint scaling configuration object.
Definition: EclHysteresisTwoPhaseLawParams.hpp:92
void setImbibitionParams(const EffLawParams &value, const EclEpsScalingPointsInfo< Scalar > &, EclTwoPhaseSystemType)
Sets the parameters used for the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:121
Scalar krnSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition: EclHysteresisTwoPhaseLawParams.hpp:192
void setKrnSwMdc(Scalar value)
Set the saturation of the wetting phase where the last switch from the main drainage curve (MDC) to i...
Definition: EclHysteresisTwoPhaseLawParams.hpp:184
void setDeltaSwImbKrw(Scalar)
Sets the saturation value which must be added if krw is calculated using the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:202
void setDrainageParams(const EffLawParams &value)
Sets the parameters used for the drainage curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:104
void setPcSwMdc(Scalar value)
Set the saturation of the wetting phase where the last switch from the main drainage curve (MDC) to i...
Definition: EclHysteresisTwoPhaseLawParams.hpp:151
void setKrwSwMdc(Scalar)
Set the saturation of the wetting phase where the last switch from the main drainage curve (MDC) to i...
Definition: EclHysteresisTwoPhaseLawParams.hpp:166
void finalize()
Calculate all dependent quantities once the independent quantities of the parameter object have been ...
Definition: EclHysteresisTwoPhaseLawParams.hpp:78
Scalar deltaSwImbKrn() const
Returns the saturation value which must be added if krn is calculated using the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:234
Scalar deltaSwImbKrw() const
Returns the saturation value which must be added if krw is calculated using the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:213
const EclHysteresisConfig & config() const
Returns the endpoint scaling configuration object.
Definition: EclHysteresisTwoPhaseLawParams.hpp:98
Scalar krwSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition: EclHysteresisTwoPhaseLawParams.hpp:175
const EffLawParams & imbibitionParams() const
Returns the parameters used for the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:141
void update(Scalar pcSw, Scalar, Scalar krnSw)
Notify the hysteresis law that a given wetting-phase saturation has been seen.
Definition: EclHysteresisTwoPhaseLawParams.hpp:243
Default implementation for asserting finalization of parameter objects.
Definition: EnsureFinalized.hpp:47
void finalize()
Mark the object as finalized.
Definition: EnsureFinalized.hpp:75
This structure represents all values which can be possibly used as scaling points in the endpoint sca...
Definition: EclEpsScalingPoints.hpp:62