My Project
ParkerLenhard.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_PARKER_LENHARD_HPP
28 #define OPM_PARKER_LENHARD_HPP
29 
30 #include "ParkerLenhardParams.hpp"
31 
33 
34 #include <algorithm>
35 #include <cassert>
36 
37 namespace Opm {
38 
46 template <class ScalarT>
48 {
49 public:
50  typedef ScalarT Scalar;
51 
58  PLScanningCurve(Scalar Swr)
59  {
60  loopNum_ = 0;
61  prev_ = new PLScanningCurve(NULL, // prev
62  this, // next
63  -1, // loop number
64  Swr, // Sw
65  1e12, // pcnw
66  Swr, // SwMic
67  Swr); // SwMdc
68  next_ = NULL;
69 
70  Sw_ = 1.0;
71  pcnw_ = 0.0;
72  SwMic_ = 1.0;
73  SwMdc_ = 1.0;
74  }
75 
76 protected:
78  PLScanningCurve* nextSC,
79  int loopN,
80  Scalar SwReversal,
81  Scalar pcnwReversal,
82  Scalar SwMiCurve,
83  Scalar SwMdCurve)
84  {
85  prev_ = prevSC;
86  next_ = nextSC;
87  loopNum_ = loopN;
88  Sw_ = SwReversal;
89  pcnw_ = pcnwReversal;
90  SwMic_ = SwMiCurve;
91  SwMdc_ = SwMdCurve;
92  }
93 
94 public:
101  {
102  if (loopNum_ == 0)
103  delete prev_;
104  if (loopNum_ >= 0)
105  delete next_;
106  }
107 
113  { return prev_; }
114 
120  { return next_; }
121 
130  void setNext(Scalar SwReversal,
131  Scalar pcnwReversal,
132  Scalar SwMiCurve,
133  Scalar SwMdCurve)
134  {
135  // if next_ is NULL, delete does nothing, so
136  // this is valid!!
137  delete next_;
138 
139  next_ = new PLScanningCurve(this, // prev
140  NULL, // next
141  loopNum() + 1,
142  SwReversal,
143  pcnwReversal,
144  SwMiCurve,
145  SwMdCurve);
146  }
147 
154  bool isValidAt_Sw(Scalar SwReversal)
155  {
156  if (isImbib())
157  // for inbibition the given saturation
158  // must be between the start of the
159  // current imbibition and the the start
160  // of the last drainage
161  return this->Sw() < SwReversal && SwReversal < prev_->Sw();
162  else
163  // for drainage the given saturation
164  // must be between the start of the
165  // last imbibition and the start
166  // of the current drainage
167  return prev_->Sw() < SwReversal && SwReversal < this->Sw();
168  }
169 
174  bool isImbib()
175  { return loopNum()%2 == 1; }
176 
181  bool isDrain()
182  { return !isImbib(); }
183 
189  int loopNum()
190  { return loopNum_; }
191 
196  Scalar Sw() const
197  { return Sw_; }
198 
202  Scalar pcnw() const
203  { return pcnw_; }
204 
209  Scalar SwMic()
210  { return SwMic_; }
211 
216  Scalar SwMdc()
217  { return SwMdc_; }
218 
219 private:
220  PLScanningCurve* prev_;
221  PLScanningCurve* next_;
222 
223  int loopNum_;
224 
225  Scalar Sw_;
226  Scalar pcnw_;
227 
228  Scalar SwMdc_;
229  Scalar SwMic_;
230 };
231 
238 template <class TraitsT, class ParamsT = ParkerLenhardParams<TraitsT> >
239 class ParkerLenhard : public TraitsT
240 {
241 public:
242  typedef TraitsT Traits;
243  typedef ParamsT Params;
244  typedef typename Traits::Scalar Scalar;
245 
247  static const int numPhases = Traits::numPhases;
248  static_assert(numPhases == 2,
249  "The Parker-Lenhard capillary pressure law only "
250  "applies to the case of two fluid phases");
251 
254  static const bool implementsTwoPhaseApi = true;
255 
258  static const bool implementsTwoPhaseSatApi = true;
259 
262  static const bool isSaturationDependent = true;
263 
266  static const bool isPressureDependent = false;
267 
270  static const bool isTemperatureDependent = false;
271 
274  static const bool isCompositionDependent = false;
275 
276  static_assert(Traits::numPhases == 2,
277  "The number of fluid phases must be two if you want to use "
278  "this material law!");
279 
280 private:
281  typedef typename ParamsT::VanGenuchten VanGenuchten;
283 
284 public:
289  static void reset(Params& params)
290  {
291  delete params.mdc(); // this will work even if mdc_ == NULL!
292  params.setMdc(new ScanningCurve(params.SwrPc()));
293  params.setCsc(params.mdc());
294  params.setPisc(NULL);
295  params.setCurrentSnr(0.0);
296  }
297 
302  template <class FluidState>
303  static void update(Params& params, const FluidState& fs)
304  {
305  Scalar Sw = scalarValue(fs.saturation(Traits::wettingPhaseIdx));
306 
307  if (Sw > 1 - 1e-5) {
308  // if the absolute saturation is almost 1,
309  // it means that we're back to the beginning
310  reset(params);
311  return;
312  }
313 
314  // find the loop number which corrosponds to the
315  // given effective saturation
316  ScanningCurve* curve = findScanningCurve_(params, Sw);
317 
318  // calculate the apparent saturation on the MIC and MDC
319  // which yield the same capillary pressure as the
320  // Sw at the current scanning curve
321  Scalar pc = pcnw<FluidState, Scalar>(params, fs);
322  Scalar Sw_mic = VanGenuchten::twoPhaseSatSw(params.micParams(), pc);
323  Scalar Sw_mdc = VanGenuchten::twoPhaseSatSw(params.mdcParams(), pc);
324 
325  curve->setNext(Sw, pc, Sw_mic, Sw_mdc);
326  if (!curve->next())
327  return;
328 
329  params.setCsc(curve);
330 
331  // if we're back on the MDC, we also have a new PISC!
332  if (params.csc() == params.mdc()) {
333  params.setPisc(params.mdc()->next());
334  params.setCurrentSnr(computeCurrentSnr_(params, Sw));
335  }
336  }
337 
342  template <class Container, class FluidState>
343  static void capillaryPressures(Container& values, const Params& params, const FluidState& fs)
344  {
345  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
346 
347  values[Traits::wettingPhaseIdx] = 0.0; // reference phase
348  values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
349  }
350 
355  template <class Container, class FluidState>
356  static void saturations(Container& /*values*/, const Params& /*params*/, const FluidState& /*fs*/)
357  { throw std::logic_error("Not implemented: ParkerLenhard::saturations()"); }
358 
363  template <class Container, class FluidState>
364  static void relativePermeabilities(Container& values, const Params& params, const FluidState& fs)
365  {
366  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
367 
368  values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
369  values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
370  }
371 
376  template <class FluidState, class Evaluation = typename FluidState::Scalar>
377  static Evaluation pcnw(const Params& params, const FluidState& fs)
378  {
379  const Evaluation& Sw =
380  decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
381 
382  return twoPhaseSatPcnw(params, Sw);
383  }
384 
385  template <class Evaluation>
386  static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
387  {
388  // calculate the current apparent saturation
389  ScanningCurve* sc = findScanningCurve_(params, scalarValue(Sw));
390 
391  // calculate the apparant saturation
392  const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
393 
394  // if the apparent saturation exceeds the 'legal' limits,
395  // we also the underlying material law decide what to do.
396  if (Sw_app > 1) {
397  return 0.0; // VanGenuchten::pcnw(params.mdcParams(), Sw_app);
398  }
399 
400  // put the apparent saturation into the main imbibition or
401  // drainage curve
402  Scalar SwAppCurSC = absoluteToApparentSw_(params, sc->Sw());
403  Scalar SwAppPrevSC = absoluteToApparentSw_(params, sc->prev()->Sw());
404  const Evaluation& pos = (Sw_app - SwAppCurSC)/(SwAppPrevSC - SwAppCurSC);
405  if (sc->isImbib()) {
406  const Evaluation& SwMic =
407  pos * (sc->prev()->SwMic() - sc->SwMic()) + sc->SwMic();
408 
409  return VanGenuchten::twoPhaseSatPcnw(params.micParams(), SwMic);
410  }
411  else { // sc->isDrain()
412  const Evaluation& SwMdc =
413  pos*(sc->prev()->SwMdc() - sc->SwMdc()) + sc->SwMdc();
414 
415  return VanGenuchten::twoPhaseSatPcnw(params.mdcParams(), SwMdc);
416  }
417  }
418 
423  template <class FluidState, class Evaluation = typename FluidState::Scalar>
424  static Evaluation Sw(const Params& /*params*/, const FluidState& /*fs*/)
425  { throw std::logic_error("Not implemented: ParkerLenhard::Sw()"); }
426 
427  template <class Evaluation>
428  static Evaluation twoPhaseSatSw(const Params& /*params*/, const Evaluation& /*pc*/)
429  { throw std::logic_error("Not implemented: ParkerLenhard::twoPhaseSatSw()"); }
430 
435  template <class FluidState, class Evaluation = typename FluidState::Scalar>
436  static Evaluation Sn(const Params& params, const FluidState& fs)
437  { return 1 - Sw<FluidState, Evaluation>(params, fs); }
438 
439  template <class Evaluation>
440  static Evaluation twoPhaseSatSn(const Params& /*params*/, const Evaluation& /*pc*/)
441  { throw std::logic_error("Not implemented: ParkerLenhard::twoPhaseSatSn()"); }
442 
447  template <class FluidState, class Evaluation = typename FluidState::Scalar>
448  static Evaluation krw(const Params& params, const FluidState& fs)
449  {
450  const Evaluation& Sw =
451  decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
452 
453  return twoPhaseSatKrw(params, Sw);
454  }
455 
456  template <class Evaluation>
457  static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
458  {
459  // for the effective permeability we only use Land's law and
460  // the relative permeability of the main drainage curve.
461  const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
462  return VanGenuchten::twoPhaseSatKrw(params.mdcParams(), Sw_app);
463  }
464 
469  template <class FluidState, class Evaluation = typename FluidState::Scalar>
470  static Evaluation krn(const Params& params, const FluidState& fs)
471  {
472  const Evaluation& Sw =
473  decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
474 
475  return twoPhaseSatKrn(params, Sw);
476  }
477 
478  template <class Evaluation>
479  static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
480  {
481  // for the effective permeability we only use Land's law and
482  // the relative permeability of the main drainage curve.
483  const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
484  return VanGenuchten::twoPhaseSatKrn(params.mdcParams(), Sw_app);
485  }
486 
490  template <class Evaluation>
491  static Evaluation absoluteToApparentSw_(const Params& params, const Evaluation& Sw)
492  {
493  return effectiveToApparentSw_(params, absoluteToEffectiveSw_(params, Sw));
494  }
495 
496 private:
506  template <class Evaluation>
507  static Evaluation absoluteToEffectiveSw_(const Params& params, const Evaluation& Sw)
508  { return (Sw - params.SwrPc())/(1 - params.SwrPc()); }
509 
519  template <class Evaluation>
520  static Evaluation effectiveToAbsoluteSw_(const Params& params, const Evaluation& Swe)
521  { return Swe*(1 - params.SwrPc()) + params.SwrPc(); }
522 
523  // return the effctive residual non-wetting saturation, given an
524  // effective wetting saturation
525  template <class Evaluation>
526  static Evaluation computeCurrentSnr_(const Params& params, const Evaluation& Sw)
527  {
528  // regularize
529  if (Sw > 1 - params.Snr())
530  return 0.0;
531  if (Sw < params.SwrPc())
532  return params.Snr();
533 
534  if (params.Snr() == 0.0)
535  return 0.0;
536 
537  // use Land's law
538  Scalar R = 1.0/params.Snr() - 1;
539  const Evaluation& curSnr = (1 - Sw)/(1 + R*(1 - Sw));
540 
541  // the current effective residual non-wetting saturation must
542  // be smaller than the residual non-wetting saturation
543  assert(curSnr <= params.Snr());
544 
545  return curSnr;
546  }
547 
548  // returns the trapped effective non-wetting saturation for a
549  // given wetting phase saturation
550  template <class Evaluation>
551  static Evaluation trappedEffectiveSn_(const Params& params, const Evaluation& Sw)
552  {
553  const Evaluation& Swe = absoluteToEffectiveSw_(params, Sw);
554  Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
555 
556  Scalar Snre = absoluteToEffectiveSw_(params, params.currentSnr());
557  return Snre*(Swe - SwePisc) / (1 - Snre - SwePisc);
558  }
559 
560  // returns the apparent saturation of the wetting phase depending
561  // on the effective saturation
562  template <class Evaluation>
563  static Evaluation effectiveToApparentSw_(const Params& params, const Evaluation& Swe)
564  {
565  if (params.pisc() == NULL ||
566  Swe <= absoluteToEffectiveSw_(params, params.pisc()->Sw()))
567  {
568  // we are on the main drainage curve, i.e. no non-wetting
569  // fluid is trapped -> apparent saturation == effective
570  // saturation
571  return Swe;
572  }
573 
574  // we are on a imbibition or drainage curve which is not the
575  // main drainage curve -> apparent saturation == effective
576  // saturation + trapped effective saturation
577  return Swe + trappedEffectiveSn_(params, Swe);
578  }
579 
580  // Returns the effective saturation to a given apparent one
581  template <class Evaluation>
582  static Evaluation apparentToEffectiveSw_(const Params& params, const Evaluation& Swapp)
583  {
584  Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
585  if (params.pisc() == NULL || Swapp <= SwePisc) {
586  // we are on the main drainage curve, i.e.
587  // no non-wetting fluid is trapped
588  // -> apparent saturation == effective saturation
589  return Swapp;
590  }
591 
592  Scalar Snre = absoluteToEffectiveSw_(params.currentSnr());
593  return
594  (Swapp*(1 - Snre - SwePisc) + Snre*SwePisc)
595  /(1 - SwePisc);
596  }
597 
598 
599  // find the loop on which the an effective
600  // saturation has to be
601  static ScanningCurve* findScanningCurve_(const Params& params, Scalar Sw)
602  {
603  if (params.pisc() == NULL || Sw <= params.pisc()->Sw()) {
604  // we don't have a PISC yet, or the effective
605  // saturation is smaller than the saturation where the
606  // PISC begins. In this case are on the MDC
607  return params.mdc();
608  }
609 
610  // if we have a primary imbibition curve, and our current
611  // effective saturation is higher than the beginning of
612  // the secondary drainage curve. this means we are on the
613  // PISC again.
614  if (params.pisc()->next() == NULL ||
615  params.pisc()->next()->Sw() < Sw)
616  {
617  return params.pisc();
618  }
619 
620  ScanningCurve* curve = params.csc()->next();
621  while (true) {
622  assert(curve != params.mdc()->prev());
623  if (curve->isValidAt_Sw(Sw)) {
624  return curve;
625  }
626  curve = curve->prev();
627  }
628  }
629 };
630 
631 } // namespace Opm
632 
633 #endif // PARKER_LENHARD_HPP
Default parameter class for the Parker-Lenhard hysteresis model.
Implementation of the van Genuchten capillary pressure - saturation relation.
Represents a scanning curve in the Parker-Lenhard hysteresis model.
Definition: ParkerLenhard.hpp:48
PLScanningCurve * prev() const
Return the previous scanning curve, i.e.
Definition: ParkerLenhard.hpp:112
PLScanningCurve * next() const
Return the next scanning curve, i.e.
Definition: ParkerLenhard.hpp:119
~PLScanningCurve()
Destructor.
Definition: ParkerLenhard.hpp:100
PLScanningCurve(Scalar Swr)
Constructs main imbibition curve.
Definition: ParkerLenhard.hpp:58
bool isDrain()
Returns true iff the scanning curve is a drainage curve.
Definition: ParkerLenhard.hpp:181
void setNext(Scalar SwReversal, Scalar pcnwReversal, Scalar SwMiCurve, Scalar SwMdCurve)
Set the next scanning curve.
Definition: ParkerLenhard.hpp:130
bool isImbib()
Returns true iff the scanning curve is a imbibition curve.
Definition: ParkerLenhard.hpp:174
Scalar Sw() const
Absolute wetting-phase saturation at the scanning curve's reversal point.
Definition: ParkerLenhard.hpp:196
bool isValidAt_Sw(Scalar SwReversal)
Returns true iff the given effective saturation Swei is within the scope of the curve,...
Definition: ParkerLenhard.hpp:154
Scalar SwMdc()
Apparent saturation of the last reversal point on the pressure MDC.
Definition: ParkerLenhard.hpp:216
Scalar pcnw() const
Capillary pressure at the last reversal point.
Definition: ParkerLenhard.hpp:202
int loopNum()
The loop number of the scanning curve.
Definition: ParkerLenhard.hpp:189
Scalar SwMic()
Apparent saturation of the last reversal point on the pressure MIC.
Definition: ParkerLenhard.hpp:209
Implements the Parker-Lenhard twophase p_c-Sw hysteresis model.
Definition: ParkerLenhard.hpp:240
static void reset(Params &params)
Resets the hysteresis model to the initial parameters on the main drainage curve.
Definition: ParkerLenhard.hpp:289
static Evaluation Sw(const Params &, const FluidState &)
Calculate the wetting phase saturations depending on the phase pressures.
Definition: ParkerLenhard.hpp:424
static Evaluation krn(const Params &params, const FluidState &fs)
The relative permeability for the non-wetting phase of the params.
Definition: ParkerLenhard.hpp:470
static const bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition: ParkerLenhard.hpp:258
static Evaluation krw(const Params &params, const FluidState &fs)
The relative permeability for the wetting phase of the medium.
Definition: ParkerLenhard.hpp:448
static void update(Params &params, const FluidState &fs)
Set the current absolute saturation for the current timestep.
Definition: ParkerLenhard.hpp:303
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
Returns the relative permeabilities of the phases dependening on the phase saturations.
Definition: ParkerLenhard.hpp:364
static const bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition: ParkerLenhard.hpp:266
static const bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition: ParkerLenhard.hpp:270
static Evaluation absoluteToApparentSw_(const Params &params, const Evaluation &Sw)
Convert an absolute wetting saturation to an apparent one.
Definition: ParkerLenhard.hpp:491
static Evaluation Sn(const Params &params, const FluidState &fs)
Calculate the non-wetting phase saturations depending on the phase pressures.
Definition: ParkerLenhard.hpp:436
static Evaluation pcnw(const Params &params, const FluidState &fs)
Returns the capillary pressure dependend on the phase saturations.
Definition: ParkerLenhard.hpp:377
static const int numPhases
The number of fluid phases.
Definition: ParkerLenhard.hpp:247
static const bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition: ParkerLenhard.hpp:262
static const bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition: ParkerLenhard.hpp:254
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
Returns the capillary pressure dependening on the phase saturations.
Definition: ParkerLenhard.hpp:343
static void saturations(Container &, const Params &, const FluidState &)
Returns the capillary pressure dependening on the phase saturations.
Definition: ParkerLenhard.hpp:356
static const bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition: ParkerLenhard.hpp:274
static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation &Sw)
The saturation-capillary pressure curve according to van Genuchten using a material law specific API.
Definition: VanGenuchten.hpp:194