27 #ifndef OPM_PARKER_LENHARD_HPP
28 #define OPM_PARKER_LENHARD_HPP
46 template <
class ScalarT>
50 typedef ScalarT Scalar;
161 return this->
Sw() < SwReversal && SwReversal < prev_->
Sw();
167 return prev_->
Sw() < SwReversal && SwReversal < this->
Sw();
238 template <
class TraitsT,
class ParamsT = ParkerLenhardParams<TraitsT> >
242 typedef TraitsT Traits;
243 typedef ParamsT Params;
244 typedef typename Traits::Scalar Scalar;
249 "The Parker-Lenhard capillary pressure law only "
250 "applies to the case of two fluid phases");
276 static_assert(Traits::numPhases == 2,
277 "The number of fluid phases must be two if you want to use "
278 "this material law!");
281 typedef typename ParamsT::VanGenuchten VanGenuchten;
293 params.setCsc(params.mdc());
294 params.setPisc(NULL);
295 params.setCurrentSnr(0.0);
302 template <
class Flu
idState>
303 static void update(Params& params,
const FluidState& fs)
305 Scalar
Sw = scalarValue(fs.saturation(Traits::wettingPhaseIdx));
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);
329 params.setCsc(curve);
332 if (params.csc() == params.mdc()) {
333 params.setPisc(params.mdc()->next());
334 params.setCurrentSnr(computeCurrentSnr_(params,
Sw));
342 template <
class Container,
class Flu
idState>
345 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
347 values[Traits::wettingPhaseIdx] = 0.0;
348 values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
355 template <
class Container,
class Flu
idState>
356 static void saturations(Container& ,
const Params& ,
const FluidState& )
357 {
throw std::logic_error(
"Not implemented: ParkerLenhard::saturations()"); }
363 template <
class Container,
class Flu
idState>
366 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
368 values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
369 values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
376 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
377 static Evaluation
pcnw(
const Params& params,
const FluidState& fs)
379 const Evaluation&
Sw =
380 decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
382 return twoPhaseSatPcnw(params,
Sw);
385 template <
class Evaluation>
386 static Evaluation twoPhaseSatPcnw(
const Params& params,
const Evaluation&
Sw)
389 ScanningCurve* sc = findScanningCurve_(params, scalarValue(
Sw));
404 const Evaluation& pos = (Sw_app - SwAppCurSC)/(SwAppPrevSC - SwAppCurSC);
406 const Evaluation& SwMic =
407 pos * (sc->prev()->SwMic() - sc->SwMic()) + sc->SwMic();
412 const Evaluation& SwMdc =
413 pos*(sc->prev()->SwMdc() - sc->SwMdc()) + sc->SwMdc();
423 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
424 static Evaluation
Sw(
const Params& ,
const FluidState& )
425 {
throw std::logic_error(
"Not implemented: ParkerLenhard::Sw()"); }
427 template <
class Evaluation>
428 static Evaluation twoPhaseSatSw(
const Params& ,
const Evaluation& )
429 {
throw std::logic_error(
"Not implemented: ParkerLenhard::twoPhaseSatSw()"); }
435 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
436 static Evaluation
Sn(
const Params& params,
const FluidState& fs)
437 {
return 1 - Sw<FluidState, Evaluation>(params, fs); }
439 template <
class Evaluation>
440 static Evaluation twoPhaseSatSn(
const Params& ,
const Evaluation& )
441 {
throw std::logic_error(
"Not implemented: ParkerLenhard::twoPhaseSatSn()"); }
447 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
448 static Evaluation
krw(
const Params& params,
const FluidState& fs)
450 const Evaluation&
Sw =
451 decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
453 return twoPhaseSatKrw(params,
Sw);
456 template <
class Evaluation>
457 static Evaluation twoPhaseSatKrw(
const Params& params,
const Evaluation&
Sw)
462 return VanGenuchten::twoPhaseSatKrw(params.mdcParams(), Sw_app);
469 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
470 static Evaluation
krn(
const Params& params,
const FluidState& fs)
472 const Evaluation&
Sw =
473 decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
475 return twoPhaseSatKrn(params,
Sw);
478 template <
class Evaluation>
479 static Evaluation twoPhaseSatKrn(
const Params& params,
const Evaluation&
Sw)
484 return VanGenuchten::twoPhaseSatKrn(params.mdcParams(), Sw_app);
490 template <
class Evaluation>
493 return effectiveToApparentSw_(params, absoluteToEffectiveSw_(params,
Sw));
506 template <
class Evaluation>
507 static Evaluation absoluteToEffectiveSw_(
const Params& params,
const Evaluation&
Sw)
508 {
return (
Sw - params.SwrPc())/(1 - params.SwrPc()); }
519 template <
class Evaluation>
520 static Evaluation effectiveToAbsoluteSw_(
const Params& params,
const Evaluation& Swe)
521 {
return Swe*(1 - params.SwrPc()) + params.SwrPc(); }
525 template <
class Evaluation>
526 static Evaluation computeCurrentSnr_(
const Params& params,
const Evaluation&
Sw)
529 if (
Sw > 1 - params.Snr())
531 if (
Sw < params.SwrPc())
534 if (params.Snr() == 0.0)
538 Scalar R = 1.0/params.Snr() - 1;
539 const Evaluation& curSnr = (1 -
Sw)/(1 + R*(1 -
Sw));
543 assert(curSnr <= params.Snr());
550 template <
class Evaluation>
551 static Evaluation trappedEffectiveSn_(
const Params& params,
const Evaluation&
Sw)
553 const Evaluation& Swe = absoluteToEffectiveSw_(params,
Sw);
554 Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
556 Scalar Snre = absoluteToEffectiveSw_(params, params.currentSnr());
557 return Snre*(Swe - SwePisc) / (1 - Snre - SwePisc);
562 template <
class Evaluation>
563 static Evaluation effectiveToApparentSw_(
const Params& params,
const Evaluation& Swe)
565 if (params.pisc() == NULL ||
566 Swe <= absoluteToEffectiveSw_(params, params.pisc()->Sw()))
577 return Swe + trappedEffectiveSn_(params, Swe);
581 template <
class Evaluation>
582 static Evaluation apparentToEffectiveSw_(
const Params& params,
const Evaluation& Swapp)
584 Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
585 if (params.pisc() == NULL || Swapp <= SwePisc) {
592 Scalar Snre = absoluteToEffectiveSw_(params.currentSnr());
594 (Swapp*(1 - Snre - SwePisc) + Snre*SwePisc)
601 static ScanningCurve* findScanningCurve_(
const Params& params, Scalar
Sw)
603 if (params.pisc() == NULL ||
Sw <= params.pisc()->Sw()) {
614 if (params.pisc()->next() == NULL ||
615 params.pisc()->next()->Sw() <
Sw)
617 return params.pisc();
620 ScanningCurve* curve = params.csc()->next();
622 assert(curve != params.mdc()->prev());
623 if (curve->isValidAt_Sw(
Sw)) {
626 curve = curve->prev();
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 ¶ms)
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 ¶ms, 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 ¶ms, const FluidState &fs)
The relative permeability for the wetting phase of the medium.
Definition: ParkerLenhard.hpp:448
static void update(Params ¶ms, const FluidState &fs)
Set the current absolute saturation for the current timestep.
Definition: ParkerLenhard.hpp:303
static void relativePermeabilities(Container &values, const Params ¶ms, 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 ¶ms, const Evaluation &Sw)
Convert an absolute wetting saturation to an apparent one.
Definition: ParkerLenhard.hpp:491
static Evaluation Sn(const Params ¶ms, const FluidState &fs)
Calculate the non-wetting phase saturations depending on the phase pressures.
Definition: ParkerLenhard.hpp:436
static Evaluation pcnw(const Params ¶ms, 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 ¶ms, 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 ¶ms, const Evaluation &Sw)
The saturation-capillary pressure curve according to van Genuchten using a material law specific API.
Definition: VanGenuchten.hpp:194