27#ifndef OPM_ECL_HYSTERESIS_TWO_PHASE_LAW_PARAMS_HPP
28#define OPM_ECL_HYSTERESIS_TWO_PHASE_LAW_PARAMS_HPP
35#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
47template <
class EffLawT>
50 using EffLawParams =
typename EffLawT::Params;
51 using Scalar =
typename EffLawParams::Traits::Scalar;
54 using Traits =
typename EffLawParams::Traits;
63 krnSwDrainRevert_ = 2.0;
64 krnSwDrainStart_ = -2.0;
71 oilWaterSystem_ =
false;
72 gasOilSystem_ =
false;
90 result.deltaSwImbKrn_ = 1.0;
92 result.initialImb_ =
true;
93 result.pcSwMic_ = 3.0;
94 result.krnSwMdc_ = 4.0;
106 if (
config().enableHysteresis()) {
107 if (
config().krHysteresisModel() == 2 ||
config().krHysteresisModel() == 3 ||
config().pcHysteresisModel() == 0) {
108 C_ = 1.0/(Sncri_ - Sncrd_ + 1.0e-12) - 1.0/(Snmaxd_ - Sncrd_);
112 updateDynamicParams_();
115 EnsureFinalized :: finalize();
121 void setConfig(std::shared_ptr<EclHysteresisConfig> value)
122 { config_ = *value; }
133 void setWagConfig(std::shared_ptr<WagHysteresisConfig::WagHysteresisConfigRecord> value)
143 {
return *wagConfig_; }
152 drainageParams_ = value;
154 oilWaterSystem_ = (twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
155 gasOilSystem_ = (twoPhaseSystem == EclTwoPhaseSystemType::GasOil);
161 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
162 Sncrd_ = info.Sgcr+info.Swl;
163 Snmaxd_ = info.Sgu+info.Swl;
164 KrndMax_ = EffLawT::twoPhaseSatKrn(
drainageParams(), 1.0-Snmaxd_);
166 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
169 KrndMax_ = EffLawT::twoPhaseSatKrn(
drainageParams(), 1.0-Snmaxd_);
172 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
174 Snmaxd_ = 1.0 - info.Swl - info.Sgl;
175 KrndMax_ = EffLawT::twoPhaseSatKrn(
drainageParams(), 1.0-Snmaxd_);
180 if (
config().pcHysteresisModel() == 0) {
181 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
183 pcmaxd_ = info.maxPcgo;
184 }
else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
186 pcmaxd_ = info.maxPcgo + info.maxPcow;
189 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
196 if (gasOilHysteresisWAG()) {
197 swatImbStart_ = Swco_;
198 swatImbStartNxt_ = -1.0;
200 krnSwDrainStart_ = Sncrd_;
201 krnSwDrainStartNxt_ = Sncrd_;
203 krnImbStartNxt_ = 0.0;
204 krnDrainStart_ = 0.0;
205 krnDrainStartNxt_ = 0.0;
208 krnSwImbStart_ = Sncrd_;
218 {
return drainageParams_; }
221 {
return drainageParams_; }
230 imbibitionParams_ = value;
232 if (!
config().enableHysteresis())
236 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
237 Sncri_ = info.Sgcr+info.Swl;
239 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
243 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
249 if (
config().pcHysteresisModel() == 0) {
250 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
252 Swmaxi_ = 1.0 - info.Sgl - info.Swl;
253 pcmaxi_ = info.maxPcgo;
254 }
else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
256 Swmaxi_ = 1.0 - info.Sgl;
257 pcmaxi_ = info.maxPcgo + info.maxPcow;
260 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
263 pcmaxi_ = info.maxPcow;
272 {
return imbibitionParams_; }
275 {
return imbibitionParams_; }
284 Scalar pcSwMic()
const
291 {
return initialImb_; }
317 { krnSwMdc_ = value; }
325 {
return krnSwMdc_; }
357 { deltaSwImbKrn_ = value; }
367 {
return deltaSwImbKrn_; }
376 Scalar Swmaxi()
const
388 Scalar SnTrapped()
const
391 if(
config().krHysteresisModel() > 1 )
394 return Sncri_ + deltaSwImbKrn_;
396 Scalar SncrtWAG()
const
397 {
return SncrtWAG_; }
399 Scalar Snmaxd()
const
403 {
return 1.0 - krnSwMdc_; }
408 Scalar krnWght()
const
409 {
return KrndHy_/KrndMax_; }
411 Scalar pcWght() const
414 return EffLawT::twoPhaseSatPcnw(
drainageParams(), 0.0)/(pcmaxi_+1e-6);
416 return pcmaxd_/(pcmaxi_+1e-6);
419 Scalar curvatureCapPrs()
const
420 {
return curvatureCapPrs_;}
422 bool gasOilHysteresisWAG()
const
423 {
return (
config().enableWagHysteresis() && gasOilSystem_ &&
wagConfig().wagGasFlag()) ; }
425 Scalar reductionDrain()
const
426 {
return std::pow(Swco_/(swatImbStart_+tolWAG_*
wagConfig().wagWaterThresholdSaturation()),
wagConfig().wagSecondaryDrainageReduction());}
428 Scalar reductionDrainNxt()
const
429 {
return std::pow(Swco_/(swatImbStartNxt_+tolWAG_*
wagConfig().wagWaterThresholdSaturation()),
wagConfig().wagSecondaryDrainageReduction());}
431 bool threePhaseState()
const
432 {
return (swatImbStart_ > (Swco_ +
wagConfig().wagWaterThresholdSaturation()) ); }
434 Scalar nState()
const
437 Scalar krnSwDrainRevert()
const
438 {
return krnSwDrainRevert_;}
440 Scalar krnDrainStart()
const
441 {
return krnDrainStart_;}
443 Scalar krnDrainStartNxt()
const
444 {
return krnDrainStartNxt_;}
446 Scalar krnImbStart()
const
447 {
return krnImbStart_;}
449 Scalar krnImbStartNxt()
const
450 {
return krnImbStartNxt_;}
452 Scalar krnSwWAG()
const
455 Scalar krnSwDrainStart()
const
456 {
return krnSwDrainStart_;}
458 Scalar krnSwDrainStartNxt()
const
459 {
return krnSwDrainStartNxt_;}
461 Scalar krnSwImbStart()
const
462 {
return krnSwImbStart_;}
464 Scalar tolWAG()
const
467 template <
class Evaluation>
468 Evaluation computeSwf(
const Evaluation& Sw)
const
470 Evaluation SgT = 1.0 - Sw - SncrtWAG();
471 Scalar SgCut =
wagConfig().wagImbCurveLinearFraction()*(Snhy()- SncrtWAG());
472 Evaluation Swf = 1.0;
477 Swf -= (Sncrd() + 0.5*( SgT + Opm::sqrt( SgT*SgT + 4.0/C*SgT)));
480 SgCut = std::max(Scalar(0.000001), SgCut);
481 Scalar SgCutValue = 0.5*( SgCut + Opm::sqrt( SgCut*SgCut + 4.0/C*SgCut));
482 Scalar SgCutSlope = SgCutValue/SgCut;
484 Swf -= (Sncrd() + SgT);
490 template <
class Evaluation>
491 Evaluation computeKrImbWAG(
const Evaluation& Sw)
const
495 Swf = computeSwf(Sw);
496 if (Swf <= krnSwDrainStart_) {
497 Evaluation Krg = EffLawT::twoPhaseSatKrn(drainageParams_, Swf);
498 Evaluation KrgImb2 = (Krg-krnDrainStart_)*reductionDrain() + krnImbStart_;
502 Evaluation Sn = Sncrd_;
503 if (Swf < 1.0-SncrtWAG_) {
505 Evaluation dd = (1.0-krnSwImbStart_ - Sncrd_) / (1.0-krnSwDrainStart_ - SncrtWAG_);
506 Sn += (1.0-Swf-SncrtWAG_)*dd;
508 Evaluation KrgDrn1 = EffLawT::twoPhaseSatKrn(drainageParams_, 1.0 - Sn);
519 bool update(Scalar pcSw, Scalar krwSw, Scalar krnSw)
521 bool updateParams =
false;
523 if (
config().pcHysteresisModel() == 0 && pcSw < pcSwMdc_) {
524 if (pcSwMdc_ == 2.0 && pcSw+1.0e-6 < Swcrd_ && oilWaterSystem_) {
531 if (initialImb_ && pcSw > pcSwMic_) {
549 if (krnSw < krnSwMdc_) {
555 if (gasOilHysteresisWAG()) {
557 wasDrain_ = isDrain_;
559 if (swatImbStartNxt_ < 0.0) {
560 swatImbStartNxt_ = std::max(Swco_, Swco_ + krnSw - krwSw);
563 if ( (swatImbStartNxt_ > Swco_ + tolWAG_) && krwSw > tolWAG_) {
564 swatImbStart_ = swatImbStartNxt_;
566 krnSwDrainStartNxt_ = krnSwWAG_;
567 krnSwDrainStart_ = krnSwDrainStartNxt_;
573 if (krnSw <= krnSwWAG_+tolWAG_) {
574 krnSwWAG_ = std::min(krnSw, krnSwWAG_);
575 krnSwDrainRevert_ = krnSwWAG_;
585 if (krnSw >= krnSwWAG_-tolWAG_) {
586 krnSwWAG_ = std::max(krnSw, krnSwWAG_);
587 krnSwDrainStartNxt_ = krnSwWAG_;
588 swatImbStartNxt_ = std::max(swatImbStartNxt_, Swco_ + krnSw - krwSw);
593 krnSwDrainStart_ = krnSwDrainStartNxt_;
594 swatImbStart_ = swatImbStartNxt_;
603 updateDynamicParams_();
608 template<
class Serializer>
612 serializer(deltaSwImbKrn_);
614 serializer(initialImb_);
615 serializer(pcSwMic_);
616 serializer(krnSwMdc_);
620 bool operator==(
const EclHysteresisTwoPhaseLawParams& rhs)
const
622 return this->deltaSwImbKrn_ == rhs.deltaSwImbKrn_ &&
623 this->Sncrt_ == rhs.Sncrt_ &&
624 this->initialImb_ == rhs.initialImb_ &&
625 this->pcSwMic_ == rhs.pcSwMic_ &&
626 this->krnSwMdc_ == rhs.krnSwMdc_ &&
627 this->KrndHy_ == rhs.KrndHy_;
631 void updateDynamicParams_()
641 if (
config().krHysteresisModel() == 0 ||
config().krHysteresisModel() == 1) {
642 Scalar krnMdcDrainage = EffLawT::twoPhaseSatKrn(
drainageParams(), krnSwMdc_);
643 Scalar SwKrnMdcImbibition = EffLawT::twoPhaseSatKrnInv(
imbibitionParams(), krnMdcDrainage);
644 deltaSwImbKrn_ = SwKrnMdcImbibition - krnSwMdc_;
645 assert(std::abs(EffLawT::twoPhaseSatKrn(
imbibitionParams(), krnSwMdc_ + deltaSwImbKrn_)
660 if (
config().krHysteresisModel() == 2 ||
config().krHysteresisModel() == 3 ||
config().pcHysteresisModel() == 0) {
661 Scalar Snhy = 1.0 - krnSwMdc_;
663 Sncrt_ = Sncrd_ + (Snhy - Sncrd_)/((1.0+
config().modParamTrapped()*(Snmaxd_-Snhy)) + C_*(Snhy - Sncrd_));
671 if (gasOilHysteresisWAG()) {
672 if (isDrain_ && krnSwMdc_ == krnSwWAG_) {
673 Scalar Snhy = 1.0 - krnSwMdc_;
680 if (isDrain_ && (1.0-krnSwDrainRevert_) > SncrtWAG_) {
681 cTransf_ = 1.0/(SncrtWAG_-Sncrd_ + 1.0e-12) - 1.0/(1.0-krnSwDrainRevert_-Sncrd_);
684 if (!wasDrain_ && isDrain_) {
685 if (threePhaseState() || nState_>1) {
687 krnDrainStart_ = EffLawT::twoPhaseSatKrn(
drainageParams(), krnSwDrainStart_);
688 krnImbStart_ = krnImbStartNxt_;
690 krnSwImbStart_ = EffLawT::twoPhaseSatKrnInv(
drainageParams(), krnImbStart_);
694 if (!wasDrain_ && !isDrain_) {
695 krnDrainStartNxt_ = EffLawT::twoPhaseSatKrn(
drainageParams(), krnSwWAG_);
696 if (threePhaseState()) {
697 krnImbStartNxt_ = computeKrImbWAG(krnSwWAG_);
700 Scalar swf = computeSwf(krnSwWAG_);
709 EclHysteresisConfig config_;
710 std::shared_ptr<WagHysteresisConfig::WagHysteresisConfigRecord> wagConfig_;
711 EffLawParams imbibitionParams_;
712 EffLawParams drainageParams_;
726 bool oilWaterSystem_;
734 Scalar deltaSwImbKrn_;
763 Scalar curvatureCapPrs_;
769 Scalar swatImbStart_;
770 Scalar swatImbStartNxt_;
772 Scalar krnSwDrainRevert_;
775 Scalar krnSwDrainStart_;
776 Scalar krnSwDrainStartNxt_;
778 Scalar krnImbStartNxt_;
779 Scalar krnDrainStart_;
780 Scalar krnDrainStartNxt_;
783 Scalar krnSwImbStart_;
Specifies the configuration used by the endpoint scaling code.
Specifies the configuration used by the ECL kr/pC hysteresis code.
Default implementation for asserting finalization of parameter objects.
Specifies the configuration used by the ECL kr/pC hysteresis code.
Definition EclHysteresisConfig.hpp:42
int pcHysteresisModel() const
Return the type of the hysteresis model which is used for capillary pressure.
Definition EclHysteresisConfig.hpp:71
int krHysteresisModel() const
Return the type of the hysteresis model which is used for relative permeability.
Definition EclHysteresisConfig.hpp:97
double curvatureCapPrs() const
Curvature parameter used for capillary pressure hysteresis.
Definition EclHysteresisConfig.hpp:113
bool enableHysteresis() const
Returns whether hysteresis is enabled.
Definition EclHysteresisConfig.hpp:53
double modParamTrapped() const
Regularisation parameter used for Killough model.
Definition EclHysteresisConfig.hpp:105
A default implementation of the parameters for the material law which implements the ECL relative per...
Definition EclHysteresisTwoPhaseLawParams.hpp:49
Scalar pcSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition EclHysteresisTwoPhaseLawParams.hpp:281
bool initialImb() const
Status of initial process.
Definition EclHysteresisTwoPhaseLawParams.hpp:290
const EclHysteresisConfig & config() const
Returns the hysteresis configuration object.
Definition EclHysteresisTwoPhaseLawParams.hpp:127
void setDeltaSwImbKrn(Scalar value)
Sets the saturation value which must be added if krn is calculated using the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:356
void setConfig(std::shared_ptr< EclHysteresisConfig > value)
Set the hysteresis configuration object.
Definition EclHysteresisTwoPhaseLawParams.hpp:121
Scalar krnSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition EclHysteresisTwoPhaseLawParams.hpp:324
bool update(Scalar pcSw, Scalar krwSw, Scalar krnSw)
Notify the hysteresis law that a given wetting-phase saturation has been seen.
Definition EclHysteresisTwoPhaseLawParams.hpp:519
void setKrnSwMdc(Scalar value)
Set the saturation of the wetting phase where the last switch from the main drainage curve (MDC) to i...
Definition EclHysteresisTwoPhaseLawParams.hpp:316
void setDeltaSwImbKrw(Scalar)
Sets the saturation value which must be added if krw is calculated using the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:334
void setWagConfig(std::shared_ptr< WagHysteresisConfig::WagHysteresisConfigRecord > value)
Set the WAG-hysteresis configuration object.
Definition EclHysteresisTwoPhaseLawParams.hpp:133
void setKrwSwMdc(Scalar)
Set the saturation of the wetting phase where the last switch from the main drainage curve (MDC) to i...
Definition EclHysteresisTwoPhaseLawParams.hpp:298
void finalize()
Calculate all dependent quantities once the independent quantities of the parameter object have been ...
Definition EclHysteresisTwoPhaseLawParams.hpp:104
void setDrainageParams(const EffLawParams &value, const EclEpsScalingPointsInfo< Scalar > &info, EclTwoPhaseSystemType twoPhaseSystem)
Sets the parameters used for the drainage curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:148
const WagHysteresisConfig::WagHysteresisConfigRecord & wagConfig() const
Returns the WAG-hysteresis configuration object.
Definition EclHysteresisTwoPhaseLawParams.hpp:142
Scalar deltaSwImbKrn() const
Returns the saturation value which must be added if krn is calculated using the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:366
void setImbibitionParams(const EffLawParams &value, const EclEpsScalingPointsInfo< Scalar > &info, EclTwoPhaseSystemType twoPhaseSystem)
Sets the parameters used for the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:226
Scalar deltaSwImbKrw() const
Returns the saturation value which must be added if krw is calculated using the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:345
const EffLawParams & imbibitionParams() const
Returns the parameters used for the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:271
Scalar krwSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition EclHysteresisTwoPhaseLawParams.hpp:307
const EffLawParams & drainageParams() const
Returns the parameters used for the drainage curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:217
Default implementation for asserting finalization of parameter objects.
Definition EnsureFinalized.hpp:47
Class for (de-)serializing.
Definition Serializer.hpp:84
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30
EclTwoPhaseSystemType
Specified which fluids are involved in a given twophase material law for endpoint scaling.
Definition EclEpsConfig.hpp:42
This structure represents all values which can be possibly used as scaling points in the endpoint sca...
Definition EclEpsScalingPoints.hpp:57
Definition WagHysteresisConfig.hpp:30