implement ECL hysteresis
so far, this is only for relperm hysteresis. the Eclipse docu seems to be incorrect (or at least vague) for capillary pressure hysteresis...
This commit is contained in:
Andreas Lauser
172
opm/material/fluidmatrixinteractions/EclHysteresisConfig.hpp
Normal file
172
opm/material/fluidmatrixinteractions/EclHysteresisConfig.hpp
Normal file
@@ -0,0 +1,172 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
Copyright (C) 2015 by Andreas Lauser
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
* \copydoc Opm::EclHysteresisConfig
|
||||
*/
|
||||
#ifndef OPM_ECL_HYSTERESIS_CONFIG_HPP
|
||||
#define OPM_ECL_HYSTERESIS_CONFIG_HPP
|
||||
|
||||
#if HAVE_OPM_PARSER
|
||||
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
||||
#endif
|
||||
|
||||
#include <opm/material/common/ErrorMacros.hpp>
|
||||
#include <opm/material/common/Exceptions.hpp>
|
||||
|
||||
#include <string>
|
||||
#include <cassert>
|
||||
#include <algorithm>
|
||||
|
||||
namespace Opm {
|
||||
/*!
|
||||
* \ingroup FluidMatrixInteractions
|
||||
*
|
||||
* \brief Specifies the configuration used by the ECL kr/pC hysteresis code
|
||||
*/
|
||||
class EclHysteresisConfig
|
||||
{
|
||||
public:
|
||||
EclHysteresisConfig()
|
||||
{
|
||||
enableHysteresis_ = false;
|
||||
pcHysteresisModel_ = 0;
|
||||
krHysteresisModel_ = 0;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Specify whether hysteresis is enabled or not.
|
||||
*/
|
||||
void setEnableHysteresis(bool yesno)
|
||||
{ enableHysteresis_ = yesno; };
|
||||
|
||||
/*!
|
||||
* \brief Returns whether hysteresis is enabled.
|
||||
*/
|
||||
bool enableHysteresis() const
|
||||
{ return enableHysteresis_; };
|
||||
|
||||
/*!
|
||||
* \brief Set the type of the hysteresis model which is used for capillary pressure.
|
||||
*
|
||||
* -1: capillary pressure hysteresis is disabled
|
||||
* 0: use the Killough model for capillary pressure hysteresis
|
||||
*/
|
||||
void setPcHysteresisModel(int value)
|
||||
{ pcHysteresisModel_ = value; };
|
||||
|
||||
/*!
|
||||
* \brief Return the type of the hysteresis model which is used for capillary pressure.
|
||||
*
|
||||
* -1: capillary pressure hysteresis is disabled
|
||||
* 0: use the Killough model for capillary pressure hysteresis
|
||||
*/
|
||||
int pcHysteresisModel() const
|
||||
{ return pcHysteresisModel_; };
|
||||
|
||||
/*!
|
||||
* \brief Set the type of the hysteresis model which is used for relative permeability.
|
||||
*
|
||||
* -1: relperm hysteresis is disabled
|
||||
* 0: use the Carlson model for relative permeability hysteresis
|
||||
*/
|
||||
void setKrHysteresisModel(int value)
|
||||
{ krHysteresisModel_ = value; };
|
||||
|
||||
/*!
|
||||
* \brief Return the type of the hysteresis model which is used for relative permeability.
|
||||
*
|
||||
* -1: relperm hysteresis is disabled
|
||||
* 0: use the Carlson model for relative permeability hysteresis
|
||||
*/
|
||||
int krHysteresisModel() const
|
||||
{ return krHysteresisModel_; };
|
||||
|
||||
#if HAVE_OPM_PARSER
|
||||
/*!
|
||||
* \brief Reads all relevant material parameters form a cell of a parsed ECL deck.
|
||||
*
|
||||
* This requires that the opm-parser module is available.
|
||||
*/
|
||||
void initFromDeck(Opm::DeckConstPtr deck)
|
||||
{
|
||||
enableHysteresis_ = false;
|
||||
|
||||
if (!deck->hasKeyword("SATOPTS"))
|
||||
return;
|
||||
|
||||
Opm::DeckItemConstPtr satoptsItem = deck->getKeyword("SATOPTS")->getRecord(0)->getItem(0);
|
||||
for (unsigned i = 0; i < satoptsItem->size(); ++i) {
|
||||
std::string satoptsValue = satoptsItem->getString(0);
|
||||
std::transform(satoptsValue.begin(),
|
||||
satoptsValue.end(),
|
||||
satoptsValue.begin(),
|
||||
::toupper);
|
||||
|
||||
if (satoptsValue == "HYSTER")
|
||||
enableHysteresis_ = true;
|
||||
}
|
||||
|
||||
// check for the (deprecated) HYST keyword
|
||||
if (deck->hasKeyword("HYST"))
|
||||
enableHysteresis_ = true;
|
||||
|
||||
if (!enableHysteresis_)
|
||||
return;
|
||||
|
||||
if (!deck->hasKeyword("EHYSTR"))
|
||||
OPM_THROW(std::runtime_error,
|
||||
"Enabling hysteresis via the HYST parameter for SATOPTS requires the "
|
||||
"presence of the EHYSTR keyword");
|
||||
|
||||
Opm::DeckKeywordConstPtr ehystrKeyword = deck->getKeyword("EHYSTR");
|
||||
if (deck->hasKeyword("NOHYKR"))
|
||||
krHysteresisModel_ = -1;
|
||||
else {
|
||||
krHysteresisModel_ = ehystrKeyword->getRecord(0)->getItem("relative_perm_hyst")->getInt(0);
|
||||
if (krHysteresisModel_ != 0)
|
||||
OPM_THROW(std::runtime_error,
|
||||
"Only the Carlson kr hystersis model (indicated by a 0 on the second item"
|
||||
" of the 'EHYSTR' keyword) is supported");
|
||||
}
|
||||
|
||||
if (deck->hasKeyword("NOHYPC"))
|
||||
pcHysteresisModel_ = -1;
|
||||
else {
|
||||
// if capillary pressure hysteresis is enabled, Eclipse always uses the
|
||||
// Killough model
|
||||
pcHysteresisModel_ = 0;
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
private:
|
||||
// enable hysteresis at all
|
||||
bool enableHysteresis_;
|
||||
|
||||
// the capillary pressure and the relperm hysteresis models to be used
|
||||
int pcHysteresisModel_;
|
||||
int krHysteresisModel_;
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
@@ -0,0 +1,267 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
Copyright (C) 2015 by Andreas Lauser
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
* \copydoc Opm::EclHysteresisTwoPhaseLaw
|
||||
*/
|
||||
#ifndef OPM_ECL_HYSTERESIS_TWO_PHASE_LAW_HPP
|
||||
#define OPM_ECL_HYSTERESIS_TWO_PHASE_LAW_HPP
|
||||
|
||||
#include "EclHysteresisTwoPhaseLawParams.hpp"
|
||||
|
||||
namespace Opm {
|
||||
/*!
|
||||
* \ingroup FluidMatrixInteractions
|
||||
*
|
||||
* \brief This material law implements the hysteresis model of the ECL file format
|
||||
*/
|
||||
template <class EffectiveLawT,
|
||||
class ParamsT = EclHysteresisTwoPhaseLawParams<EffectiveLawT> >
|
||||
class EclHysteresisTwoPhaseLaw : public EffectiveLawT::Traits
|
||||
{
|
||||
public:
|
||||
typedef EffectiveLawT EffectiveLaw;
|
||||
typedef typename EffectiveLaw::Params EffectiveLawParams;
|
||||
|
||||
typedef typename EffectiveLaw::Traits Traits;
|
||||
typedef ParamsT Params;
|
||||
typedef typename EffectiveLaw::Scalar Scalar;
|
||||
|
||||
enum { wettingPhaseIdx = Traits::wettingPhaseIdx };
|
||||
enum { nonWettingPhaseIdx = Traits::nonWettingPhaseIdx };
|
||||
|
||||
//! The number of fluid phases
|
||||
static const int numPhases = EffectiveLaw::numPhases;
|
||||
static_assert(numPhases == 2,
|
||||
"The endpoint scaling applies to the nested twophase laws, not to "
|
||||
"the threephase one!");
|
||||
|
||||
//! Specify whether this material law implements the two-phase
|
||||
//! convenience API
|
||||
static const bool implementsTwoPhaseApi = true;
|
||||
|
||||
static_assert(EffectiveLaw::implementsTwoPhaseApi,
|
||||
"The material laws put into EclEpsTwoPhaseLaw must implement the "
|
||||
"two-phase material law API!");
|
||||
|
||||
//! Specify whether this material law implements the two-phase
|
||||
//! convenience API which only depends on the phase saturations
|
||||
static const bool implementsTwoPhaseSatApi = true;
|
||||
|
||||
static_assert(EffectiveLaw::implementsTwoPhaseSatApi,
|
||||
"The material laws put into EclEpsTwoPhaseLaw must implement the "
|
||||
"two-phase material law saturation API!");
|
||||
|
||||
//! Specify whether the quantities defined by this material law
|
||||
//! are saturation dependent
|
||||
static const bool isSaturationDependent = true;
|
||||
|
||||
//! Specify whether the quantities defined by this material law
|
||||
//! are dependent on the absolute pressure
|
||||
static const bool isPressureDependent = false;
|
||||
|
||||
//! Specify whether the quantities defined by this material law
|
||||
//! are temperature dependent
|
||||
static const bool isTemperatureDependent = false;
|
||||
|
||||
//! Specify whether the quantities defined by this material law
|
||||
//! are dependent on the phase composition
|
||||
static const bool isCompositionDependent = false;
|
||||
|
||||
/*!
|
||||
* \brief The capillary pressure-saturation curves depending on absolute saturations.
|
||||
*
|
||||
* \param values A random access container which stores the
|
||||
* relative pressure of each fluid phase.
|
||||
* \param params The parameter object expressing the coefficients
|
||||
* required by the van Genuchten law.
|
||||
* \param fs The fluid state for which the capillary pressure
|
||||
* ought to be calculated
|
||||
*/
|
||||
template <class Container, class FluidState>
|
||||
static void capillaryPressures(Container &values, const Params ¶ms, const FluidState &fs)
|
||||
{
|
||||
OPM_THROW(NotAvailable,
|
||||
"The capillaryPressures(fs) method is not yet implemented");
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief The relative permeability-saturation curves depending on absolute saturations.
|
||||
*
|
||||
* \param values A random access container which stores the
|
||||
* relative permeability of each fluid phase.
|
||||
* \param params The parameter object expressing the coefficients
|
||||
* required by the van Genuchten law.
|
||||
* \param fs The fluid state for which the relative permeabilities
|
||||
* ought to be calculated
|
||||
*/
|
||||
template <class Container, class FluidState>
|
||||
static void relativePermeabilities(Container &values, const Params ¶ms, const FluidState &fs)
|
||||
{
|
||||
OPM_THROW(NotAvailable,
|
||||
"The pcnw(fs) method is not yet implemented");
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief The capillary pressure-saturation curve.
|
||||
*
|
||||
*
|
||||
* \param params A object that stores the appropriate coefficients
|
||||
* for the respective law.
|
||||
*
|
||||
* \return Capillary pressure [Pa] calculated by specific
|
||||
* constitutive relation (e.g. Brooks & Corey, van
|
||||
* Genuchten, linear...)
|
||||
*/
|
||||
template <class FluidState, class Evaluation = typename FluidState::Scalar>
|
||||
static Evaluation pcnw(const Params ¶ms, const FluidState &fs)
|
||||
{
|
||||
OPM_THROW(NotAvailable,
|
||||
"The pcnw(fs) method is not yet implemented");
|
||||
}
|
||||
|
||||
template <class Evaluation>
|
||||
static Evaluation twoPhaseSatPcnw(const Params ¶ms, const Evaluation& Sw)
|
||||
{
|
||||
// TODO: capillary pressure hysteresis
|
||||
return EffectiveLaw::twoPhaseSatPcnw(params.drainageParams(), Sw);
|
||||
/*
|
||||
if (!params.config().enableHysteresis() || params.config().pcHysteresisModel() < 0)
|
||||
return EffectiveLaw::twoPhaseSatPcnw(params.drainageParams(), Sw);
|
||||
|
||||
if (Sw < params.SwMdc())
|
||||
return EffectiveLaw::twoPhaseSatPcnw(params.drainageParams(), Sw);
|
||||
|
||||
const Evaluation& SwEff = Sw;
|
||||
|
||||
//return EffectiveLaw::twoPhaseSatPcnw(params.imbibitionParams(), SwEff);
|
||||
return EffectiveLaw::twoPhaseSatPcnw(params.drainageParams(), SwEff);
|
||||
*/
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief The saturation-capillary pressure curves.
|
||||
*/
|
||||
template <class Container, class FluidState>
|
||||
static void saturations(Container &values, const Params ¶ms, const FluidState &fs)
|
||||
{
|
||||
OPM_THROW(NotAvailable,
|
||||
"The saturations(fs) method is not yet implemented");
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Calculate wetting liquid phase saturation given that
|
||||
* the rest of the fluid state has been initialized
|
||||
*/
|
||||
template <class FluidState, class Evaluation = typename FluidState::Scalar>
|
||||
static Evaluation Sw(const Params ¶ms, const FluidState &fs)
|
||||
{
|
||||
OPM_THROW(NotAvailable,
|
||||
"The Sw(fs) method is not yet implemented");
|
||||
}
|
||||
|
||||
template <class Evaluation>
|
||||
static Evaluation twoPhaseSatSw(const Params ¶ms, const Evaluation& pc)
|
||||
{
|
||||
OPM_THROW(NotAvailable,
|
||||
"The twoPhaseSatSw(pc) method is not yet implemented");
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Calculate non-wetting liquid phase saturation given that
|
||||
* the rest of the fluid state has been initialized
|
||||
*/
|
||||
template <class FluidState, class Evaluation = typename FluidState::Scalar>
|
||||
static Evaluation Sn(const Params ¶ms, const FluidState &fs)
|
||||
{
|
||||
OPM_THROW(NotAvailable,
|
||||
"The Sn(pc) method is not yet implemented");
|
||||
}
|
||||
|
||||
template <class Evaluation>
|
||||
static Evaluation twoPhaseSatSn(const Params ¶ms, const Evaluation& pc)
|
||||
{
|
||||
OPM_THROW(NotAvailable,
|
||||
"The twoPhaseSatSn(pc) method is not yet implemented");
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief The relative permeability for the wetting phase.
|
||||
*
|
||||
* \param params A container object that is populated with the appropriate coefficients for the respective law.
|
||||
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
|
||||
* is constructed accordingly. Afterwards the values are set there, too.
|
||||
* \return Relative permeability of the wetting phase calculated as implied by EffectiveLaw e.g. Brooks & Corey, van Genuchten, linear... .
|
||||
*
|
||||
*/
|
||||
template <class FluidState, class Evaluation = typename FluidState::Scalar>
|
||||
static Evaluation krw(const Params ¶ms, const FluidState &fs)
|
||||
{
|
||||
OPM_THROW(NotAvailable,
|
||||
"The krw(fs) method is not yet implemented");
|
||||
}
|
||||
|
||||
template <class Evaluation>
|
||||
static Evaluation twoPhaseSatKrw(const Params ¶ms, const Evaluation& Sw)
|
||||
{
|
||||
|
||||
// if no relperm hysteresis is enabled, use the drainage curve
|
||||
if (!params.config().enableHysteresis() || params.config().krHysteresisModel() < 0)
|
||||
return EffectiveLaw::twoPhaseSatKrw(params.drainageParams(), Sw);
|
||||
|
||||
// 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.krwSwMdc())
|
||||
return EffectiveLaw::twoPhaseSatKrw(params.drainageParams(), Sw);
|
||||
|
||||
return EffectiveLaw::twoPhaseSatKrw(params.imbibitionParams(),
|
||||
Sw + params.deltaSwImbKrw());
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief The relative permeability of the non-wetting phase.
|
||||
*/
|
||||
template <class FluidState, class Evaluation = typename FluidState::Scalar>
|
||||
static Evaluation krn(const Params ¶ms, const FluidState &fs)
|
||||
{
|
||||
OPM_THROW(NotAvailable,
|
||||
"The krn(fs) method is not yet implemented");
|
||||
}
|
||||
|
||||
template <class Evaluation>
|
||||
static Evaluation twoPhaseSatKrn(const Params ¶ms, const Evaluation& Sw)
|
||||
{
|
||||
// if no relperm hysteresis is enabled, use the drainage curve
|
||||
if (!params.config().enableHysteresis() || params.config().krHysteresisModel() < 0)
|
||||
return EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Sw);
|
||||
|
||||
// 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())
|
||||
return EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Sw);
|
||||
|
||||
return EffectiveLaw::twoPhaseSatKrn(params.imbibitionParams(),
|
||||
Sw + params.deltaSwImbKrn());
|
||||
}
|
||||
};
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
@@ -0,0 +1,366 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
Copyright (C) 2015 by Andreas Lauser
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
* \copydoc Opm::EclHysteresisTwoPhaseLawParams
|
||||
*/
|
||||
#ifndef OPM_ECL_HYSTERESIS_TWO_PHASE_LAW_PARAMS_HPP
|
||||
#define OPM_ECL_HYSTERESIS_TWO_PHASE_LAW_PARAMS_HPP
|
||||
|
||||
#include "EclHysteresisConfig.hpp"
|
||||
#include "EclEpsScalingPoints.hpp"
|
||||
|
||||
#if HAVE_OPM_PARSER
|
||||
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
||||
#endif
|
||||
|
||||
#include <string>
|
||||
#include <memory>
|
||||
#include <cassert>
|
||||
#include <algorithm>
|
||||
|
||||
namespace Opm {
|
||||
/*!
|
||||
* \ingroup FluidMatrixInteractions
|
||||
*
|
||||
* \brief A default implementation of the parameters for the material law which
|
||||
* implements the ECL relative permeability and capillary pressure hysteresis
|
||||
*/
|
||||
template <class EffLawT>
|
||||
class EclHysteresisTwoPhaseLawParams
|
||||
{
|
||||
typedef typename EffLawT::Params EffLawParams;
|
||||
typedef typename EffLawParams::Traits::Scalar Scalar;
|
||||
|
||||
public:
|
||||
typedef typename EffLawParams::Traits Traits;
|
||||
|
||||
EclHysteresisTwoPhaseLawParams()
|
||||
{
|
||||
pcSwMdc_ = 1.0;
|
||||
krnSwMdc_ = 1.0;
|
||||
krwSwMdc_ = 1.0;
|
||||
|
||||
#ifndef NDEBUG
|
||||
finalized_ = false;
|
||||
#endif
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Calculate all dependent quantities once the independent
|
||||
* quantities of the parameter object have been set.
|
||||
*/
|
||||
void finalize()
|
||||
{
|
||||
if (config().enableHysteresis()) {
|
||||
//C_ = 1.0/(Sncri_ - Sncrd_) + 1.0/(Snmaxd_ - Sncrd_);
|
||||
|
||||
updateDynamicParams_();
|
||||
}
|
||||
|
||||
#ifndef NDEBUG
|
||||
finalized_ = true;
|
||||
#endif
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Set the endpoint scaling configuration object.
|
||||
*/
|
||||
void setConfig(std::shared_ptr<EclHysteresisConfig> value)
|
||||
{ config_ = value; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the endpoint scaling configuration object.
|
||||
*/
|
||||
const EclHysteresisConfig& config() const
|
||||
{ return *config_; }
|
||||
|
||||
/*!
|
||||
* \brief Sets the parameters used for the drainage curve
|
||||
*/
|
||||
void setDrainageParams(std::shared_ptr<EffLawParams> value,
|
||||
const EclEpsScalingPointsInfo<Scalar>& info,
|
||||
EclTwoPhaseSystemType twoPhaseSystem)
|
||||
|
||||
{
|
||||
drainageParams_ = value;
|
||||
|
||||
#if 0
|
||||
if (twoPhaseSystem == EclGasOilSystem) {
|
||||
Sncrd_ = info.Sgcr;
|
||||
Snmaxd_ = 1 - info.Sogcr;
|
||||
}
|
||||
else {
|
||||
assert(twoPhaseSystem == EclOilWaterSystem);
|
||||
Sncrd_ = info.Sowcr;
|
||||
Snmaxd_ = 1 - info.Swcr;
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Sets the parameters used for the drainage curve
|
||||
*/
|
||||
const EffLawParams& drainageParams() const
|
||||
{ return *drainageParams_; }
|
||||
|
||||
/*!
|
||||
* \brief Sets the parameters used for the imbibition curve
|
||||
*/
|
||||
void setImbibitionParams(std::shared_ptr<EffLawParams> value,
|
||||
const EclEpsScalingPointsInfo<Scalar>& info,
|
||||
EclTwoPhaseSystemType twoPhaseSystem)
|
||||
{
|
||||
imbibitionParams_ = value;
|
||||
|
||||
/*
|
||||
if (twoPhaseSystem == EclGasOilSystem) {
|
||||
Sncri_ = info.Sgcr;
|
||||
}
|
||||
else {
|
||||
assert(twoPhaseSystem == EclOilWaterSystem);
|
||||
Sncri_ = info.Sowcr;
|
||||
}
|
||||
*/
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Sets the parameters used for the imbibition curve
|
||||
*/
|
||||
const EffLawParams& imbibitionParams() const
|
||||
{ return *imbibitionParams_; }
|
||||
|
||||
/*!
|
||||
* \brief Set the saturation of the wetting phase where the last switch from the main
|
||||
* drainage curve (MDC) to imbibition happend on the capillary pressure curve.
|
||||
*/
|
||||
void setPcSwMdc(Scalar value)
|
||||
{ pcSwMdc_ = value; };
|
||||
|
||||
/*!
|
||||
* \brief Set the saturation of the wetting phase where the last switch from the main
|
||||
* drainage curve to imbibition happend on the capillary pressure curve.
|
||||
*/
|
||||
Scalar pcSwMdc() const
|
||||
{ return pcSwMdc_; };
|
||||
|
||||
/*!
|
||||
* \brief Set the saturation of the wetting phase where the last switch from the main
|
||||
* drainage curve (MDC) to imbibition happend on the relperm curve for the
|
||||
* wetting phase.
|
||||
*/
|
||||
void setKrwSwMdc(Scalar value)
|
||||
{ krwSwMdc_ = value; };
|
||||
|
||||
/*!
|
||||
* \brief Set the saturation of the wetting phase where the last switch from the main
|
||||
* drainage curve to imbibition happend on the relperm curve for the
|
||||
* wetting phase.
|
||||
*/
|
||||
Scalar krwSwMdc() const
|
||||
{ return krwSwMdc_; };
|
||||
|
||||
/*!
|
||||
* \brief Set the saturation of the wetting phase where the last switch from the main
|
||||
* drainage curve (MDC) to imbibition happend on the relperm curve for the
|
||||
* non-wetting phase.
|
||||
*/
|
||||
void setKrnSwMdc(Scalar value)
|
||||
{ krnSwMdc_ = value; };
|
||||
|
||||
/*!
|
||||
* \brief Set the saturation of the wetting phase where the last switch from the main
|
||||
* drainage curve to imbibition happend on the relperm curve for the
|
||||
* non-wetting phase.
|
||||
*/
|
||||
Scalar krnSwMdc() const
|
||||
{ return krnSwMdc_; };
|
||||
|
||||
/*!
|
||||
* \brief Sets the saturation value which must be added if krw is calculated using
|
||||
* the imbibition curve.
|
||||
*
|
||||
* This means that krw(Sw) = krw_drainage(Sw) if Sw < SwMdc and
|
||||
* krw(Sw) = krw_imbibition(Sw + Sw_shift,krw) else
|
||||
*/
|
||||
void setDeltaSwImbKrw(Scalar value)
|
||||
{ deltaSwImbKrw_ = value; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the saturation value which must be added if krw is calculated using
|
||||
* the imbibition curve.
|
||||
*
|
||||
* This means that krw(Sw) = krw_drainage(Sw) if Sw < SwMdc and
|
||||
* krw(Sw) = krw_imbibition(Sw + Sw_shift,krw) else
|
||||
*/
|
||||
Scalar deltaSwImbKrw() const
|
||||
{ return deltaSwImbKrw_; }
|
||||
|
||||
/*!
|
||||
* \brief Sets the saturation value which must be added if krn is calculated using
|
||||
* the imbibition curve.
|
||||
*
|
||||
* This means that krn(Sw) = krn_drainage(Sw) if Sw < SwMdc and
|
||||
* krn(Sw) = krn_imbibition(Sw + Sw_shift,krn) else
|
||||
*/
|
||||
void setDeltaSwImbKrn(Scalar value)
|
||||
{ deltaSwImbKrn_ = value; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the saturation value which must be added if krn is calculated using
|
||||
* the imbibition curve.
|
||||
*
|
||||
* This means that krn(Sw) = krn_drainage(Sw) if Sw < SwMdc and
|
||||
* krn(Sw) = krn_imbibition(Sw + Sw_shift,krn) else
|
||||
*/
|
||||
Scalar deltaSwImbKrn() const
|
||||
{ return deltaSwImbKrn_; }
|
||||
|
||||
#if 0
|
||||
/*!
|
||||
* \brief Sets the "trapped" non-wetting phase saturation.
|
||||
*
|
||||
* This quantity is used for capillary pressure hysteresis. How it should be
|
||||
* calculated depends on the hysteresis model.
|
||||
*/
|
||||
void setSncrt(Scalar value)
|
||||
{ Sncrt_ = value; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the "trapped" non-wetting phase saturation.
|
||||
*
|
||||
* This quantity is used for capillary pressure hysteresis. How it should be
|
||||
* calculated depends on the hysteresis model.
|
||||
*/
|
||||
Scalar Sncrt() const
|
||||
{ return Sncrt_; }
|
||||
#endif
|
||||
|
||||
/*!
|
||||
* \brief Notify the hysteresis law that a given wetting-phase saturation has been seen
|
||||
*
|
||||
* This updates the scanning curves and the imbibition<->drainage reversal points as
|
||||
* appropriate.
|
||||
*/
|
||||
void update(Scalar pcSw, Scalar krwSw, Scalar krnSw)
|
||||
{
|
||||
bool updateParams = false;
|
||||
if (pcSw < pcSwMdc_) {
|
||||
pcSwMdc_ = pcSw;
|
||||
updateParams = true;
|
||||
}
|
||||
|
||||
if (krwSw < krwSwMdc_) {
|
||||
krwSwMdc_ = krwSw;
|
||||
updateParams = true;
|
||||
}
|
||||
|
||||
if (krnSw < krnSwMdc_) {
|
||||
krnSwMdc_ = krnSw;
|
||||
updateParams = true;
|
||||
}
|
||||
|
||||
if (updateParams)
|
||||
updateDynamicParams_();
|
||||
}
|
||||
|
||||
private:
|
||||
|
||||
#ifndef NDEBUG
|
||||
void assertFinalized_() const
|
||||
{ assert(finalized_); }
|
||||
|
||||
bool finalized_;
|
||||
#else
|
||||
void assertFinalized_() const
|
||||
{ }
|
||||
#endif
|
||||
|
||||
void updateDynamicParams_()
|
||||
{
|
||||
// calculate the saturation deltas for the relative permeabilities
|
||||
Scalar krwMdcDrainage = EffLawT::twoPhaseSatKrw(drainageParams(), krwSwMdc_);
|
||||
Scalar SwKrwMdcImbibition = EffLawT::twoPhaseSatKrwInv(imbibitionParams(), krwMdcDrainage);
|
||||
deltaSwImbKrw_ = SwKrwMdcImbibition - krwSwMdc_;
|
||||
|
||||
Scalar krnMdcDrainage = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_);
|
||||
Scalar SwKrnMdcImbibition = EffLawT::twoPhaseSatKrnInv(imbibitionParams(), krnMdcDrainage);
|
||||
deltaSwImbKrn_ = SwKrnMdcImbibition - krnSwMdc_;
|
||||
|
||||
Scalar pcMdcDrainage = EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_);
|
||||
Scalar SwPcMdcImbibition = EffLawT::twoPhaseSatPcnwInv(imbibitionParams(), pcMdcDrainage);
|
||||
deltaSwImbPc_ = SwPcMdcImbibition - pcSwMdc_;
|
||||
|
||||
// assert(std::abs(EffLawT::twoPhaseSatPcnw(imbibitionParams(), pcSwMdc_ + deltaSwImbPc_)
|
||||
// - EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_)) < 1e-8);
|
||||
assert(std::abs(EffLawT::twoPhaseSatKrn(imbibitionParams(), krnSwMdc_ + deltaSwImbKrn_)
|
||||
- EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_)) < 1e-8);
|
||||
assert(std::abs(EffLawT::twoPhaseSatKrw(imbibitionParams(), krwSwMdc_ + deltaSwImbKrw_)
|
||||
- EffLawT::twoPhaseSatKrw(drainageParams(), krwSwMdc_)) < 1e-8);
|
||||
|
||||
#if 0
|
||||
Scalar Snhy = 1.0 - SwMdc_;
|
||||
|
||||
Sncrt_ = Sncrd_ + (Snhy - Sncrd_)/(1 + C_*(Snhy - Sncrd_));
|
||||
#endif
|
||||
}
|
||||
|
||||
std::shared_ptr<EffLawParams> effectiveLawParams_;
|
||||
|
||||
std::shared_ptr<EclHysteresisConfig> config_;
|
||||
std::shared_ptr<EffLawParams> imbibitionParams_;
|
||||
std::shared_ptr<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 krnSwMdc_;
|
||||
Scalar pcSwMdc_;
|
||||
|
||||
// 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
|
||||
// the capillary pressure
|
||||
Scalar deltaSwImbKrw_;
|
||||
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
|
||||
// Sncri_: critical non-wetting phase saturation for the imbibition curve
|
||||
// Snmaxd_: non-wetting phase saturation where the non-wetting relperm reaches its
|
||||
// maximum on the drainage curve
|
||||
// C_: factor required to calculate the trapped non-wetting phase saturation using
|
||||
// the Killough approach
|
||||
//Scalar Sncrd_;
|
||||
//Scalar Sncri_;
|
||||
//Scalar Snmaxd_;
|
||||
//Scalar C_;
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
@@ -44,6 +44,7 @@
|
||||
#include <opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp>
|
||||
#include <opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp>
|
||||
#include <opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp>
|
||||
#include <opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLaw.hpp>
|
||||
|
||||
// include the helper classes to construct traits
|
||||
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
|
||||
@@ -383,6 +384,14 @@ int main(int argc, char **argv)
|
||||
testTwoPhaseApi<MaterialLaw, TwoPhaseFluidState>();
|
||||
testTwoPhaseSatApi<MaterialLaw, TwoPhaseFluidState>();
|
||||
}
|
||||
{
|
||||
typedef Opm::BrooksCorey<TwoPhaseTraits> RawMaterialLaw;
|
||||
typedef Opm::EclHysteresisTwoPhaseLaw<RawMaterialLaw> MaterialLaw;
|
||||
testGenericApi<MaterialLaw, TwoPhaseFluidState>();
|
||||
testTwoPhaseApi<MaterialLaw, TwoPhaseFluidState>();
|
||||
testTwoPhaseSatApi<MaterialLaw, TwoPhaseFluidState>();
|
||||
}
|
||||
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user