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

Merge pull request #169 from andlaus/implement_JFUNC

Browse Source 
implement ECL-like Leverett capillary pressure scaling
This commit is contained in:
Andreas Lauser
2016-10-21 09:47:20 +02:00
committed by GitHub
parent 0a09945fb6 ab2650923d
commit b7dd411c1c
4 changed files with 268 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files

105
opm/material/fluidmatrixinteractions/EclEpsConfig.hpp
View File

@@ -33,6 +33,8 @@
#include <opm/parser/eclipse/Deck/DeckKeyword.hpp>
#include <opm/parser/eclipse/Deck/DeckRecord.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>
#endif
#include <string>
@@ -64,6 +66,7 @@ public:
// by default, do not scale anything
enableSatScaling_ = false;
enablePcScaling_ = false;
enableLeverettScaling_ = false;
enableKrwScaling_ = false;
enableKrnScaling_ = false;
enableThreePointKrSatScaling_ = false;
@@ -131,6 +134,58 @@ public:
bool enablePcScaling() const
{ return enablePcScaling_; }
/*!
* \brief Specify whether the Leverett capillary pressure scaling is enabled.
*
* If setting this to true, Leverett capillary pressure scaling will be used instead
* of the normal capillary pressure scaling and the value of enablePcScaling() will
* not matter anymore.
*/
void setEnableLeverettScaling(bool yesno)
{ enableLeverettScaling_ = yesno; }
/*!
* \brief Specify the conversion factor between deck and SI pressure units
*
* i.e. for METRIC decks this should be set to 10^5, for FIELD decks to 6894.7573,
* and so on. This factor is required to be set if Leverett pressure scaling is
* enabled because the ECL Leverett function implementation "misuses" the pressure
* column of the underlying tables (SWOF, SGOF, SWFN, etc.) for the (dimensionless)
* J-Function. Since opm-parser automatically applies the conversion factor for deck
* to SI pressure values to the column, opm-material needs this value in order to
* restore the value of the J-function from the "pressure" value.
*
* Note that this is a BIG *FAT* *_HACK_*
*/
void setDeckPressureConversionFactor(double pToSiFactor)
{ pToSiFactor_ = pToSiFactor; }
/*!
* \brief Returns the conversion factor between deck and SI pressure units
*
* i.e. for METRIC decks this should be set to 10^5, for FIELD decks to 6894.7573,
* and so on. This factor is required to be set if Leverett pressure scaling is
* enabled because the ECL Leverett function implementation "misuses" the pressure
* column of the underlying tables (SWOF, SGOF, SWFN, etc.) for the (dimensionless)
* J-Function. Since opm-parser automatically applies the conversion factor for deck
* to SI pressure values to the column, opm-material needs this value in order to
* restore the value of the J-function from the "pressure" value.
*
* Note that this is a BIG *FAT* *_HACK_*
*/
double deckPressureConversionFactor() const
{ return pToSiFactor_; }
/*!
* \brief Returns whether the Leverett capillary pressure scaling is enabled.
*
* If this returns true, Leverett capillary pressure scaling will be used instead of
* the normal capillary pressure scaling and the value of enablePcScaling() does not
* matter anymore.
*/
bool enableLeverettScaling() const
{ return enableLeverettScaling_; }
#if HAVE_OPM_PARSER
/*!
* \brief Reads all relevant material parameters form a cell of a parsed ECL deck.
@@ -148,6 +203,7 @@ public:
enableSatScaling_ = false;
enableThreePointKrSatScaling_ = false;
enablePcScaling_ = false;
enableLeverettScaling_ = false;
enableKrwScaling_ = false;
enableKrnScaling_ = false;
return;
@@ -177,16 +233,59 @@ public:
else
enableThreePointKrSatScaling_ = false;
// determine the conversion factor for deck pressures to Pascals
pToSiFactor_ = deck.getActiveUnitSystem().getDimension("Pressure").getSIScaling();
auto& props = eclState.get3DProperties();
// check if we are supposed to scale the Y axis of the capillary pressure
if (twoPhaseSystemType == EclOilWaterSystem)
if (twoPhaseSystemType == EclOilWaterSystem) {
enablePcScaling_ =
props.hasDeckDoubleGridProperty("PCW")
|| props.hasDeckDoubleGridProperty("SWATINIT");
// check if Leverett capillary pressure scaling is requested
if (deck.hasKeyword("JFUNC")) {
const auto& jfuncKeyword = deck.getKeyword("JFUNC");
std::string flagString =
jfuncKeyword.getRecord(0).getItem("FLAG").getTrimmedString(0);
std::transform(flagString.begin(),
flagString.end(),
flagString.begin(),
::toupper);
if (flagString == "BOTH" || flagString == "WATER")
enableLeverettScaling_ = true;
}
if (enablePcScaling_ && enableLeverettScaling_)
OPM_THROW(std::runtime_error,
"Capillary pressure scaling and the Leverett scaling function are "
"mutually exclusive: The deck contains the PCW property and the "
"JFUNC keyword applies to the water phase.");
}
else {
assert(twoPhaseSystemType == EclGasOilSystem);
enablePcScaling_ = props.hasDeckDoubleGridProperty("PCG");
// check if Leverett capillary pressure scaling is requested
if (deck.hasKeyword("JFUNC")) {
const auto& jfuncKeyword = deck.getKeyword("JFUNC");
std::string flagString =
jfuncKeyword.getRecord(0).getItem("FLAG").getTrimmedString(0);
std::transform(flagString.begin(),
flagString.end(),
flagString.begin(),
::toupper);
if (flagString == "BOTH" || flagString == "GAS")
enableLeverettScaling_ = true;
}
if (enablePcScaling_ && enableLeverettScaling_)
OPM_THROW(std::runtime_error,
"Capillary pressure scaling and the Leverett scaling function are "
"mutually exclusive: The deck contains the PCG property and the "
"JFUNC keyword applies to the gas phase.");
}
// check if we are supposed to scale the Y axis of the wetting phase relperm
@@ -208,6 +307,9 @@ public:
#endif
private:
// conversion factor for deck pressure units to SI quantities
double pToSiFactor_;
// enable scaling of the input saturations (i.e., rescale the x-Axis)
bool enableSatScaling_;
@@ -221,6 +323,7 @@ private:
// enable the scaling of the capillary pressure and relative permeability outputs
// (i.e., rescale the y-Axis)
bool enablePcScaling_;
bool enableLeverettScaling_;
bool enableKrwScaling_;
bool enableKrnScaling_;
};

130
opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp
View File

@@ -43,6 +43,8 @@
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#endif
#include <opm/material/common/Means.hpp>
#include <array>
#include <string>
#include <iostream>
@@ -87,6 +89,18 @@ public:
retrieveGridPropertyData_(&krw, eclState, kwPrefix+"KRW");
retrieveGridPropertyData_(&kro, eclState, kwPrefix+"KRO");
retrieveGridPropertyData_(&krg, eclState, kwPrefix+"KRG");
// _may_ be needed to calculate the Leverett capillary pressure scaling factor
const auto& ecl3dProps = eclState.get3DProperties();
poro = &ecl3dProps.getDoubleGridProperty("PORO").getData();
if (ecl3dProps.hasDeckDoubleGridProperty("PERMX")) {
permx = &ecl3dProps.getDoubleGridProperty("PERMX").getData();
permy = permx;
}
if (ecl3dProps.hasDeckDoubleGridProperty("PERMY"))
permy = &ecl3dProps.getDoubleGridProperty("PERMY").getData();
}
#endif
@@ -105,6 +119,9 @@ public:
const DoubleData* krw;
const DoubleData* kro;
const DoubleData* krg;
const DoubleData* poro;
const DoubleData* permx;
const DoubleData* permy;
private:
#if HAVE_OPM_PARSER
@@ -153,6 +170,11 @@ struct EclEpsScalingPointsInfo
Scalar maxPcow; // maximum capillary pressure of the oil-water system
Scalar maxPcgo; // maximum capillary pressure of the gas-oil system
// the Leverett capillary pressure scaling factors. (those only make sense for the
// scaled points, for the unscaled ones they are 1.0.)
Scalar pcowLeverettFactor;
Scalar pcgoLeverettFactor;
// maximum relative permabilities
Scalar maxKrw; // maximum relative permability of water
Scalar maxKrow; // maximum relative permability of oil in the oil-water system
@@ -175,6 +197,8 @@ struct EclEpsScalingPointsInfo
<< " Sogu: " << Sogu << "\n"
<< " maxPcow: " << maxPcow << "\n"
<< " maxPcgo: " << maxPcgo << "\n"
<< " pcowLeverettFactor: " << pcowLeverettFactor << "\n"
<< " pcgoLeverettFactor: " << pcgoLeverettFactor << "\n"
<< " maxKrw: " << maxKrw << "\n"
<< " maxKrg: " << maxKrg << "\n"
<< " maxKrow: " << maxKrow << "\n"
@@ -257,15 +281,20 @@ struct EclEpsScalingPointsInfo
throw std::domain_error("No valid saturation keyword family specified");
}
// there are no "unscaled" Leverett factors, so we just set them to 1.0
pcowLeverettFactor = 1.0;
pcgoLeverettFactor = 1.0;
}
#endif
/*!
* \brief Extract the values of the scaled scaling parameters.
*
* I.e., the values which are "seen" by the physical model.
*/
void extractScaled(const EclEpsGridProperties& epsProperties, unsigned cartesianCellIdx)
void extractScaled(const Opm::Deck& deck,
const Opm::EclipseState& /*eclState*/,
const EclEpsGridProperties& epsProperties,
unsigned cartesianCellIdx)
{
// overwrite the unscaled values with the values for the cell if it is
// explicitly specified by the corresponding keyword.
@@ -287,7 +316,75 @@ struct EclEpsScalingPointsInfo
// quite likely that's wrong!
extractGridPropertyValue_(maxKrow, epsProperties.kro, cartesianCellIdx);
extractGridPropertyValue_(maxKrog, epsProperties.kro, cartesianCellIdx);
// compute the Leverett capillary pressure scaling factors if applicable. note
// that this needs to be done using non-SI units to make it correspond to the
// documentation.
pcowLeverettFactor = 1.0;
pcgoLeverettFactor = 1.0;
if (deck.hasKeyword("JFUNC")) {
const auto& jfuncKeyword = deck.getKeyword("JFUNC");
const auto& jfuncRecord = jfuncKeyword.getRecord(0);
const std::string& jfuncFlag = jfuncRecord.getItem("FLAG").template get<std::string>(0);
const std::string& jfuncDir =
jfuncRecord.getItem("DIRECTION").template get<std::string>(0);
Scalar perm;
if (jfuncDir == "X")
perm =
(*epsProperties.permx)[cartesianCellIdx];
else if (jfuncDir == "Y")
perm =
(*epsProperties.permy)[cartesianCellIdx];
else if (jfuncDir == "XY")
// TODO: verify that this really is the arithmetic mean. (the
// documentation just says that the "average" should be used, IMO the
// harmonic mean would be more appropriate because that's what's usually
// applied when calculating the fluxes.)
perm =
Opm::arithmeticMean((*epsProperties.permx)[cartesianCellIdx],
(*epsProperties.permy)[cartesianCellIdx]);
else
OPM_THROW(std::runtime_error, "Illegal direction indicator for the JFUNC "
"keyword ('"<<jfuncDir<<"')");
// convert permeability from m^2 to mD
perm *= 1.01325e15;
Scalar poro = (*epsProperties.poro)[cartesianCellIdx];
Scalar alpha = jfuncRecord.getItem("ALPHA_FACTOR").template get<double>(0);
Scalar beta = jfuncRecord.getItem("BETA_FACTOR").template get<double>(0);
// the part of the Leverett capillary pressure which does not depend on
// surface tension.
Scalar commonFactor = std::pow(poro, alpha)/std::pow(perm, beta);
// multiply the documented constant by 10^5 because we want the pressures
// in [Pa], not in [bar]
const Scalar Uconst = 0.318316 * 1e5;
// compute the oil-water Leverett factor.
if (jfuncFlag == "WATER" || jfuncFlag == "BOTH") {
// note that we use the surface tension in terms of [dyn/cm]
Scalar gamma =
jfuncRecord.getItem("OW_SURFACE_TENSION").getSIDouble(0)
* (/*dyn/N=*/1e5 * /*m/cm=*/1e-2);
pcowLeverettFactor = commonFactor*gamma*Uconst;
}
// compute the gas-oil Leverett factor.
if (jfuncFlag == "GAS" || jfuncFlag == "BOTH") {
// note that we use the surface tension in terms of [dyn/cm]
Scalar gamma =
jfuncRecord.getItem("GO_SURFACE_TENSION").getSIDouble(0)
* (/*dyn/N=*/1e5 * /*m/cm=*/1e-2);
pcgoLeverettFactor = commonFactor*gamma*Uconst;
}
}
}
#endif
private:
#if HAVE_OPM_PARSER
@@ -538,7 +635,10 @@ public:
saturationKrnPoints_[0] = epsInfo.Swl + epsInfo.Sgl;
}
maxPcnw_ = epsInfo.maxPcow;
if (config.enableLeverettScaling())
maxPcnwOrLeverettFactor_ = epsInfo.pcowLeverettFactor;
else
maxPcnwOrLeverettFactor_ = epsInfo.maxPcow;
maxKrw_ = epsInfo.maxKrw;
maxKrn_ = epsInfo.maxKrow;
}
@@ -574,7 +674,11 @@ public:
saturationKrnPoints_[0] = 1.0 - epsInfo.Sgu;
}
maxPcnw_ = epsInfo.maxPcgo;
if (config.enableLeverettScaling())
maxPcnwOrLeverettFactor_ = epsInfo.pcgoLeverettFactor;
else
maxPcnwOrLeverettFactor_ = epsInfo.maxPcgo;
maxKrw_ = epsInfo.maxKrog;
maxKrn_ = epsInfo.maxKrg;
}
@@ -620,13 +724,25 @@ public:
* \brief Sets the maximum capillary pressure
*/
void setMaxPcnw(Scalar value)
{ maxPcnw_ = value; }
{ maxPcnwOrLeverettFactor_ = value; }
/*!
* \brief Returns the maximum capillary pressure
*/
Scalar maxPcnw() const
{ return maxPcnw_; }
{ return maxPcnwOrLeverettFactor_; }
/*!
* \brief Sets the Leverett scaling factor for capillary pressure
*/
void setLeverettFactor(Scalar value)
{ maxPcnwOrLeverettFactor_ = value; }
/*!
* \brief Returns the Leverett scaling factor for capillary pressure
*/
Scalar leverettFactor() const
{ return maxPcnwOrLeverettFactor_; }
/*!
* \brief Sets the maximum wetting phase relative permeability
@@ -661,7 +777,7 @@ public:
private:
// The the points used for the "y-axis" scaling of capillary pressure
Scalar maxPcnw_;
Scalar maxPcnwOrLeverettFactor_;
// The the points used for the "y-axis" scaling of wetting phase relative permability
Scalar maxKrw_;

35
opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp
View File

@@ -468,21 +468,40 @@ private:
template <class Evaluation>
static Evaluation unscaledToScaledPcnw_(const Params &params, const Evaluation& unscaledPcnw)
{
if (!params.config().enablePcScaling())
return unscaledPcnw;
if (params.config().enableLeverettScaling()) {
// HACK: we need to divide the "unscaled capillary pressure" by the
// conversion factor for deck to SI pressures because this quantity is
// actually *not* a pressure but the (dimensionless) value of the
// J-function. Since opm-parser currently always treads this column as a
// pressure, it multiplies it with the conversion factor and we have to
// revert this here.
Scalar alpha = params.scaledPoints().leverettFactor();
return unscaledPcnw*(alpha/params.config().deckPressureConversionFactor());
}
else if (params.config().enablePcScaling()) {
Scalar alpha = params.scaledPoints().maxPcnw()/params.unscaledPoints().maxPcnw();
return unscaledPcnw*alpha;
}
Scalar alpha = params.scaledPoints().maxPcnw()/params.unscaledPoints().maxPcnw();
return unscaledPcnw*alpha;
return unscaledPcnw;
}
template <class Evaluation>
static Evaluation scaledToUnscaledPcnw_(const Params &params, const Evaluation& scaledPcnw)
{
if (!params.config().enablePcScaling())
return scaledPcnw;
if (params.config().enableLeverettScaling()) {
// HACK: we need to multiply the "scaled capillary pressure" by the conversion
// factor for deck to SI pressures because this quantity is actually *not* a
// pressure but the (dimensionless) value of the J-function.
Scalar alpha = params.scaledPoints().leverettFactor();
return scaledPcnw*(params.config().deckPressureConversionFactor()/alpha);
}
else if (params.config().enablePcScaling()) {
Scalar alpha = params.unscaledPoints().maxPcnw()/params.scaledPoints().maxPcnw();
return scaledPcnw/alpha;
}
Scalar alpha = params.unscaledPoints().maxPcnw()/params.scaledPoints().maxPcnw();
return scaledPcnw/alpha;
return scaledPcnw;
}
/*!

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

@@ -274,7 +274,6 @@ private:
oilWaterEclEpsConfig_-> initFromDeck(deck, eclState, Opm::EclOilWaterSystem);
enableEndPointScaling_ = deck.hasKeyword("ENDSCALE");
}
void readGlobalHysteresisOptions_(const Opm::Deck& deck)
@@ -379,12 +378,16 @@ private:
readGasOilScaledPoints_(gasOilScaledInfoVector,
gasOilScaledPointsVector,
gasOilConfig,
deck,
eclState,
epsGridProperties,
elemIdx,
cartElemIdx);
readOilWaterScaledPoints_(oilWaterScaledEpsInfoDrainage_,
oilWaterScaledEpsPointsDrainage,
oilWaterConfig,
deck,
eclState,
epsGridProperties,
elemIdx,
cartElemIdx);
@@ -393,12 +396,16 @@ private:
readGasOilScaledPoints_(gasOilScaledImbInfoVector,
gasOilScaledImbPointsVector,
gasOilConfig,
deck,
eclState,
epsImbGridProperties,
elemIdx,
cartElemIdx);
readOilWaterScaledPoints_(oilWaterScaledImbInfoVector,
oilWaterScaledImbPointsVector,
oilWaterConfig,
deck,
eclState,
epsImbGridProperties,
elemIdx,
cartElemIdx);
@@ -771,6 +778,8 @@ private:
void readGasOilScaledPoints_(InfoContainer& destInfo,
PointsContainer& destPoints,
std::shared_ptr<EclEpsConfig> config,
const Opm::Deck& deck,
const Opm::EclipseState& eclState,
const EclEpsGridProperties& epsGridProperties,
unsigned elemIdx,
unsigned cartElemIdx)
@@ -778,7 +787,7 @@ private:
unsigned satnumIdx = static_cast<unsigned>((*epsGridProperties.satnum)[cartElemIdx]) - 1; // ECL uses Fortran indices!
destInfo[elemIdx] = std::make_shared<EclEpsScalingPointsInfo<Scalar> >(unscaledEpsInfo_[satnumIdx]);
destInfo[elemIdx]->extractScaled(epsGridProperties, cartElemIdx);
destInfo[elemIdx]->extractScaled(deck, eclState, epsGridProperties, cartElemIdx);
destPoints[elemIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
destPoints[elemIdx]->init(*destInfo[elemIdx], *config, EclGasOilSystem);
@@ -788,6 +797,8 @@ private:
void readOilWaterScaledPoints_(InfoContainer& destInfo,
PointsContainer& destPoints,
std::shared_ptr<EclEpsConfig> config,
const Opm::Deck& deck,
const Opm::EclipseState& eclState,
const EclEpsGridProperties& epsGridProperties,
unsigned elemIdx,
unsigned cartElemIdx)
@@ -795,7 +806,7 @@ private:
unsigned satnumIdx = static_cast<unsigned>((*epsGridProperties.satnum)[cartElemIdx]) - 1; // ECL uses Fortran indices!
destInfo[elemIdx] = std::make_shared<EclEpsScalingPointsInfo<Scalar> >(unscaledEpsInfo_[satnumIdx]);
destInfo[elemIdx]->extractScaled(epsGridProperties, cartElemIdx);
destInfo[elemIdx]->extractScaled(deck, eclState, epsGridProperties, cartElemIdx);
destPoints[elemIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
destPoints[elemIdx]->init(*destInfo[elemIdx], *config, EclOilWaterSystem);