OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Some support for WAG hysteresis (material).

Browse Source
This commit is contained in:
Ove Sævareid
2023-04-05 20:33:41 +02:00
committed by Tor Harald Sandve
parent 0fc0b6f4fa
commit 34a9d9c949
9 changed files with 334 additions and 16 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp
View File

@@ -452,10 +452,9 @@ public:
changed = changed || oilchanged;
bool gaschanged = params.gasOilParams().update(/*pcSw=*/ 1.0 - Swco - Sg,
/*krwSw=*/ 1.0 - Swco - Sg,
/*krnSw=*/ 1.0 - Swco - Sg);
boool gaschanged = params.gasOilParams().update(/*pcSw=*/ 1.0 - Swco - Sg,
/*krwSw=*/ 1.0 - std::max(Swco, Sw) - Sg, //(Three phase behavior ...)
/*krnSw=*/ 1.0 - Swco - Sg);
changed = changed || gaschanged;
}
else {

9
opm/material/fluidmatrixinteractions/EclHysteresisConfig.hpp
View File

@@ -113,6 +113,12 @@ public:
double curvatureCapPrs() const
{ return curvatureCapPrs_; }
/*!
* \brief Returns whether hysteresis is enabled.
*/
bool enableWagHysteresis() const
{ return enableWagHyst_; }
#if HAVE_ECL_INPUT
/*!
* \brief Reads all relevant material parameters form a cell of a parsed ECL deck.
@@ -131,6 +137,9 @@ private:
int krHysteresisModel_{-1};
double modParamTrapped_{};
double curvatureCapPrs_{};
// WAG hysteresis
bool enableWagHyst_{false};
};
} // namespace Opm

37
opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLaw.hpp
View File

@@ -297,6 +297,43 @@ public:
if (!params.config().enableHysteresis() || params.config().krHysteresisModel() < 0)
return EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Sw);
// If WAG hysteresis is enabled, the convential hysteresis model is ignored.
// (Two-phase model, non-wetting: only gas in oil.)
if (params.gasOilHysteresisWAG()) {
// Primary drainage
if (Sw <= params.krnSwMdc()+params.tolWAG() && params.nState()==1) {
Evaluation krn = EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Sw);
return krn;
}
// Imbibition or reversion to two-phase drainage retracing imb curve
// (Shift along primary drainage curve.)
if (params.nState()==1) {
Evaluation Swf = params.computeSwf(Sw);
Evaluation krn = EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Swf);
return krn;
}
// Three-phase drainage along current secondary drainage curve
if (Sw <= params.krnSwDrainRevert()+params.tolWAG() /*&& params.nState()>=1 */) {
Evaluation Krg = EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Sw);
Evaluation KrgDrain2 = (Krg-params.krnDrainStart())*params.reductionDrain() + params.krnImbStart();
return KrgDrain2;
}
// Subsequent imbibition: Scanning curve derived from previous secondary drainage
if (Sw >= params.krnSwWAG()-params.tolWAG() /*&& Sw > params.krnSwDrainRevert() && params.nState()>=1 */) {
Evaluation KrgImb2 = params.computeKrImbWAG(Sw);
return KrgImb2;
}
else {/* Sw < params.krnSwWAG() */ // Reversion along "next" drainage curve
Evaluation Krg = EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Sw);
Evaluation KrgDrainNxt = (Krg-params.krnDrainStartNxt())*params.reductionDrainNxt() + params.krnImbStartNxt();
return KrgDrainNxt;
}
}
// if it is enabled, use either the drainage or the imbibition curve. if the
// imbibition curve is used, the saturation must be shifted.
if (Sw <= params.krnSwMdc())

268
opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLawParams.hpp
View File

@@ -59,16 +59,28 @@ public:
// using the drainage curve before the first saturation update.
pcSwMdc_ = 2.0;
krnSwMdc_ = 2.0;
krnSwDrainRevert_ = 2.0;
krnSwDrainStart_ = -2.0;
krnSwWAG_ = 2.0;
// krwSwMdc_ = 2.0;
pcSwMic_ = -1.0;
initialImb_ = false;
oilWaterSystem_ = false;
gasOilSystem_ = false;
pcmaxd_ = 0.0;
pcmaxi_ = 0.0;
deltaSwImbKrn_ = 0.0;
// deltaSwImbKrw_ = 0.0;
Swco_ = 0.0;
swatImbStart_ = 0.0;
isDrain_ = true;
cTransf_ = 0.0;
tolWAG_ = 0.001;
nState_ = 0;
}
static EclHysteresisTwoPhaseLawParams serializationTestObject()
@@ -103,17 +115,32 @@ public:
}
/*!
* \brief Set the endpoint scaling configuration object.
* \brief Set the hysteresis configuration object.
*/
void setConfig(std::shared_ptr<EclHysteresisConfig> value)
{ config_ = *value; }
/*!
* \brief Returns the endpoint scaling configuration object.
* \brief Returns the hysteresis configuration object.
*/
const EclHysteresisConfig& config() const
{ return config_; }
/*!
* \brief Set the WAG-hysteresis configuration object.
*/
void setWagConfig(std::shared_ptr<WagHysteresisConfig::WagHysteresisConfigRecord> value)
{
wagConfig_ = value;
cTransf_ = wagConfig().wagLandsParam();
}
/*!
* \brief Returns the WAG-hysteresis configuration object.
*/
const WagHysteresisConfig::WagHysteresisConfigRecord& wagConfig() const
{ return *wagConfig_; }
/*!
* \brief Sets the parameters used for the drainage curve
*/
@@ -124,11 +151,12 @@ public:
drainageParams_ = value;
oilWaterSystem_ = (twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
gasOilSystem_ = (twoPhaseSystem == EclTwoPhaseSystemType::GasOil);
if (!config().enableHysteresis())
return;
if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().pcHysteresisModel() == 0) {
if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().pcHysteresisModel() == 0 || gasOilHysteresisWAG()) {
Swco_ = info.Swl;
if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
Sncrd_ = info.Sgcr+info.Swl;
Snmaxd_ = info.Sgu+info.Swl;
@@ -162,6 +190,24 @@ public:
pcmaxd_ = -17.0; // At this point 'info.maxPcow' holds pre-swatinit value ...;
}
}
// For WAG hysteresis, assume initial state along primary drainage curve.
if (gasOilHysteresisWAG()) {
swatImbStart_ = Swco_;
swatImbStartNxt_ = -1.0; // Trigger check for saturation gt Swco at first update ...
cTransf_ = wagConfig().wagLandsParam();
krnSwDrainStart_ = Sncrd_;
krnSwDrainStartNxt_ = Sncrd_;
krnImbStart_ = 0.0;
krnImbStartNxt_ = 0.0;
krnDrainStart_ = 0.0;
krnDrainStartNxt_ = 0.0;
isDrain_ = true;
wasDrain_ = true;
krnSwImbStart_ = Sncrd_;
SncrtWAG_ = Sncrd_;
nState_ = 1;
}
}
/*!
@@ -346,6 +392,8 @@ public:
else // For Carlson we use the shift to compute it from the critial saturation
return Sncri_ + deltaSwImbKrn_;
}
Scalar SncrtWAG() const
{ return SncrtWAG_; }
Scalar Snmaxd() const
{ return Snmaxd_; }
@@ -353,6 +401,9 @@ public:
Scalar Snhy() const
{ return 1.0 - krnSwMdc_; }
Scalar Swco() const
{ return Swco_; }
Scalar krnWght() const
{ return KrndHy_/KrndMax_; }
@@ -367,6 +418,96 @@ public:
Scalar curvatureCapPrs() const
{ return curvatureCapPrs_;}
bool gasOilHysteresisWAG() const
{ return (config().enableWagHysteresis() && gasOilSystem_ && wagConfig().wagGasFlag()) ; }
Scalar reductionDrain() const
{ return std::pow(Swco_/(swatImbStart_+tolWAG_*wagConfig().wagWaterThresholdSaturation()), wagConfig().wagSecondaryDrainageReduction());}
Scalar reductionDrainNxt() const
{ return std::pow(Swco_/(swatImbStartNxt_+tolWAG_*wagConfig().wagWaterThresholdSaturation()), wagConfig().wagSecondaryDrainageReduction());}
bool threePhaseState() const
{ return (swatImbStart_ > (Swco_ + wagConfig().wagWaterThresholdSaturation()) ); }
Scalar nState() const
{ return nState_;}
Scalar krnSwDrainRevert() const
{ return krnSwDrainRevert_;}
Scalar krnDrainStart() const
{ return krnDrainStart_;}
Scalar krnDrainStartNxt() const
{ return krnDrainStartNxt_;}
Scalar krnImbStart() const
{ return krnImbStart_;}
Scalar krnImbStartNxt() const
{ return krnImbStartNxt_;}
Scalar krnSwWAG() const
{ return krnSwWAG_;}
Scalar krnSwDrainStart() const
{ return krnSwDrainStart_;}
Scalar krnSwDrainStartNxt() const
{ return krnSwDrainStartNxt_;}
Scalar krnSwImbStart() const
{ return krnSwImbStart_;}
Scalar tolWAG() const
{ return tolWAG_;}
template <class Evaluation>
Evaluation computeSwf(const Evaluation& Sw) const
{
Evaluation SgT = 1.0 - Sw - SncrtWAG(); // Sg-Sg_crit_trapped
Scalar SgCut = wagConfig().wagImbCurveLinearFraction()*(Snhy()- SncrtWAG());
Evaluation Swf = 1.0;
//Scalar C = wagConfig().wagLandsParam();
Scalar C = cTransf_;
if (SgT > SgCut) {
Swf -= (Sncrd() + 0.5*( SgT + Opm::sqrt( SgT*SgT + 4.0/C*SgT))); // 1-Sgf
}
else {
SgCut = std::max(Scalar(0.000001), SgCut);
Scalar SgCutValue = 0.5*( SgCut + Opm::sqrt( SgCut*SgCut + 4.0/C*SgCut));
Scalar SgCutSlope = SgCutValue/SgCut;
SgT *= SgCutSlope;
Swf -= (Sncrd() + SgT);
}
return Swf;
}
template <class Evaluation>
Evaluation computeKrImbWAG(const Evaluation& Sw) const
{
Evaluation Swf = Sw;
if (nState_ <= 2) // Skipping for "higher order" curves seems consistent with benchmark, further investigations needed ...
Swf = computeSwf(Sw);
if (Swf <= krnSwDrainStart_) { // Use secondary drainage curve
Evaluation Krg = EffLawT::twoPhaseSatKrn(drainageParams_, Swf);
Evaluation KrgImb2 = (Krg-krnDrainStart_)*reductionDrain() + krnImbStart_;
return KrgImb2;
}
else { // Fallback to primary drainage curve
Evaluation Sn = Sncrd_;
if (Swf < 1.0-SncrtWAG_) {
// Notation: Sn.. = Sg.. + Swco
Evaluation dd = (1.0-krnSwImbStart_ - Sncrd_) / (1.0-krnSwDrainStart_ - SncrtWAG_);
Sn += (1.0-Swf-SncrtWAG_)*dd;
}
Evaluation KrgDrn1 = EffLawT::twoPhaseSatKrn(drainageParams_, 1.0 - Sn);
return KrgDrn1;
}
}
/*!
* \brief Notify the hysteresis law that a given wetting-phase saturation has been seen
@@ -374,7 +515,7 @@ public:
* This updates the scanning curves and the imbibition<->drainage reversal points as
* appropriate.
*/
bool update(Scalar pcSw, Scalar /* krwSw */, Scalar krnSw)
void update(Scalar pcSw, Scalar krwSw, Scalar krnSw)
{
bool updateParams = false;
@@ -410,6 +551,51 @@ public:
updateParams = true;
}
if (gasOilHysteresisWAG()) {
wasDrain_ = isDrain_;
if (swatImbStartNxt_ < 0.0) { // Initial check ...
swatImbStartNxt_ = std::max(Swco_, Swco_ + krnSw - krwSw);
if (swatImbStartNxt_ > Swco_ + tolWAG_) {
swatImbStart_ = swatImbStartNxt_;
krnSwWAG_ = krnSw;
krnSwDrainStartNxt_ = krnSwWAG_;
krnSwDrainStart_ = krnSwDrainStartNxt_;
wasDrain_ = false; // Signal start from threephase state ...
}
}
if (isDrain_) {
if (krnSw <= krnSwWAG_+tolWAG_) { // continue along drainage curve
krnSwWAG_ = std::min(krnSw, krnSwWAG_);
krnSwDrainRevert_ = krnSwWAG_;
updateParams = true;
}
else { // start new imbibition curve
isDrain_ = false;
krnSwWAG_ = krnSw;
updateParams = true;
}
}
else {
if (krnSw >= krnSwWAG_-tolWAG_) { // continue along imbibition curve
krnSwWAG_ = std::max(krnSw, krnSwWAG_);
krnSwDrainStartNxt_ = krnSwWAG_;
swatImbStartNxt_ = std::max(swatImbStartNxt_, Swco_ + krnSw - krwSw);
updateParams = true;
}
else { // start new drainage curve
isDrain_ = true;
krnSwDrainStart_ = krnSwDrainStartNxt_;
swatImbStart_ = swatImbStartNxt_;
krnSwWAG_ = krnSw;
updateParams = true;
}
}
}
if (updateParams)
updateDynamicParams_();
@@ -477,16 +663,55 @@ private:
Sncrt_ = Sncrd_;
}
}
if (gasOilHysteresisWAG()) {
if (isDrain_ && krnSwMdc_ == krnSwWAG_) {
Scalar Snhy = 1.0 - krnSwMdc_;
SncrtWAG_ = Sncrd_;
if (Snhy > Sncrd_) {
SncrtWAG_ += (Snhy - Sncrd_)/(1.0+config().modParamTrapped()*(Snmaxd_-Snhy) + wagConfig().wagLandsParam()*(Snhy - Sncrd_));
}
}
if (isDrain_ && krnSwDrainRevert_ > SncrtWAG_) { //Reversal from drain to imb
cTransf_ = 1.0/(SncrtWAG_-Sncrd_ + 1.0e-12) - 1.0/(1.0-krnSwDrainRevert_-Sncrd_);
}
if (!wasDrain_ && isDrain_) { // Start of new drainage cycle
if (threePhaseState() || nState_>1) { // Never return to primary (two-phase) state after leaving
nState_ += 1;
krnDrainStart_ = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwDrainStart_);
krnImbStart_ = krnImbStartNxt_;
// Scanning shift for primary drainage
krnSwImbStart_ = EffLawT::twoPhaseSatKrnInv(drainageParams(), krnImbStart_);
}
}
if (!wasDrain_ && !isDrain_) { //Moving along current imb curve
krnDrainStartNxt_ = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwWAG_);
if (threePhaseState()) {
krnImbStartNxt_ = computeKrImbWAG(krnSwWAG_);
}
else {
Scalar swf = computeSwf(krnSwWAG_);
krnImbStartNxt_ = EffLawT::twoPhaseSatKrn(drainageParams(), swf);
}
}
}
}
EclHysteresisConfig config_;
std::shared_ptr<WagHysteresisConfig::WagHysteresisConfigRecord> wagConfig_;
EffLawParams imbibitionParams_;
EffLawParams drainageParams_;
// largest wettinging phase saturation which is on the main-drainage curve. These are
// three different values because the sourounding code can choose to use different
// definitions for the saturations for different quantities
//Scalar krwSwMdc_;
Scalar krwSwMdc_;
Scalar krnSwMdc_;
Scalar pcSwMdc_;
@@ -496,6 +721,8 @@ private:
bool initialImb_;
bool oilWaterSystem_;
bool gasOilSystem_;
// offsets added to wetting phase saturation uf using the imbibition curves need to
// be used to calculate the wetting phase relperm, the non-wetting phase relperm and
@@ -504,9 +731,6 @@ private:
Scalar deltaSwImbKrn_;
//Scalar deltaSwImbPc_;
// trapped non-wetting phase saturation
Scalar Sncrt_;
// the following uses the conventions of the Eclipse technical description:
//
// Sncrd_: critical non-wetting phase saturation for the drainage curve
@@ -518,6 +742,7 @@ private:
// maximum on the drainage curve
// C_: factor required to calculate the trapped non-wetting phase saturation using
// the Killough approach
// Swcod_: connate water saturation value used for wag hysteresis (2. drainage)
Scalar Sncrd_;
Scalar Sncri_;
Scalar Swcri_;
@@ -533,6 +758,31 @@ private:
Scalar pcmaxi_; // max pc for imb
Scalar curvatureCapPrs_; // curvature parameter used for capillary pressure hysteresis
Scalar Sncrt_; // trapped non-wetting phase saturation
// Used for WAG hysteresis
Scalar Swco_; // Connate water.
Scalar swatImbStart_; // Water saturation at start of current drainage curve (end of previous imb curve).
Scalar swatImbStartNxt_; // Water saturation at start of next drainage curve (end of current imb curve).
Scalar krnSwWAG_; // Saturation value after latest completed timestep.
Scalar krnSwDrainRevert_; // Saturation value at end of current drainage curve.
Scalar cTransf_; // Modified Lands constant used for free gas calculations to obtain consistent scanning curve
// when reversion to imb occurs above historical maximum gas saturation (i.e. Sw > krwSwMdc_).
Scalar krnSwDrainStart_; // Saturation value at start of current drainage curve (end of previous imb curve).
Scalar krnSwDrainStartNxt_; // Saturation value at start of current drainage curve (end of previous imb curve).
Scalar krnImbStart_; // Relperm at start of current drainage curve (end of previous imb curve).
Scalar krnImbStartNxt_; // Relperm at start of next drainage curve (end of current imb curve).
Scalar krnDrainStart_; // Primary (input) relperm evaluated at start of current drainage curve.
Scalar krnDrainStartNxt_; // Primary (input) relperm evaluated at start of next drainage curve.
bool isDrain_; // Status is either drainage or imbibition
bool wasDrain_; // Previous status.
Scalar krnSwImbStart_; // Saturation value where primary drainage relperm equals krnImbStart_
int nState_; // Number of cycles. Primary cycle is nState_=1.
Scalar SncrtWAG_;
Scalar tolWAG_;
};
} // namespace Opm

5
opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp
View File

@@ -35,6 +35,7 @@
#include <opm/material/fluidmatrixinteractions/SatCurveMultiplexer.hpp>
#include <opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp>
#include <opm/input/eclipse/EclipseState/WagHysteresisConfig.hpp>
#include <opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLaw.hpp>
#include <opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
@@ -170,7 +171,7 @@ private:
std::shared_ptr<GasOilTwoPhaseHystParams> getGasOilParams();
std::shared_ptr<OilWaterTwoPhaseHystParams> getOilWaterParams();
std::shared_ptr<GasWaterTwoPhaseHystParams> getGasWaterParams();
void setConfig();
void setConfig(unsigned satRegionIdx);
void setDrainageParamsOilGas(unsigned elemIdx, unsigned satRegionIdx);
void setDrainageParamsOilWater(unsigned elemIdx, unsigned satRegionIdx);
void setDrainageParamsGasWater(unsigned elemIdx, unsigned satRegionIdx);
@@ -375,6 +376,8 @@ private:
bool enableEndPointScaling_;
std::shared_ptr<EclHysteresisConfig> hysteresisConfig_;
std::vector<std::shared_ptr<WagHysteresisConfig::WagHysteresisConfigRecord>> wagHystersisConfig_;
std::shared_ptr<EclEpsConfig> oilWaterEclEpsConfig_;
std::vector<EclEpsScalingPointsInfo<Scalar>> unscaledEpsInfo_;

2
src/opm/material/fluidmatrixinteractions/EclHysteresisConfig.cpp
View File

@@ -41,6 +41,8 @@ void EclHysteresisConfig::initFromState(const Runspec& runspec)
pcHysteresisModel_ = runspec.hysterPar().pcHysteresisModel();
modParamTrapped_ = runspec.hysterPar().modParamTrapped();
curvatureCapPrs_ = runspec.hysterPar().curvatureCapPrs();
enableWagHyst_ = runspec.hysterPar().activeWag();
}
} // namespace Opm

11
src/opm/material/fluidmatrixinteractions/EclMaterialLawManager.cpp
View File

@@ -101,6 +101,17 @@ initFromState(const EclipseState& eclState)
this->unscaledEpsInfo_[satRegionIdx]
.extractUnscaled(rtep, rfunc, satRegionIdx);
}
// WAG hysteresis parameters per SATNUM.
if (eclState.runspec().hysterPar().activeWag()) {
if (numSatRegions != eclState.getWagHysteresis().size())
throw std::runtime_error("Inconsistent Wag-hysteresis data");
wagHystersisConfig_.resize(numSatRegions);
for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) {
wagHystersisConfig_[satRegionIdx] = std::make_shared<WagHysteresisConfig::
WagHysteresisConfigRecord >(eclState.getWagHysteresis()[satRegionIdx]);
}
}
}
template<class TraitsT>

9
src/opm/material/fluidmatrixinteractions/EclMaterialLawManagerHystParams.cpp
View File

@@ -78,11 +78,18 @@ getGasWaterParams()
template <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::HystParams::
setConfig()
setConfig(unsigned satRegionIdx)
{
this->gasOilParams_->setConfig(this->parent_.hysteresisConfig_);
this->oilWaterParams_->setConfig(this->parent_.hysteresisConfig_);
this->gasWaterParams_->setConfig(this->parent_.hysteresisConfig_);
if (this->parent_.hysteresisConfig_->enableWagHysteresis()) {
this->gasOilParams_->setWagConfig(this->parent_.wagHystersisConfig_[satRegionIdx]);
this->oilWaterParams_->setWagConfig(this->parent_.wagHystersisConfig_[satRegionIdx]);
this->gasWaterParams_->setWagConfig(this->parent_.wagHystersisConfig_[satRegionIdx]);
}
} // namespace Opm
template <class Traits>

2
src/opm/material/fluidmatrixinteractions/EclMaterialLawManagerInitParams.cpp
View File

@@ -65,7 +65,7 @@ run() {
unsigned satRegionIdx = satRegion_(*satnumArray[i], elemIdx);
//unsigned satNumCell = this->parent_.satnumRegionArray_[elemIdx];
HystParams hystParams {*this};
hystParams.setConfig();
hystParams.setConfig(satRegionIdx);
hystParams.setDrainageParamsOilGas(elemIdx, satRegionIdx);
hystParams.setDrainageParamsOilWater(elemIdx, satRegionIdx);
hystParams.setDrainageParamsGasWater(elemIdx, satRegionIdx);