51 using EffLawParams =
typename EffLawT::Params;
52 using Scalar =
typename EffLawParams::Traits::Scalar;
55 using Traits =
typename EffLawParams::Traits;
60 result.deltaSwImbKrn_ = 1.0;
64 result.initialImb_ =
true;
65 result.pcSwMic_ = 3.0;
66 result.krnSwMdc_ = 4.0;
67 result.krwSwMdc_ = 4.5;
80 if (
config().enableHysteresis()) {
81 if (
config().krHysteresisModel() == 2 ||
config().krHysteresisModel() == 3 ||
config().krHysteresisModel() == 4 ||
config().pcHysteresisModel() == 0) {
82 C_ = 1.0/(Sncri_ - Sncrd_ + 1.0e-12) - 1.0/(Snmaxd_ - Sncrd_);
85 if (
config().krHysteresisModel() == 4) {
86 Cw_ = 1.0/(Swcri_ - Swcrd_ + 1.0e-12) - 1.0/(Swmaxd_ - Swcrd_);
87 Krwd_sncri_ = EffLawT::twoPhaseSatKrw(
drainageParams(), 1 - Sncri());
89 updateDynamicParams_();
91 EnsureFinalized :: finalize();
97 void setConfig(std::shared_ptr<EclHysteresisConfig> value)
109 void setWagConfig(std::shared_ptr<WagHysteresisConfig::WagHysteresisConfigRecord> value)
119 {
return *wagConfig_; }
128 drainageParams_ = value;
130 oilWaterSystem_ = (twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
131 gasOilSystem_ = (twoPhaseSystem == EclTwoPhaseSystemType::GasOil);
137 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
138 Sncrd_ = info.Sgcr + info.Swl;
140 Snmaxd_ = info.Sgu + info.Swl;
141 KrndMax_ = EffLawT::twoPhaseSatKrn(
drainageParams(), 1.0-Snmaxd_);
143 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
147 KrndMax_ = EffLawT::twoPhaseSatKrn(
drainageParams(), 1.0-Snmaxd_);
150 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
153 Snmaxd_ = 1.0 - info.Swl - info.Sgl;
154 KrndMax_ = EffLawT::twoPhaseSatKrn(
drainageParams(), 1.0-Snmaxd_);
158 if (
config().krHysteresisModel() == 4) {
160 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
161 Swmaxd_ = 1.0 - info.Sgl - info.Swl;
164 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
169 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
176 if (
config().pcHysteresisModel() == 0) {
177 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
178 pcmaxd_ = info.maxPcgo;
179 }
else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
180 pcmaxd_ = info.maxPcgo + info.maxPcow;
183 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
189 if (gasOilHysteresisWAG()) {
190 swatImbStart_ = Swco_;
191 swatImbStartNxt_ = -1.0;
193 krnSwDrainStart_ = Sncrd_;
194 krnSwDrainStartNxt_ = Sncrd_;
196 krnImbStartNxt_ = 0.0;
197 krnDrainStart_ = 0.0;
198 krnDrainStartNxt_ = 0.0;
201 krnSwImbStart_ = Sncrd_;
211 {
return drainageParams_; }
214 {
return drainageParams_; }
223 imbibitionParams_ = value;
225 if (!
config().enableHysteresis())
229 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
230 Sncri_ = info.Sgcr + info.Swl;
233 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
238 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
244 if (
config().pcHysteresisModel() == 0) {
245 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
246 Swmaxi_ = 1.0 - info.Sgl - info.Swl;
247 pcmaxi_ = info.maxPcgo;
248 }
else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
249 Swmaxi_ = 1.0 - info.Sgl;
250 pcmaxi_ = info.maxPcgo + info.maxPcow;
253 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
255 pcmaxi_ = info.maxPcow;
259 if (
config().krHysteresisModel() == 4) {
269 {
return imbibitionParams_; }
272 {
return imbibitionParams_; }
281 Scalar pcSwMic()
const
288 {
return initialImb_; }
296 { krwSwMdc_ = value; };
304 {
return krwSwMdc_; };
312 { krnSwMdc_ = value; }
320 {
return krnSwMdc_; }
350 { deltaSwImbKrn_ = value; }
360 {
return deltaSwImbKrn_; }
369 Scalar Swmaxi()
const
384 Scalar SnTrapped(
bool maximumTrapping)
const
386 if(!maximumTrapping && isDrain_)
390 if(
config().krHysteresisModel() > 1 )
393 return Sncri_ + deltaSwImbKrn_;
396 Scalar SnStranded(Scalar sg, Scalar krg)
const {
397 const Scalar sn = EffLawT::twoPhaseSatKrnInv(drainageParams_, krg);
398 return sg - (1.0 - sn) + Sncrd_;
401 Scalar SwTrapped()
const
403 if(
config().krHysteresisModel() == 0 ||
config().krHysteresisModel() == 2)
406 if(
config().krHysteresisModel() == 1 ||
config().krHysteresisModel() == 3)
410 if(
config().krHysteresisModel() == 4 )
418 Scalar SncrtWAG()
const
419 {
return SncrtWAG_; }
421 Scalar Snmaxd()
const
424 Scalar Swmaxd()
const
428 {
return 1.0 - krnSwMdc_; }
431 {
return krwSwMdc_; }
436 Scalar krnWght()
const
437 {
return KrndHy_/KrndMax_; }
439 Scalar krwWght()
const
442 Scalar deltaKrw = Krwi_snrmax() - Krwd_sncri();
443 Scalar Krwi_snr = Krwd_sncrt() + deltaKrw * (Sncrt() / max(1e-12, Sncri()));
444 return (Krwi_snr - KrwdHy()) / ( Krwi_snrmax() - Krwi_snmax());
447 Scalar krwdMax()
const
451 Scalar KrwdHy()
const
457 Scalar Krwd_sncri()
const
462 Scalar Krwi_snmax()
const
467 Scalar Krwi_snrmax()
const
472 Scalar Krwd_sncrt()
const
477 Scalar pcWght() const
480 return EffLawT::twoPhaseSatPcnw(
drainageParams(), 0.0)/(pcmaxi_+1e-6);
482 return pcmaxd_/(pcmaxi_+1e-6);
485 Scalar curvatureCapPrs()
const
486 {
return curvatureCapPrs_;}
488 bool gasOilHysteresisWAG()
const
489 {
return (
config().enableWagHysteresis() && gasOilSystem_ &&
wagConfig().wagGasFlag()) ; }
491 Scalar reductionDrain()
const
492 {
return std::pow(Swco_/(swatImbStart_+tolWAG_*
wagConfig().wagWaterThresholdSaturation()),
wagConfig().wagSecondaryDrainageReduction());}
494 Scalar reductionDrainNxt()
const
495 {
return std::pow(Swco_/(swatImbStartNxt_+tolWAG_*
wagConfig().wagWaterThresholdSaturation()),
wagConfig().wagSecondaryDrainageReduction());}
497 bool threePhaseState()
const
498 {
return (swatImbStart_ > (Swco_ +
wagConfig().wagWaterThresholdSaturation()) ); }
500 Scalar nState()
const
503 Scalar krnSwDrainRevert()
const
504 {
return krnSwDrainRevert_;}
506 Scalar krnDrainStart()
const
507 {
return krnDrainStart_;}
509 Scalar krnDrainStartNxt()
const
510 {
return krnDrainStartNxt_;}
512 Scalar krnImbStart()
const
513 {
return krnImbStart_;}
515 Scalar krnImbStartNxt()
const
516 {
return krnImbStartNxt_;}
518 Scalar krnSwWAG()
const
521 Scalar krnSwDrainStart()
const
522 {
return krnSwDrainStart_;}
524 Scalar krnSwDrainStartNxt()
const
525 {
return krnSwDrainStartNxt_;}
527 Scalar krnSwImbStart()
const
528 {
return krnSwImbStart_;}
530 Scalar tolWAG()
const
533 template <
class Evaluation>
534 Evaluation computeSwf(
const Evaluation& Sw)
const
536 Evaluation SgT = 1.0 - Sw - SncrtWAG();
537 Scalar SgCut =
wagConfig().wagImbCurveLinearFraction()*(Snhy()- SncrtWAG());
538 Evaluation Swf = 1.0;
543 Swf -= (Sncrd() + 0.5*( SgT + Opm::sqrt( SgT*SgT + 4.0/C*SgT)));
546 SgCut = std::max(Scalar(0.000001), SgCut);
547 Scalar SgCutValue = 0.5*( SgCut + Opm::sqrt( SgCut*SgCut + 4.0/C*SgCut));
548 Scalar SgCutSlope = SgCutValue/SgCut;
550 Swf -= (Sncrd() + SgT);
556 template <
class Evaluation>
557 Evaluation computeKrImbWAG(
const Evaluation& Sw)
const
561 Swf = computeSwf(Sw);
562 if (Swf <= krnSwDrainStart_) {
563 Evaluation Krg = EffLawT::twoPhaseSatKrn(drainageParams_, Swf);
564 Evaluation KrgImb2 = (Krg-krnDrainStart_)*reductionDrain() + krnImbStart_;
568 Evaluation Sn = Sncrd_;
569 if (Swf < 1.0-SncrtWAG_) {
571 Evaluation dd = (1.0-krnSwImbStart_ - Sncrd_) / (1.0-krnSwDrainStart_ - SncrtWAG_);
572 Sn += (1.0-Swf-SncrtWAG_)*dd;
574 Evaluation KrgDrn1 = EffLawT::twoPhaseSatKrn(drainageParams_, 1.0 - Sn);
585 bool update(Scalar pcSw, Scalar krwSw, Scalar krnSw)
587 bool updateParams =
false;
589 if (
config().pcHysteresisModel() == 0 && pcSw < pcSwMdc_) {
590 if (pcSwMdc_ == 2.0 && pcSw+1.0e-6 < Swcrd_ && oilWaterSystem_) {
597 if (initialImb_ && pcSw > pcSwMic_) {
602 if (krnSw < krnSwMdc_) {
605 if (
config().krHysteresisModel() == 4) {
610 if (krwSw > krwSwMdc_) {
616 if (!gasOilHysteresisWAG()) {
617 this->isDrain_ = (krnSw <= this->krnSwMdc_);
619 wasDrain_ = isDrain_;
621 if (swatImbStartNxt_ < 0.0) {
622 swatImbStartNxt_ = std::max(Swco_, Swco_ + krnSw - krwSw);
625 if ( (swatImbStartNxt_ > Swco_ + tolWAG_) && krwSw > tolWAG_) {
626 swatImbStart_ = swatImbStartNxt_;
628 krnSwDrainStartNxt_ = krnSwWAG_;
629 krnSwDrainStart_ = krnSwDrainStartNxt_;
635 if (krnSw <= krnSwWAG_+tolWAG_) {
636 krnSwWAG_ = std::min(krnSw, krnSwWAG_);
637 krnSwDrainRevert_ = krnSwWAG_;
647 if (krnSw >= krnSwWAG_-tolWAG_) {
648 krnSwWAG_ = std::max(krnSw, krnSwWAG_);
649 krnSwDrainStartNxt_ = krnSwWAG_;
650 swatImbStartNxt_ = std::max(swatImbStartNxt_, Swco_ + krnSw - krwSw);
655 krnSwDrainStart_ = krnSwDrainStartNxt_;
656 swatImbStart_ = swatImbStartNxt_;
665 updateDynamicParams_();
670 template<
class Serializer>
674 serializer(deltaSwImbKrn_);
678 serializer(initialImb_);
679 serializer(pcSwMic_);
680 serializer(krnSwMdc_);
681 serializer(krwSwMdc_);
686 bool operator==(
const EclHysteresisTwoPhaseLawParams& rhs)
const
688 return this->deltaSwImbKrn_ == rhs.deltaSwImbKrn_ &&
690 this->Sncrt_ == rhs.Sncrt_ &&
691 this->Swcrt_ == rhs.Swcrt_ &&
692 this->initialImb_ == rhs.initialImb_ &&
693 this->pcSwMic_ == rhs.pcSwMic_ &&
694 this->krnSwMdc_ == rhs.krnSwMdc_ &&
695 this->krwSwMdc_ == rhs.krwSwMdc_ &&
696 this->KrndHy_ == rhs.KrndHy_ &&
697 this->KrwdHy_ == rhs.KrwdHy_;
701 void updateDynamicParams_()
710 if (
config().krHysteresisModel() == 0 ||
config().krHysteresisModel() == 1) {
711 Scalar krnMdcDrainage = EffLawT::twoPhaseSatKrn(
drainageParams(), krnSwMdc_);
712 Scalar SwKrnMdcImbibition = EffLawT::twoPhaseSatKrnInv(
imbibitionParams(), krnMdcDrainage);
713 deltaSwImbKrn_ = SwKrnMdcImbibition - krnSwMdc_;
720 if (
config().krHysteresisModel() == 2 ||
721 config().krHysteresisModel() == 3 ||
722 config().krHysteresisModel() == 4 ||
723 config().pcHysteresisModel() == 0)
725 const Scalar snhy = 1.0 - krnSwMdc_;
727 Sncrt_ = Sncrd_ + (snhy - Sncrd_) /
728 ((1.0 +
config().modParamTrapped()*(Snmaxd_ - snhy)) + C_ * (snhy - Sncrd_));
735 if (
config().krHysteresisModel() == 4) {
736 Scalar swhy = krnSwMdc_;
737 if (swhy >= Swcrd_) {
738 Swcrt_ = Swcrd_ + (swhy - Swcrd_) /
739 ((1.0 +
config().modParamTrapped() * (Swmaxd_ - swhy)) + Cw_ * (swhy - Swcrd_));
743 Krwd_sncrt_ = EffLawT::twoPhaseSatKrw(
drainageParams(), 1 - Sncrt());
747 if (gasOilHysteresisWAG()) {
748 if (isDrain_ && krnSwMdc_ == krnSwWAG_) {
749 Scalar snhy = 1.0 - krnSwMdc_;
752 SncrtWAG_ += (snhy - Sncrd_) /
754 wagConfig().wagLandsParam() * (snhy - Sncrd_));
758 if (isDrain_ && (1.0 - krnSwDrainRevert_) > SncrtWAG_) {
759 cTransf_ = 1.0 / (SncrtWAG_ - Sncrd_ + 1.0e-12) - 1.0 / (1.0 - krnSwDrainRevert_ - Sncrd_);
762 if (!wasDrain_ && isDrain_) {
763 if (threePhaseState() || nState_ > 1) {
765 krnDrainStart_ = EffLawT::twoPhaseSatKrn(
drainageParams(), krnSwDrainStart_);
766 krnImbStart_ = krnImbStartNxt_;
768 krnSwImbStart_ = EffLawT::twoPhaseSatKrnInv(
drainageParams(), krnImbStart_);
772 if (!wasDrain_ && !isDrain_) {
773 krnDrainStartNxt_ = EffLawT::twoPhaseSatKrn(
drainageParams(), krnSwWAG_);
774 if (threePhaseState()) {
775 krnImbStartNxt_ = computeKrImbWAG(krnSwWAG_);
778 Scalar swf = computeSwf(krnSwWAG_);
787 EclHysteresisConfig config_{};
788 std::shared_ptr<WagHysteresisConfig::WagHysteresisConfigRecord> wagConfig_{};
789 EffLawParams imbibitionParams_{};
790 EffLawParams drainageParams_{};
795 Scalar krwSwMdc_{-2.0};
796 Scalar krnSwMdc_{2.0};
797 Scalar pcSwMdc_{2.0};
800 Scalar pcSwMic_{1.0};
802 bool initialImb_{
false};
804 bool oilWaterSystem_{
false};
805 bool gasOilSystem_{
false};
812 Scalar deltaSwImbKrn_{};
846 Scalar Krwd_sncri_{};
847 Scalar Krwi_snmax_{};
848 Scalar Krwi_snrmax_{};
849 Scalar Krwd_sncrt_{};
855 Scalar curvatureCapPrs_{};
861 Scalar swatImbStart_{};
862 Scalar swatImbStartNxt_{};
863 Scalar krnSwWAG_{2.0};
864 Scalar krnSwDrainRevert_{2.0};
867 Scalar krnSwDrainStart_{-2.0};
868 Scalar krnSwDrainStartNxt_{};
869 Scalar krnImbStart_{};
870 Scalar krnImbStartNxt_{};
871 Scalar krnDrainStart_{};
872 Scalar krnDrainStartNxt_{};
875 Scalar krnSwImbStart_{};
880 Scalar tolWAG_{0.001};