26 #ifndef OPM_PARKER_LENHARD_HPP
27 #define OPM_PARKER_LENHARD_HPP
45 template <
class ScalarT>
160 return this->
Sw() < Sw && Sw < prev_->
Sw();
166 return prev_->
Sw() < Sw && Sw < this->
Sw();
237 template <
class TraitsT,
class ParamsT = ParkerLenhardParams<TraitsT> >
247 static_assert(numPhases == 2,
248 "The Parker-Lenhard capillary pressure law only "
249 "applies to the case of two fluid phases");
253 static const bool implementsTwoPhaseApi =
true;
257 static const bool implementsTwoPhaseSatApi =
true;
261 static const bool isSaturationDependent =
true;
265 static const bool isPressureDependent =
false;
269 static const bool isTemperatureDependent =
false;
273 static const bool isCompositionDependent =
false;
275 static_assert(Traits::numPhases == 2,
276 "The number of fluid phases must be two if you want to use "
277 "this material law!");
288 static void reset(Params ¶ms)
291 params.setMdc(
new ScanningCurve(params.SwrPc()));
292 params.setCsc(params.mdc());
293 params.setPisc(NULL);
294 params.setCurrentSnr(0.0);
301 template <
class Flu
idState>
302 static void update(Params ¶ms,
const FluidState &fs)
306 Scalar Sw = FsToolbox::value(fs.saturation(Traits::wettingPhaseIdx));
317 ScanningCurve *curve = findScanningCurve_(params, Sw);
322 Scalar pc = pcnw<FluidState, Scalar>(params, fs);
323 Scalar Sw_mic = VanGenuchten::twoPhaseSatSw(params.micParams(), pc);
324 Scalar Sw_mdc = VanGenuchten::twoPhaseSatSw(params.mdcParams(), pc);
326 curve->setNext(Sw, pc, Sw_mic, Sw_mdc);
330 params.setCsc(curve);
333 if (params.csc() == params.mdc()) {
334 params.setPisc(params.mdc()->next());
335 params.setCurrentSnr(computeCurrentSnr_(params, Sw));
343 template <
class Container,
class Flu
idState>
344 static void capillaryPressures(Container &values,
const Params ¶ms,
const FluidState &fs)
346 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
348 values[Traits::wettingPhaseIdx] = 0.0;
349 values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
356 template <
class Container,
class Flu
idState>
357 static void saturations(Container &,
const Params &,
const FluidState &)
358 { OPM_THROW(std::logic_error,
"Not implemented: ParkerLenhard::saturations()"); }
364 template <
class Container,
class Flu
idState>
365 static void relativePermeabilities(Container &values,
const Params ¶ms,
const FluidState &fs)
367 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
369 values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
370 values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
377 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
378 static Evaluation pcnw(
const Params ¶ms,
const FluidState &fs)
380 typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
382 const Evaluation& Sw =
383 FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
385 return twoPhaseSatPcnw(params, Sw);
388 template <
class Evaluation>
389 static Evaluation twoPhaseSatPcnw(
const Params ¶ms,
const Evaluation& Sw)
391 typedef MathToolbox<Evaluation> Toolbox;
394 ScanningCurve *sc = findScanningCurve_(params, Toolbox::value(Sw));
397 const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
407 Scalar SwAppCurSC = absoluteToApparentSw_(params, sc->Sw());
408 Scalar SwAppPrevSC = absoluteToApparentSw_(params, sc->prev()->Sw());
409 const Evaluation& pos = (Sw_app - SwAppCurSC)/(SwAppPrevSC - SwAppCurSC);
411 const Evaluation& SwMic =
412 pos * (sc->prev()->SwMic() - sc->SwMic()) + sc->SwMic();
414 return VanGenuchten::twoPhaseSatPcnw(params.micParams(), SwMic);
417 const Evaluation& SwMdc =
418 pos*(sc->prev()->SwMdc() - sc->SwMdc()) + sc->SwMdc();
420 return VanGenuchten::twoPhaseSatPcnw(params.mdcParams(), SwMdc);
428 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
429 static Evaluation Sw(
const Params &,
const FluidState &)
430 { OPM_THROW(std::logic_error,
"Not implemented: ParkerLenhard::Sw()"); }
432 template <
class Evaluation>
433 static Evaluation twoPhaseSatSw(
const Params &,
const Evaluation& )
434 { OPM_THROW(std::logic_error,
"Not implemented: ParkerLenhard::twoPhaseSatSw()"); }
440 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
441 static Evaluation Sn(
const Params ¶ms,
const FluidState &fs)
442 {
return 1 - Sw<FluidState, Evaluation>(params, fs); }
444 template <
class Evaluation>
445 static Evaluation twoPhaseSatSn(
const Params &,
const Evaluation& )
446 { OPM_THROW(std::logic_error,
"Not implemented: ParkerLenhard::twoPhaseSatSn()"); }
452 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
453 static Evaluation krw(
const Params ¶ms,
const FluidState &fs)
455 typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
457 const Evaluation& Sw =
458 FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
460 return twoPhaseSatKrw(params, Sw);
463 template <
class Evaluation>
464 static Evaluation twoPhaseSatKrw(
const Params ¶ms,
const Evaluation& Sw)
468 const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
469 return VanGenuchten::twoPhaseSatKrw(params.mdcParams(), Sw_app);
476 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
477 static Evaluation krn(
const Params ¶ms,
const FluidState &fs)
479 typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
481 const Evaluation& Sw =
482 FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
484 return twoPhaseSatKrn(params, Sw);
487 template <
class Evaluation>
488 static Evaluation twoPhaseSatKrn(
const Params ¶ms,
const Evaluation& Sw)
492 const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
493 return VanGenuchten::twoPhaseSatKrn(params.mdcParams(), Sw_app);
499 template <
class Evaluation>
500 static Evaluation absoluteToApparentSw_(
const Params ¶ms,
const Evaluation& Sw)
502 return effectiveToApparentSw_(params, absoluteToEffectiveSw_(params, Sw));
515 template <
class Evaluation>
516 static Evaluation absoluteToEffectiveSw_(
const Params ¶ms,
const Evaluation& Sw)
517 {
return (Sw - params.SwrPc())/(1 - params.SwrPc()); }
528 template <
class Evaluation>
529 static Evaluation effectiveToAbsoluteSw_(
const Params ¶ms,
const Evaluation& Swe)
530 {
return Swe*(1 - params.SwrPc()) + params.SwrPc(); }
534 template <
class Evaluation>
535 static Evaluation computeCurrentSnr_(
const Params ¶ms,
const Evaluation& Sw)
538 if (Sw > 1 - params.Snr())
540 if (Sw < params.SwrPc())
543 if (params.Snr() == 0.0)
547 Scalar R = 1.0/params.Snr() - 1;
548 const Evaluation& curSnr = (1 - Sw)/(1 + R*(1 - Sw));
552 assert(curSnr <= params.Snr());
559 template <
class Evaluation>
560 static Evaluation trappedEffectiveSn_(
const Params ¶ms,
const Evaluation& Sw)
562 const Evaluation& Swe = absoluteToEffectiveSw_(params, Sw);
563 Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
565 Scalar Snre = absoluteToEffectiveSw_(params, params.currentSnr());
566 return Snre*(Swe - SwePisc) / (1 - Snre - SwePisc);
571 template <
class Evaluation>
572 static Evaluation effectiveToApparentSw_(
const Params ¶ms,
const Evaluation& Swe)
574 if (params.pisc() == NULL ||
575 Swe <= absoluteToEffectiveSw_(params, params.pisc()->Sw()))
586 return Swe + trappedEffectiveSn_(params, Swe);
590 template <
class Evaluation>
591 static Evaluation apparentToEffectiveSw_(
const Params ¶ms,
const Evaluation& Swapp)
593 Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
594 if (params.pisc() == NULL || Swapp <= SwePisc) {
601 Scalar Snre = absoluteToEffectiveSw_(params.currentSnr());
603 (Swapp*(1 - Snre - SwePisc) + Snre*SwePisc)
610 static ScanningCurve *findScanningCurve_(
const Params ¶ms, Scalar Sw)
612 if (params.pisc() == NULL || Sw <= params.pisc()->Sw()) {
623 if (params.pisc()->next() == NULL ||
624 params.pisc()->next()->Sw() < Sw)
626 return params.pisc();
629 ScanningCurve *curve = params.csc()->next();
631 assert(curve != params.mdc()->prev());
632 if (curve->isValidAt_Sw(Sw)) {
635 curve = curve->prev();
642 #endif // PARKER_LENHARD_HPP
Scalar SwMic()
Apparent saturation of the last reversal point on the pressure MIC.
Definition: ParkerLenhard.hpp:208
static const int numPhases
The number of fluid phases.
Definition: ParkerLenhard.hpp:246
Definition: Air_Mesitylene.hpp:31
Implements the Parker-Lenhard twophase p_c-Sw hysteresis model. This class adheres to the twophase ca...
Definition: ParkerLenhard.hpp:238
Implementation of the van Genuchten capillary pressure - saturation relation.
Definition: VanGenuchten.hpp:54
Implementation of the van Genuchten capillary pressure - saturation relation.
void setNext(Scalar Sw, Scalar pcnw, Scalar SwMic, Scalar SwMdc)
Set the next scanning curve.
Definition: ParkerLenhard.hpp:129
Traits::Scalar Scalar
Definition: ParkerLenhard.hpp:243
PLScanningCurve(Scalar Swr)
Constructs main imbibition curve.
Definition: ParkerLenhard.hpp:57
ScalarT Scalar
Definition: ParkerLenhard.hpp:49
Scalar pcnw() const
Capillary pressure at the last reversal point.
Definition: ParkerLenhard.hpp:201
bool isDrain()
Returns true iff the scanning curve is a drainage curve.
Definition: ParkerLenhard.hpp:180
int loopNum()
The loop number of the scanning curve.
Definition: ParkerLenhard.hpp:188
Default parameter class for the Parker-Lenhard hysteresis model.
~PLScanningCurve()
Destructor. After it was called all references to the next() curve are invalid!
Definition: ParkerLenhard.hpp:99
bool isValidAt_Sw(Scalar Sw)
Returns true iff the given effective saturation Swei is within the scope of the curve, i.e. whether Swei is part of the curve's domain and the curve thus applies to Swi.
Definition: ParkerLenhard.hpp:153
Scalar Sw() const
Absolute wetting-phase saturation at the scanning curve's reversal point.
Definition: ParkerLenhard.hpp:195
Represents a scanning curve in the Parker-Lenhard hysteresis model.
Definition: ParkerLenhard.hpp:46
PLScanningCurve * prev() const
Return the previous scanning curve, i.e. the curve with one less reversal than the current one...
Definition: ParkerLenhard.hpp:111
Scalar SwMdc()
Apparent saturation of the last reversal point on the pressure MDC.
Definition: ParkerLenhard.hpp:215
PLScanningCurve(PLScanningCurve *prev, PLScanningCurve *next, int loopN, Scalar Sw, Scalar pcnw, Scalar SwMic, Scalar SwMdc)
Definition: ParkerLenhard.hpp:76
ParamsT Params
Definition: ParkerLenhard.hpp:242
PLScanningCurve * next() const
Return the next scanning curve, i.e. the curve with one more reversal than the current one...
Definition: ParkerLenhard.hpp:118
TraitsT Traits
Definition: ParkerLenhard.hpp:241
bool isImbib()
Returns true iff the scanning curve is a imbibition curve.
Definition: ParkerLenhard.hpp:173