27#ifndef OPM_PARKER_LENHARD_HPP
28#define OPM_PARKER_LENHARD_HPP
46template <
class ScalarT>
50 typedef ScalarT Scalar;
163 return this->
Sw() < SwReversal && SwReversal < prev_->
Sw();
169 return prev_->
Sw() < SwReversal && SwReversal < this->
Sw();
240template <
class TraitsT,
class ParamsT = ParkerLenhardParams<TraitsT> >
244 typedef TraitsT Traits;
245 typedef ParamsT Params;
246 typedef typename Traits::Scalar Scalar;
251 "The Parker-Lenhard capillary pressure law only "
252 "applies to the case of two fluid phases");
278 static_assert(Traits::numPhases == 2,
279 "The number of fluid phases must be two if you want to use "
280 "this material law!");
283 typedef typename ParamsT::VanGenuchten VanGenuchten;
295 params.setCsc(params.mdc());
296 params.setPisc(
nullptr);
297 params.setCurrentSnr(0.0);
304 template <
class Flu
idState>
305 static void update(Params& params,
const FluidState& fs)
307 Scalar
Sw = scalarValue(fs.saturation(Traits::wettingPhaseIdx));
323 Scalar pc = pcnw<FluidState, Scalar>(params, fs);
324 Scalar Sw_mic = VanGenuchten::twoPhaseSatSw(params.micParams(), pc);
325 Scalar Sw_mdc = VanGenuchten::twoPhaseSatSw(params.mdcParams(), pc);
331 params.setCsc(curve);
334 if (params.csc() == params.mdc()) {
335 params.setPisc(params.mdc()->next());
336 params.setCurrentSnr(computeCurrentSnr_(params,
Sw));
344 template <
class Container,
class Flu
idState>
347 typedef typename std::remove_reference<
decltype(values[0])>::type Evaluation;
349 values[Traits::wettingPhaseIdx] = 0.0;
350 values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
357 template <
class Container,
class Flu
idState>
358 static void saturations(Container& ,
const Params& ,
const FluidState& )
359 {
throw std::logic_error(
"Not implemented: ParkerLenhard::saturations()"); }
365 template <
class Container,
class Flu
idState>
368 typedef typename std::remove_reference<
decltype(values[0])>::type Evaluation;
370 values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
371 values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
378 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
379 static Evaluation
pcnw(
const Params& params,
const FluidState& fs)
381 const Evaluation&
Sw =
382 decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
384 return twoPhaseSatPcnw(params,
Sw);
387 template <
class Evaluation>
388 static Evaluation twoPhaseSatPcnw(
const Params& params,
const Evaluation&
Sw)
391 ScanningCurve* sc = findScanningCurve_(params, scalarValue(
Sw));
406 const Evaluation& pos = (Sw_app - SwAppCurSC)/(SwAppPrevSC - SwAppCurSC);
408 const Evaluation& SwMic =
409 pos * (sc->prev()->SwMic() - sc->SwMic()) + sc->SwMic();
414 const Evaluation& SwMdc =
415 pos*(sc->prev()->SwMdc() - sc->SwMdc()) + sc->SwMdc();
425 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
426 static Evaluation
Sw(
const Params& ,
const FluidState& )
427 {
throw std::logic_error(
"Not implemented: ParkerLenhard::Sw()"); }
429 template <
class Evaluation>
430 static Evaluation twoPhaseSatSw(
const Params& ,
const Evaluation& )
431 {
throw std::logic_error(
"Not implemented: ParkerLenhard::twoPhaseSatSw()"); }
437 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
438 static Evaluation
Sn(
const Params& params,
const FluidState& fs)
439 {
return 1 - Sw<FluidState, Evaluation>(params, fs); }
441 template <
class Evaluation>
442 static Evaluation twoPhaseSatSn(
const Params& ,
const Evaluation& )
443 {
throw std::logic_error(
"Not implemented: ParkerLenhard::twoPhaseSatSn()"); }
449 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
450 static Evaluation
krw(
const Params& params,
const FluidState& fs)
452 const Evaluation&
Sw =
453 decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
455 return twoPhaseSatKrw(params,
Sw);
458 template <
class Evaluation>
459 static Evaluation twoPhaseSatKrw(
const Params& params,
const Evaluation&
Sw)
464 return VanGenuchten::twoPhaseSatKrw(params.mdcParams(), Sw_app);
471 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
472 static Evaluation
krn(
const Params& params,
const FluidState& fs)
474 const Evaluation&
Sw =
475 decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
477 return twoPhaseSatKrn(params,
Sw);
480 template <
class Evaluation>
481 static Evaluation twoPhaseSatKrn(
const Params& params,
const Evaluation&
Sw)
486 return VanGenuchten::twoPhaseSatKrn(params.mdcParams(), Sw_app);
492 template <
class Evaluation>
495 return effectiveToApparentSw_(params, absoluteToEffectiveSw_(params,
Sw));
508 template <
class Evaluation>
509 static Evaluation absoluteToEffectiveSw_(
const Params& params,
const Evaluation&
Sw)
510 {
return (
Sw - params.SwrPc())/(1 - params.SwrPc()); }
521 template <
class Evaluation>
522 static Evaluation effectiveToAbsoluteSw_(
const Params& params,
const Evaluation& Swe)
523 {
return Swe*(1 - params.SwrPc()) + params.SwrPc(); }
527 template <
class Evaluation>
528 static Evaluation computeCurrentSnr_(
const Params& params,
const Evaluation&
Sw)
531 if (
Sw > 1 - params.Snr())
533 if (
Sw < params.SwrPc())
536 if (params.Snr() == 0.0)
540 Scalar R = 1.0/params.Snr() - 1;
541 const Evaluation& curSnr = (1 -
Sw)/(1 + R*(1 -
Sw));
545 assert(curSnr <= params.Snr());
552 template <
class Evaluation>
553 static Evaluation trappedEffectiveSn_(
const Params& params,
const Evaluation&
Sw)
555 const Evaluation& Swe = absoluteToEffectiveSw_(params,
Sw);
556 Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
558 Scalar Snre = absoluteToEffectiveSw_(params, params.currentSnr());
559 return Snre*(Swe - SwePisc) / (1 - Snre - SwePisc);
564 template <
class Evaluation>
565 static Evaluation effectiveToApparentSw_(
const Params& params,
const Evaluation& Swe)
567 if (params.pisc() ==
nullptr ||
568 Swe <= absoluteToEffectiveSw_(params, params.pisc()->Sw()))
579 return Swe + trappedEffectiveSn_(params, Swe);
583 template <
class Evaluation>
584 static Evaluation apparentToEffectiveSw_(
const Params& params,
const Evaluation& Swapp)
586 Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
587 if (params.pisc() ==
nullptr || Swapp <= SwePisc) {
594 Scalar Snre = absoluteToEffectiveSw_(params.currentSnr());
596 (Swapp*(1 - Snre - SwePisc) + Snre*SwePisc)
603 static ScanningCurve* findScanningCurve_(
const Params& params, Scalar
Sw)
605 if (params.pisc() ==
nullptr ||
Sw <= params.pisc()->Sw()) {
616 if (params.pisc()->next() ==
nullptr ||
617 params.pisc()->next()->Sw() <
Sw)
619 return params.pisc();
622 ScanningCurve* curve = params.csc()->next();
624 assert(curve != params.mdc()->prev());
625 if (curve->isValidAt_Sw(
Sw)) {
628 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 * next() const
Return the next scanning curve, i.e.
Definition ParkerLenhard.hpp:121
~PLScanningCurve()
Destructor.
Definition ParkerLenhard.hpp:102
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:183
void setNext(Scalar SwReversal, Scalar pcnwReversal, Scalar SwMiCurve, Scalar SwMdCurve)
Set the next scanning curve.
Definition ParkerLenhard.hpp:132
bool isImbib()
Returns true iff the scanning curve is a imbibition curve.
Definition ParkerLenhard.hpp:176
Scalar Sw() const
Absolute wetting-phase saturation at the scanning curve's reversal point.
Definition ParkerLenhard.hpp:198
bool isValidAt_Sw(Scalar SwReversal)
Returns true iff the given effective saturation Swei is within the scope of the curve,...
Definition ParkerLenhard.hpp:156
Scalar SwMdc()
Apparent saturation of the last reversal point on the pressure MDC.
Definition ParkerLenhard.hpp:218
Scalar pcnw() const
Capillary pressure at the last reversal point.
Definition ParkerLenhard.hpp:204
PLScanningCurve * prev() const
Return the previous scanning curve, i.e.
Definition ParkerLenhard.hpp:114
int loopNum()
The loop number of the scanning curve.
Definition ParkerLenhard.hpp:191
Scalar SwMic()
Apparent saturation of the last reversal point on the pressure MIC.
Definition ParkerLenhard.hpp:211
Implements the Parker-Lenhard twophase p_c-Sw hysteresis model.
Definition ParkerLenhard.hpp:242
static void reset(Params ¶ms)
Resets the hysteresis model to the initial parameters on the main drainage curve.
Definition ParkerLenhard.hpp:291
static Evaluation Sw(const Params &, const FluidState &)
Calculate the wetting phase saturations depending on the phase pressures.
Definition ParkerLenhard.hpp:426
static Evaluation krn(const Params ¶ms, const FluidState &fs)
The relative permeability for the non-wetting phase of the params.
Definition ParkerLenhard.hpp:472
static const bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition ParkerLenhard.hpp:260
static Evaluation krw(const Params ¶ms, const FluidState &fs)
The relative permeability for the wetting phase of the medium.
Definition ParkerLenhard.hpp:450
static void update(Params ¶ms, const FluidState &fs)
Set the current absolute saturation for the current timestep.
Definition ParkerLenhard.hpp:305
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:366
static const bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition ParkerLenhard.hpp:268
static const bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition ParkerLenhard.hpp:272
static Evaluation absoluteToApparentSw_(const Params ¶ms, const Evaluation &Sw)
Convert an absolute wetting saturation to an apparent one.
Definition ParkerLenhard.hpp:493
static Evaluation Sn(const Params ¶ms, const FluidState &fs)
Calculate the non-wetting phase saturations depending on the phase pressures.
Definition ParkerLenhard.hpp:438
static Evaluation pcnw(const Params ¶ms, const FluidState &fs)
Returns the capillary pressure dependend on the phase saturations.
Definition ParkerLenhard.hpp:379
static const int numPhases
The number of fluid phases.
Definition ParkerLenhard.hpp:249
static const bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition ParkerLenhard.hpp:264
static const bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition ParkerLenhard.hpp:256
static void capillaryPressures(Container &values, const Params ¶ms, const FluidState &fs)
Returns the capillary pressure dependening on the phase saturations.
Definition ParkerLenhard.hpp:345
static void saturations(Container &, const Params &, const FluidState &)
Returns the capillary pressure dependening on the phase saturations.
Definition ParkerLenhard.hpp:358
static const bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition ParkerLenhard.hpp:276
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
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30