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/EclipseState/EclipseState.hpp>
35 #endif
36 
37 #include <cassert>
38 #include <cmath>
39 #include <memory>
40 
42 
43 namespace Opm {
50 template <class EffLawT>
52 {
53  using EffLawParams = typename EffLawT::Params;
54  using Scalar = typename EffLawParams::Traits::Scalar;
55 
56 public:
57  using Traits = typename EffLawParams::Traits;
58 
60  {
61  // These are initialized to two (even though they represent saturations)
62  // to signify that the values are larger than physically possible, and force
63  // using the drainage curve before the first saturation update.
64  pcSwMdc_ = 2.0;
65  krnSwMdc_ = 2.0;
66  // krwSwMdc_ = 2.0;
67 
68  pcSwMic_ = -1.0;
69  initialImb_ = false;
70  oilWaterSystem_ = false;
71  pcmaxd_ = 0.0;
72  pcmaxi_ = 0.0;
73 
74  deltaSwImbKrn_ = 0.0;
75  // deltaSwImbKrw_ = 0.0;
76  }
77 
82  void finalize()
83  {
84  if (config().enableHysteresis()) {
85  if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().pcHysteresisModel() == 0) {
86  C_ = 1.0/(Sncri_ - Sncrd_ + 1.0e-12) - 1.0/(Snmaxd_ - Sncrd_);
87  curvatureCapPrs_ = config().curvatureCapPrs();
88  }
89 
90  updateDynamicParams_();
91  }
92 
94  }
95 
99  void setConfig(std::shared_ptr<EclHysteresisConfig> value)
100  { config_ = value; }
101 
106  { return *config_; }
107 
111  void setDrainageParams(const EffLawParams& value,
113  EclTwoPhaseSystemType twoPhaseSystem)
114  {
115  drainageParams_ = value;
116 
117  oilWaterSystem_ = (twoPhaseSystem == EclOilWaterSystem);
118 
119  if (!config().enableHysteresis())
120  return;
121 
122  if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().pcHysteresisModel() == 0) {
123  if (twoPhaseSystem == EclGasOilSystem) {
124  Sncrd_ = info.Sgcr+info.Swl;
125  Snmaxd_ = info.Sgu+info.Swl;
126  KrndMax_ = EffLawT::twoPhaseSatKrn(drainageParams(), 1.0-Snmaxd_);
127  }
128  else {
129  assert(twoPhaseSystem == EclOilWaterSystem);
130  Sncrd_ = info.Sowcr;
131  Snmaxd_ = 1.0 - info.Swl - info.Sgl;
132  KrndMax_ = EffLawT::twoPhaseSatKrn(drainageParams(), 1.0-Snmaxd_);
133  }
134  }
135 
136  // Additional Killough hysteresis model for pc
137  if (config().pcHysteresisModel() == 0) {
138  if (twoPhaseSystem == EclGasOilSystem) {
139  Swcrd_ = info.Sogcr;
140  pcmaxd_ = info.maxPcgo;
141  }
142  else {
143  assert(twoPhaseSystem == EclOilWaterSystem);
144  Swcrd_ = info.Swcr;
145  pcmaxd_ = -17.0; // At this point 'info.maxPcow' holds pre-swatinit value ...;
146  }
147  }
148  }
149 
153  const EffLawParams& drainageParams() const
154  { return drainageParams_; }
155 
156  EffLawParams& drainageParams()
157  { return drainageParams_; }
158 
162  void setImbibitionParams(const EffLawParams& value,
164  EclTwoPhaseSystemType twoPhaseSystem)
165  {
166  imbibitionParams_ = value;
167 
168  if (!config().enableHysteresis())
169  return;
170 
171  // Killough hysteresis model for nonw kr
172  if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().pcHysteresisModel() == 0) {
173  if (twoPhaseSystem == EclGasOilSystem) {
174  Sncri_ = info.Sgcr+info.Swl;
175  }
176  else {
177  assert(twoPhaseSystem == EclOilWaterSystem);
178  Sncri_ = info.Sowcr;
179  }
180  }
181 
182  // Killough hysteresis model for pc
183  if (config().pcHysteresisModel() == 0) {
184  if (twoPhaseSystem == EclGasOilSystem) {
185  Swcri_ = info.Sogcr;
186  Swmaxi_ = 1.0 - info.Sgl - info.Swl;
187  pcmaxi_ = info.maxPcgo;
188  }
189  else {
190  assert(twoPhaseSystem == EclOilWaterSystem);
191  Swcri_ = info.Swcr;
192  Swmaxi_ = info.Swu;
193  pcmaxi_ = info.maxPcow;
194  }
195  }
196  }
197 
201  const EffLawParams& imbibitionParams() const
202  { return imbibitionParams_; }
203 
204  EffLawParams& imbibitionParams()
205  { return imbibitionParams_; }
206 
211  Scalar pcSwMdc() const
212  { return pcSwMdc_; }
213 
214  Scalar pcSwMic() const
215  { return pcSwMic_; }
216 
220  bool initialImb() const
221  { return initialImb_; }
222 
228  void setKrwSwMdc(Scalar /* value */)
229  {}
230  // { krwSwMdc_ = value; };
231 
237  Scalar krwSwMdc() const
238  { return 2.0; }
239  // { return krwSwMdc_; };
240 
246  void setKrnSwMdc(Scalar value)
247  { krnSwMdc_ = value; }
248 
254  Scalar krnSwMdc() const
255  { return krnSwMdc_; }
256 
264  void setDeltaSwImbKrw(Scalar /* value */)
265  {}
266  // { deltaSwImbKrw_ = value; }
267 
275  Scalar deltaSwImbKrw() const
276  { return 0.0; }
277 // { return deltaSwImbKrw_; }
278 
286  void setDeltaSwImbKrn(Scalar value)
287  { deltaSwImbKrn_ = value; }
288 
296  Scalar deltaSwImbKrn() const
297  { return deltaSwImbKrn_; }
298 
299 
300  Scalar Swcri() const
301  { return Swcri_; }
302 
303  Scalar Swcrd() const
304  { return Swcrd_; }
305 
306  Scalar Swmaxi() const
307  { return Swmaxi_; }
308 
309  Scalar Sncri() const
310  { return Sncri_; }
311 
312  Scalar Sncrd() const
313  { return Sncrd_; }
314 
315  Scalar Sncrt() const
316  { return Sncrt_; }
317 
318  Scalar Snmaxd() const
319  { return Snmaxd_; }
320 
321  Scalar Snhy() const
322  { return 1.0 - krnSwMdc_; }
323 
324  Scalar krnWght() const
325  { return KrndHy_/KrndMax_; }
326 
327  Scalar pcWght() const // Aligning pci and pcd at Swir
328  {
329  if (pcmaxd_ < 0.0)
330  return EffLawT::twoPhaseSatPcnw(drainageParams(), 0.0)/(pcmaxi_+1e-6);
331  else
332  return pcmaxd_/(pcmaxi_+1e-6);
333  }
334 
335  Scalar curvatureCapPrs() const
336  { return curvatureCapPrs_;}
337 
338 
345  void update(Scalar pcSw, Scalar /* krwSw */, Scalar krnSw)
346  {
347  bool updateParams = false;
348 
349  if (config().pcHysteresisModel() == 0 && pcSw < pcSwMdc_) {
350  if (pcSwMdc_ == 2.0 && pcSw+1.0e-6 < Swcrd_ && oilWaterSystem_) {
351  initialImb_ = true;
352  }
353  pcSwMdc_ = pcSw;
354  updateParams = true;
355  }
356 
357  if (initialImb_ && pcSw > pcSwMic_) {
358  pcSwMic_ = pcSw;
359  updateParams = true;
360  }
361 
362 /*
363  // This is quite hacky: Eclipse says that it only uses relperm hysteresis for the
364  // wetting phase (indicated by '0' for the second item of the EHYSTER keyword),
365  // even though this makes about zero sense: one would expect that hysteresis can
366  // be limited to the oil phase, but the oil phase is the wetting phase of the
367  // gas-oil twophase system whilst it is non-wetting for water-oil.
368  if (krwSw < krwSwMdc_)
369  {
370  krwSwMdc_ = krwSw;
371  updateParams = true;
372  }
373 */
374 
375  if (krnSw < krnSwMdc_) {
376  krnSwMdc_ = krnSw;
377  KrndHy_ = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_);
378  updateParams = true;
379  }
380 
381  if (updateParams)
382  updateDynamicParams_();
383  }
384 
385 private:
386  void updateDynamicParams_()
387  {
388  // HACK: Eclipse seems to disable the wetting-phase relperm even though this is
389  // quite pointless from the physical POV. (see comment above)
390 /*
391  // calculate the saturation deltas for the relative permeabilities
392  Scalar krwMdcDrainage = EffLawT::twoPhaseSatKrw(drainageParams(), krwSwMdc_);
393  Scalar SwKrwMdcImbibition = EffLawT::twoPhaseSatKrwInv(imbibitionParams(), krwMdcDrainage);
394  deltaSwImbKrw_ = SwKrwMdcImbibition - krwSwMdc_;
395 */
396  if (config().krHysteresisModel() == 0 || config().krHysteresisModel() == 1) {
397  Scalar krnMdcDrainage = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_);
398  Scalar SwKrnMdcImbibition = EffLawT::twoPhaseSatKrnInv(imbibitionParams(), krnMdcDrainage);
399  deltaSwImbKrn_ = SwKrnMdcImbibition - krnSwMdc_;
400  assert(std::abs(EffLawT::twoPhaseSatKrn(imbibitionParams(), krnSwMdc_ + deltaSwImbKrn_)
401  - EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_)) < 1e-8);
402  }
403 
404  // Scalar pcMdcDrainage = EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_);
405  // Scalar SwPcMdcImbibition = EffLawT::twoPhaseSatPcnwInv(imbibitionParams(), pcMdcDrainage);
406  // deltaSwImbPc_ = SwPcMdcImbibition - pcSwMdc_;
407 
408 // assert(std::abs(EffLawT::twoPhaseSatPcnw(imbibitionParams(), pcSwMdc_ + deltaSwImbPc_)
409 // - EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_)) < 1e-8);
410 // assert(std::abs(EffLawT::twoPhaseSatKrn(imbibitionParams(), krnSwMdc_ + deltaSwImbKrn_)
411 // - EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_)) < 1e-8);
412 // assert(std::abs(EffLawT::twoPhaseSatKrw(imbibitionParams(), krwSwMdc_ + deltaSwImbKrw_)
413 // - EffLawT::twoPhaseSatKrw(drainageParams(), krwSwMdc_)) < 1e-8);
414 
415  if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().pcHysteresisModel() == 0) {
416  Scalar Snhy = 1.0 - krnSwMdc_;
417  if (Snhy > Sncrd_)
418  Sncrt_ = Sncrd_ + (Snhy - Sncrd_)/((1.0+config().modParamTrapped()*(Snmaxd_-Snhy)) + C_*(Snhy - Sncrd_));
419  else
420  {
421  Sncrt_ = Sncrd_;
422  }
423  }
424  }
425 
426  std::shared_ptr<EclHysteresisConfig> config_;
427  EffLawParams imbibitionParams_;
428  EffLawParams drainageParams_;
429 
430  // largest wettinging phase saturation which is on the main-drainage curve. These are
431  // three different values because the sourounding code can choose to use different
432  // definitions for the saturations for different quantities
433  //Scalar krwSwMdc_;
434  Scalar krnSwMdc_;
435  Scalar pcSwMdc_;
436 
437  // largest wettinging phase saturation along main imbibition curve
438  Scalar pcSwMic_;
439  // Initial process is imbibition (for initial saturations at or below critical drainage saturation)
440  bool initialImb_;
441 
442  bool oilWaterSystem_;
443 
444  // offsets added to wetting phase saturation uf using the imbibition curves need to
445  // be used to calculate the wetting phase relperm, the non-wetting phase relperm and
446  // the capillary pressure
447  //Scalar deltaSwImbKrw_;
448  Scalar deltaSwImbKrn_;
449  //Scalar deltaSwImbPc_;
450 
451  // trapped non-wetting phase saturation
452  Scalar Sncrt_;
453 
454  // the following uses the conventions of the Eclipse technical description:
455  //
456  // Sncrd_: critical non-wetting phase saturation for the drainage curve
457  // Sncri_: critical non-wetting phase saturation for the imbibition curve
458  // Swcri_: critical wetting phase saturation for the imbibition curve
459  // Swcrd_: critical wetting phase saturation for the drainage curve
460  // Swmaxi_; maximum wetting phase saturation for the imbibition curve
461  // Snmaxd_: non-wetting phase saturation where the non-wetting relperm reaches its
462  // maximum on the drainage curve
463  // C_: factor required to calculate the trapped non-wetting phase saturation using
464  // the Killough approach
465  Scalar Sncrd_;
466  Scalar Sncri_;
467  Scalar Swcri_;
468  Scalar Swcrd_;
469  Scalar Swmaxi_;
470  Scalar Snmaxd_;
471  Scalar C_;
472 
473  Scalar KrndMax_; // Krn_drain(Snmaxd_)
474  Scalar KrndHy_; // Krn_drain(1-krnSwMdc_)
475 
476  Scalar pcmaxd_; // max pc for drain
477  Scalar pcmaxi_; // max pc for imb
478 
479  Scalar curvatureCapPrs_; // curvature parameter used for capillary pressure hysteresis
480 };
481 
482 } // namespace Opm
483 
484 #endif
EclTwoPhaseSystemType
Specified which fluids are involved in a given twophase material law for endpoint scaling.
Definition: EclEpsConfig.hpp:43
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:41
int pcHysteresisModel() const
Return the type of the hysteresis model which is used for capillary pressure.
Definition: EclHysteresisConfig.hpp:77
int krHysteresisModel() const
Return the type of the hysteresis model which is used for relative permeability.
Definition: EclHysteresisConfig.hpp:103
double curvatureCapPrs() const
Curvature parameter used for capillary pressure hysteresis.
Definition: EclHysteresisConfig.hpp:119
bool enableHysteresis() const
Returns whether hysteresis is enabled.
Definition: EclHysteresisConfig.hpp:59
A default implementation of the parameters for the material law which implements the ECL relative per...
Definition: EclHysteresisTwoPhaseLawParams.hpp:52
Scalar pcSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition: EclHysteresisTwoPhaseLawParams.hpp:211
bool initialImb() const
Status of initial process.
Definition: EclHysteresisTwoPhaseLawParams.hpp:220
const EffLawParams & drainageParams() const
Returns the parameters used for the drainage curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:153
void setDeltaSwImbKrn(Scalar value)
Sets the saturation value which must be added if krn is calculated using the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:286
void setConfig(std::shared_ptr< EclHysteresisConfig > value)
Set the endpoint scaling configuration object.
Definition: EclHysteresisTwoPhaseLawParams.hpp:99
Scalar krnSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition: EclHysteresisTwoPhaseLawParams.hpp:254
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:246
void setDeltaSwImbKrw(Scalar)
Sets the saturation value which must be added if krw is calculated using the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:264
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:228
void finalize()
Calculate all dependent quantities once the independent quantities of the parameter object have been ...
Definition: EclHysteresisTwoPhaseLawParams.hpp:82
void setDrainageParams(const EffLawParams &value, const EclEpsScalingPointsInfo< Scalar > &info, EclTwoPhaseSystemType twoPhaseSystem)
Sets the parameters used for the drainage curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:111
Scalar deltaSwImbKrn() const
Returns the saturation value which must be added if krn is calculated using the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:296
void setImbibitionParams(const EffLawParams &value, const EclEpsScalingPointsInfo< Scalar > &info, EclTwoPhaseSystemType twoPhaseSystem)
Sets the parameters used for the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:162
Scalar deltaSwImbKrw() const
Returns the saturation value which must be added if krw is calculated using the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:275
const EclHysteresisConfig & config() const
Returns the endpoint scaling configuration object.
Definition: EclHysteresisTwoPhaseLawParams.hpp:105
Scalar krwSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition: EclHysteresisTwoPhaseLawParams.hpp:237
const EffLawParams & imbibitionParams() const
Returns the parameters used for the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:201
void update(Scalar pcSw, Scalar, Scalar krnSw)
Notify the hysteresis law that a given wetting-phase saturation has been seen.
Definition: EclHysteresisTwoPhaseLawParams.hpp:345
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:60