51 using EffLawParams =
typename EffLawT::Params;
52 using Scalar =
typename EffLawParams::Traits::Scalar;
55 using Traits =
typename EffLawParams::Traits;
65 krnSwDrainRevert_ = 2.0;
66 krnSwDrainStart_ = -2.0;
71 oilWaterSystem_ =
false;
72 gasOilSystem_ =
false;
90 result.deltaSwImbKrn_ = 1.0;
94 result.initialImb_ =
true;
95 result.pcSwMic_ = 3.0;
96 result.krnSwMdc_ = 4.0;
97 result.krwSwMdc_ = 4.5;
110 if (
config().enableHysteresis()) {
111 if (
config().krHysteresisModel() == 2 ||
config().krHysteresisModel() == 3 ||
config().krHysteresisModel() == 4 ||
config().pcHysteresisModel() == 0) {
112 C_ = 1.0/(Sncri_ - Sncrd_ + 1.0e-12) - 1.0/(Snmaxd_ - Sncrd_);
115 if (
config().krHysteresisModel() == 4) {
116 Cw_ = 1.0/(Swcri_ - Swcrd_ + 1.0e-12) - 1.0/(Swmaxd_ - Swcrd_);
118 updateDynamicParams_();
121 EnsureFinalized :: finalize();
127 void setConfig(std::shared_ptr<EclHysteresisConfig> value)
128 { config_ = *value; }
139 void setWagConfig(std::shared_ptr<WagHysteresisConfig::WagHysteresisConfigRecord> value)
149 {
return *wagConfig_; }
158 drainageParams_ = value;
160 oilWaterSystem_ = (twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
161 gasOilSystem_ = (twoPhaseSystem == EclTwoPhaseSystemType::GasOil);
167 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
168 Sncrd_ = info.Sgcr+info.Swl;
169 Swcrd_ = info.Sogcr + info.Swl;
170 Snmaxd_ = info.Sgu+info.Swl;
171 KrndMax_ = EffLawT::twoPhaseSatKrn(
drainageParams(), 1.0-Snmaxd_);
173 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
177 KrndMax_ = EffLawT::twoPhaseSatKrn(
drainageParams(), 1.0-Snmaxd_);
180 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
183 Snmaxd_ = 1.0 - info.Swl - info.Sgl;
184 KrndMax_ = EffLawT::twoPhaseSatKrn(
drainageParams(), 1.0-Snmaxd_);
188 if (
config().krHysteresisModel() == 4) {
190 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
191 Swmaxd_ = 1.0 - info.Sgl - info.Swl;
194 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
199 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
206 if (
config().pcHysteresisModel() == 0) {
207 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
208 pcmaxd_ = info.maxPcgo;
209 }
else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
210 pcmaxd_ = info.maxPcgo + info.maxPcow;
213 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
219 if (gasOilHysteresisWAG()) {
220 swatImbStart_ = Swco_;
221 swatImbStartNxt_ = -1.0;
223 krnSwDrainStart_ = Sncrd_;
224 krnSwDrainStartNxt_ = Sncrd_;
226 krnImbStartNxt_ = 0.0;
227 krnDrainStart_ = 0.0;
228 krnDrainStartNxt_ = 0.0;
231 krnSwImbStart_ = Sncrd_;
241 {
return drainageParams_; }
244 {
return drainageParams_; }
253 imbibitionParams_ = value;
255 if (!
config().enableHysteresis())
259 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
260 Sncri_ = info.Sgcr+info.Swl;
263 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
268 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
275 if (
config().pcHysteresisModel() == 0) {
276 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
277 Swmaxi_ = 1.0 - info.Sgl - info.Swl;
278 pcmaxi_ = info.maxPcgo;
279 }
else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
280 Swmaxi_ = 1.0 - info.Sgl;
281 pcmaxi_ = info.maxPcgo + info.maxPcow;
284 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
286 pcmaxi_ = info.maxPcow;
295 {
return imbibitionParams_; }
298 {
return imbibitionParams_; }
307 Scalar pcSwMic()
const
314 {
return initialImb_; }
322 { krwSwMdc_ = value; };
330 {
return krwSwMdc_; };
338 { krnSwMdc_ = value; }
346 {
return krnSwMdc_; }
376 { deltaSwImbKrn_ = value; }
386 {
return deltaSwImbKrn_; }
395 Scalar Swmaxi()
const
410 Scalar SnTrapped()
const
417 if(
config().krHysteresisModel() > 1 )
420 return Sncri_ + deltaSwImbKrn_;
423 Scalar SwTrapped()
const
425 if(
config().krHysteresisModel() == 0 ||
config().krHysteresisModel() == 2)
428 if(
config().krHysteresisModel() == 1 ||
config().krHysteresisModel() == 3)
432 if(
config().krHysteresisModel() == 4 )
440 Scalar SncrtWAG()
const
441 {
return SncrtWAG_; }
443 Scalar Snmaxd()
const
446 Scalar Swmaxd()
const
450 {
return 1.0 - krnSwMdc_; }
453 {
return krwSwMdc_; }
458 Scalar krnWght()
const
459 {
return KrndHy_/KrndMax_; }
461 Scalar krwWght()
const
463 return KrwdHy_/KrwdMax_; }
465 Scalar krwdMax()
const
469 Scalar pcWght() const
472 return EffLawT::twoPhaseSatPcnw(
drainageParams(), 0.0)/(pcmaxi_+1e-6);
474 return pcmaxd_/(pcmaxi_+1e-6);
477 Scalar curvatureCapPrs()
const
478 {
return curvatureCapPrs_;}
480 bool gasOilHysteresisWAG()
const
481 {
return (
config().enableWagHysteresis() && gasOilSystem_ &&
wagConfig().wagGasFlag()) ; }
483 Scalar reductionDrain()
const
484 {
return std::pow(Swco_/(swatImbStart_+tolWAG_*
wagConfig().wagWaterThresholdSaturation()),
wagConfig().wagSecondaryDrainageReduction());}
486 Scalar reductionDrainNxt()
const
487 {
return std::pow(Swco_/(swatImbStartNxt_+tolWAG_*
wagConfig().wagWaterThresholdSaturation()),
wagConfig().wagSecondaryDrainageReduction());}
489 bool threePhaseState()
const
490 {
return (swatImbStart_ > (Swco_ +
wagConfig().wagWaterThresholdSaturation()) ); }
492 Scalar nState()
const
495 Scalar krnSwDrainRevert()
const
496 {
return krnSwDrainRevert_;}
498 Scalar krnDrainStart()
const
499 {
return krnDrainStart_;}
501 Scalar krnDrainStartNxt()
const
502 {
return krnDrainStartNxt_;}
504 Scalar krnImbStart()
const
505 {
return krnImbStart_;}
507 Scalar krnImbStartNxt()
const
508 {
return krnImbStartNxt_;}
510 Scalar krnSwWAG()
const
513 Scalar krnSwDrainStart()
const
514 {
return krnSwDrainStart_;}
516 Scalar krnSwDrainStartNxt()
const
517 {
return krnSwDrainStartNxt_;}
519 Scalar krnSwImbStart()
const
520 {
return krnSwImbStart_;}
522 Scalar tolWAG()
const
525 template <
class Evaluation>
526 Evaluation computeSwf(
const Evaluation& Sw)
const
528 Evaluation SgT = 1.0 - Sw - SncrtWAG();
529 Scalar SgCut =
wagConfig().wagImbCurveLinearFraction()*(Snhy()- SncrtWAG());
530 Evaluation Swf = 1.0;
535 Swf -= (Sncrd() + 0.5*( SgT + Opm::sqrt( SgT*SgT + 4.0/C*SgT)));
538 SgCut = std::max(Scalar(0.000001), SgCut);
539 Scalar SgCutValue = 0.5*( SgCut + Opm::sqrt( SgCut*SgCut + 4.0/C*SgCut));
540 Scalar SgCutSlope = SgCutValue/SgCut;
542 Swf -= (Sncrd() + SgT);
548 template <
class Evaluation>
549 Evaluation computeKrImbWAG(
const Evaluation& Sw)
const
553 Swf = computeSwf(Sw);
554 if (Swf <= krnSwDrainStart_) {
555 Evaluation Krg = EffLawT::twoPhaseSatKrn(drainageParams_, Swf);
556 Evaluation KrgImb2 = (Krg-krnDrainStart_)*reductionDrain() + krnImbStart_;
560 Evaluation Sn = Sncrd_;
561 if (Swf < 1.0-SncrtWAG_) {
563 Evaluation dd = (1.0-krnSwImbStart_ - Sncrd_) / (1.0-krnSwDrainStart_ - SncrtWAG_);
564 Sn += (1.0-Swf-SncrtWAG_)*dd;
566 Evaluation KrgDrn1 = EffLawT::twoPhaseSatKrn(drainageParams_, 1.0 - Sn);
577 bool update(Scalar pcSw, Scalar krwSw, Scalar krnSw)
579 bool updateParams =
false;
581 if (
config().pcHysteresisModel() == 0 && pcSw < pcSwMdc_) {
582 if (pcSwMdc_ == 2.0 && pcSw+1.0e-6 < Swcrd_ && oilWaterSystem_) {
589 if (initialImb_ && pcSw > pcSwMic_) {
594 if (krnSw < krnSwMdc_) {
600 if (krwSw > krwSwMdc_) {
608 if (!gasOilHysteresisWAG()) {
609 if (krnSw < krnSwMdc_) {
615 wasDrain_ = isDrain_;
617 if (swatImbStartNxt_ < 0.0) {
618 swatImbStartNxt_ = std::max(Swco_, Swco_ + krnSw - krwSw);
621 if ( (swatImbStartNxt_ > Swco_ + tolWAG_) && krwSw > tolWAG_) {
622 swatImbStart_ = swatImbStartNxt_;
624 krnSwDrainStartNxt_ = krnSwWAG_;
625 krnSwDrainStart_ = krnSwDrainStartNxt_;
631 if (krnSw <= krnSwWAG_+tolWAG_) {
632 krnSwWAG_ = std::min(krnSw, krnSwWAG_);
633 krnSwDrainRevert_ = krnSwWAG_;
643 if (krnSw >= krnSwWAG_-tolWAG_) {
644 krnSwWAG_ = std::max(krnSw, krnSwWAG_);
645 krnSwDrainStartNxt_ = krnSwWAG_;
646 swatImbStartNxt_ = std::max(swatImbStartNxt_, Swco_ + krnSw - krwSw);
651 krnSwDrainStart_ = krnSwDrainStartNxt_;
652 swatImbStart_ = swatImbStartNxt_;
661 updateDynamicParams_();
666 template<
class Serializer>
670 serializer(deltaSwImbKrn_);
674 serializer(initialImb_);
675 serializer(pcSwMic_);
676 serializer(krnSwMdc_);
677 serializer(krwSwMdc_);
682 bool operator==(
const EclHysteresisTwoPhaseLawParams& rhs)
const
684 return this->deltaSwImbKrn_ == rhs.deltaSwImbKrn_ &&
686 this->Sncrt_ == rhs.Sncrt_ &&
687 this->Swcrt_ == rhs.Swcrt_ &&
688 this->initialImb_ == rhs.initialImb_ &&
689 this->pcSwMic_ == rhs.pcSwMic_ &&
690 this->krnSwMdc_ == rhs.krnSwMdc_ &&
691 this->krwSwMdc_ == rhs.krwSwMdc_ &&
692 this->KrndHy_ == rhs.KrndHy_ &&
693 this->KrwdHy_ == rhs.KrwdHy_;
697 void updateDynamicParams_()
706 if (
config().krHysteresisModel() == 0 ||
config().krHysteresisModel() == 1) {
707 Scalar krnMdcDrainage = EffLawT::twoPhaseSatKrn(
drainageParams(), krnSwMdc_);
708 Scalar SwKrnMdcImbibition = EffLawT::twoPhaseSatKrnInv(
imbibitionParams(), krnMdcDrainage);
709 deltaSwImbKrn_ = SwKrnMdcImbibition - krnSwMdc_;
716 if (
config().krHysteresisModel() == 2 ||
config().krHysteresisModel() == 3 ||
config().krHysteresisModel() == 4 ||
config().pcHysteresisModel() == 0) {
717 Scalar Snhy = 1.0 - krnSwMdc_;
719 Sncrt_ = Sncrd_ + (Snhy - Sncrd_)/((1.0+
config().modParamTrapped()*(Snmaxd_-Snhy)) + C_*(Snhy - Sncrd_));
727 if (
config().krHysteresisModel() == 4) {
728 Scalar Swhy = krwSwMdc_;
729 if (Swhy >= Swcrd_) {
730 Swcrt_ = Swcrd_ + (Swhy - Swcrd_)/((1.0+
config().modParamTrapped()*(Swmaxd_-Swhy)) + Cw_*(Swhy - Swcrd_));
738 if (gasOilHysteresisWAG()) {
739 if (isDrain_ && krnSwMdc_ == krnSwWAG_) {
740 Scalar Snhy = 1.0 - krnSwMdc_;
747 if (isDrain_ && (1.0-krnSwDrainRevert_) > SncrtWAG_) {
748 cTransf_ = 1.0/(SncrtWAG_-Sncrd_ + 1.0e-12) - 1.0/(1.0-krnSwDrainRevert_-Sncrd_);
751 if (!wasDrain_ && isDrain_) {
752 if (threePhaseState() || nState_>1) {
754 krnDrainStart_ = EffLawT::twoPhaseSatKrn(
drainageParams(), krnSwDrainStart_);
755 krnImbStart_ = krnImbStartNxt_;
757 krnSwImbStart_ = EffLawT::twoPhaseSatKrnInv(
drainageParams(), krnImbStart_);
761 if (!wasDrain_ && !isDrain_) {
762 krnDrainStartNxt_ = EffLawT::twoPhaseSatKrn(
drainageParams(), krnSwWAG_);
763 if (threePhaseState()) {
764 krnImbStartNxt_ = computeKrImbWAG(krnSwWAG_);
767 Scalar swf = computeSwf(krnSwWAG_);
776 EclHysteresisConfig config_{};
777 std::shared_ptr<WagHysteresisConfig::WagHysteresisConfigRecord> wagConfig_{};
778 EffLawParams imbibitionParams_{};
779 EffLawParams drainageParams_{};
793 bool oilWaterSystem_;
801 Scalar deltaSwImbKrn_;
836 Scalar curvatureCapPrs_{};
843 Scalar swatImbStart_;
844 Scalar swatImbStartNxt_{};
846 Scalar krnSwDrainRevert_;
849 Scalar krnSwDrainStart_;
850 Scalar krnSwDrainStartNxt_{};
851 Scalar krnImbStart_{};
852 Scalar krnImbStartNxt_{};
853 Scalar krnDrainStart_{};
854 Scalar krnDrainStartNxt_{};
857 Scalar krnSwImbStart_{};