Loading...
Searching...
No Matches
ParkerLenhard.hpp
Go to the documentation of this file.
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef OPM_PARKER_LENHARD_HPP
#define OPM_PARKER_LENHARD_HPP
#include "ParkerLenhardParams.hpp"
#include <opm/material/fluidmatrixinteractions/VanGenuchten.hpp>
#include <algorithm>
#include <cassert>
namespace Opm {
template <class ScalarT>
class PLScanningCurve
{
public:
typedef ScalarT Scalar;
PLScanningCurve(Scalar Swr)
{
loopNum_ = 0;
prev_ = new PLScanningCurve(nullptr, // prev
this, // next
-1, // loop number
Swr, // Sw
1e12, // pcnw
Swr, // SwMic
Swr); // SwMdc
next_ = nullptr;
Sw_ = 1.0;
pcnw_ = 0.0;
SwMic_ = 1.0;
SwMdc_ = 1.0;
}
protected:
PLScanningCurve(PLScanningCurve* prevSC,
PLScanningCurve* nextSC,
int loopN,
Scalar SwReversal,
Scalar pcnwReversal,
Scalar SwMiCurve,
Scalar SwMdCurve)
{
prev_ = prevSC;
next_ = nextSC;
loopNum_ = loopN;
Sw_ = SwReversal;
pcnw_ = pcnwReversal;
SwMic_ = SwMiCurve;
SwMdc_ = SwMdCurve;
}
public:
~PLScanningCurve()
{
if (loopNum_ == 0)
delete prev_;
if (loopNum_ >= 0)
delete next_;
}
PLScanningCurve* prev() const
{ return prev_; }
PLScanningCurve* next() const
{ return next_; }
void setNext(Scalar SwReversal,
Scalar pcnwReversal,
Scalar SwMiCurve,
Scalar SwMdCurve)
{
// if next_ is nullptr, delete does nothing, so
// this is valid!!
delete next_;
next_ = new PLScanningCurve(this, // prev
nullptr, // next
loopNum() + 1,
SwReversal,
pcnwReversal,
SwMiCurve,
SwMdCurve);
}
bool isValidAt_Sw(Scalar SwReversal)
{
if (isImbib())
// for inbibition the given saturation
// must be between the start of the
// current imbibition and the the start
// of the last drainage
return this->Sw() < SwReversal && SwReversal < prev_->Sw();
else
// for drainage the given saturation
// must be between the start of the
// last imbibition and the start
// of the current drainage
return prev_->Sw() < SwReversal && SwReversal < this->Sw();
}
bool isImbib()
{ return loopNum()%2 == 1; }
bool isDrain()
{ return !isImbib(); }
int loopNum()
{ return loopNum_; }
Scalar Sw() const
{ return Sw_; }
Scalar pcnw() const
{ return pcnw_; }
Scalar SwMic()
{ return SwMic_; }
Scalar SwMdc()
{ return SwMdc_; }
private:
PLScanningCurve* prev_;
PLScanningCurve* next_;
int loopNum_;
Scalar Sw_;
Scalar pcnw_;
Scalar SwMdc_;
Scalar SwMic_;
};
template <class TraitsT, class ParamsT = ParkerLenhardParams<TraitsT> >
class ParkerLenhard : public TraitsT
{
public:
typedef TraitsT Traits;
typedef ParamsT Params;
typedef typename Traits::Scalar Scalar;
static const int numPhases = Traits::numPhases;
static_assert(numPhases == 2,
"The Parker-Lenhard capillary pressure law only "
"applies to the case of two fluid phases");
static const bool implementsTwoPhaseApi = true;
static const bool implementsTwoPhaseSatApi = true;
static const bool isSaturationDependent = true;
static const bool isPressureDependent = false;
static const bool isTemperatureDependent = false;
static const bool isCompositionDependent = false;
static_assert(Traits::numPhases == 2,
"The number of fluid phases must be two if you want to use "
"this material law!");
private:
typedef typename ParamsT::VanGenuchten VanGenuchten;
typedef PLScanningCurve<Scalar> ScanningCurve;
public:
static void reset(Params& params)
{
delete params.mdc(); // this will work even if mdc_ == nullptr!
params.setMdc(new ScanningCurve(params.SwrPc()));
params.setCsc(params.mdc());
params.setPisc(nullptr);
params.setCurrentSnr(0.0);
}
template <class FluidState>
static void update(Params& params, const FluidState& fs)
{
Scalar Sw = scalarValue(fs.saturation(Traits::wettingPhaseIdx));
if (Sw > 1 - 1e-5) {
// if the absolute saturation is almost 1,
// it means that we're back to the beginning
reset(params);
return;
}
// find the loop number which corrosponds to the
// given effective saturation
ScanningCurve* curve = findScanningCurve_(params, Sw);
// calculate the apparent saturation on the MIC and MDC
// which yield the same capillary pressure as the
// Sw at the current scanning curve
Scalar pc = pcnw<FluidState, Scalar>(params, fs);
Scalar Sw_mic = VanGenuchten::twoPhaseSatSw(params.micParams(), pc);
Scalar Sw_mdc = VanGenuchten::twoPhaseSatSw(params.mdcParams(), pc);
curve->setNext(Sw, pc, Sw_mic, Sw_mdc);
if (!curve->next())
return;
params.setCsc(curve);
// if we're back on the MDC, we also have a new PISC!
if (params.csc() == params.mdc()) {
params.setPisc(params.mdc()->next());
params.setCurrentSnr(computeCurrentSnr_(params, Sw));
}
}
template <class Container, class FluidState>
static void capillaryPressures(Container& values, const Params& params, const FluidState& fs)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
}
template <class Container, class FluidState>
static void saturations(Container& /*values*/, const Params& /*params*/, const FluidState& /*fs*/)
{ throw std::logic_error("Not implemented: ParkerLenhard::saturations()"); }
template <class Container, class FluidState>
static void relativePermeabilities(Container& values, const Params& params, const FluidState& fs)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
}
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params& params, const FluidState& fs)
{
const Evaluation& Sw =
decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatPcnw(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
{
// calculate the current apparent saturation
ScanningCurve* sc = findScanningCurve_(params, scalarValue(Sw));
// calculate the apparant saturation
const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
// if the apparent saturation exceeds the 'legal' limits,
// we also the underlying material law decide what to do.
if (Sw_app > 1) {
return 0.0; // VanGenuchten::pcnw(params.mdcParams(), Sw_app);
}
// put the apparent saturation into the main imbibition or
// drainage curve
Scalar SwAppCurSC = absoluteToApparentSw_(params, sc->Sw());
Scalar SwAppPrevSC = absoluteToApparentSw_(params, sc->prev()->Sw());
const Evaluation& pos = (Sw_app - SwAppCurSC)/(SwAppPrevSC - SwAppCurSC);
if (sc->isImbib()) {
const Evaluation& SwMic =
pos * (sc->prev()->SwMic() - sc->SwMic()) + sc->SwMic();
return VanGenuchten::twoPhaseSatPcnw(params.micParams(), SwMic);
}
else { // sc->isDrain()
const Evaluation& SwMdc =
pos*(sc->prev()->SwMdc() - sc->SwMdc()) + sc->SwMdc();
return VanGenuchten::twoPhaseSatPcnw(params.mdcParams(), SwMdc);
}
}
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params& /*params*/, const FluidState& /*fs*/)
{ throw std::logic_error("Not implemented: ParkerLenhard::Sw()"); }
template <class Evaluation>
static Evaluation twoPhaseSatSw(const Params& /*params*/, const Evaluation& /*pc*/)
{ throw std::logic_error("Not implemented: ParkerLenhard::twoPhaseSatSw()"); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params& params, const FluidState& fs)
{ return 1 - Sw<FluidState, Evaluation>(params, fs); }
template <class Evaluation>
static Evaluation twoPhaseSatSn(const Params& /*params*/, const Evaluation& /*pc*/)
{ throw std::logic_error("Not implemented: ParkerLenhard::twoPhaseSatSn()"); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params& params, const FluidState& fs)
{
const Evaluation& Sw =
decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
{
// for the effective permeability we only use Land's law and
// the relative permeability of the main drainage curve.
const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
return VanGenuchten::twoPhaseSatKrw(params.mdcParams(), Sw_app);
}
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params& params, const FluidState& fs)
{
const Evaluation& Sw =
decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrn(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
{
// for the effective permeability we only use Land's law and
// the relative permeability of the main drainage curve.
const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
return VanGenuchten::twoPhaseSatKrn(params.mdcParams(), Sw_app);
}
template <class Evaluation>
static Evaluation absoluteToApparentSw_(const Params& params, const Evaluation& Sw)
{
return effectiveToApparentSw_(params, absoluteToEffectiveSw_(params, Sw));
}
private:
template <class Evaluation>
static Evaluation absoluteToEffectiveSw_(const Params& params, const Evaluation& Sw)
{ return (Sw - params.SwrPc())/(1 - params.SwrPc()); }
template <class Evaluation>
static Evaluation effectiveToAbsoluteSw_(const Params& params, const Evaluation& Swe)
{ return Swe*(1 - params.SwrPc()) + params.SwrPc(); }
// return the effctive residual non-wetting saturation, given an
// effective wetting saturation
template <class Evaluation>
static Evaluation computeCurrentSnr_(const Params& params, const Evaluation& Sw)
{
// regularize
if (Sw > 1 - params.Snr())
return 0.0;
if (Sw < params.SwrPc())
return params.Snr();
if (params.Snr() == 0.0)
return 0.0;
// use Land's law
Scalar R = 1.0/params.Snr() - 1;
const Evaluation& curSnr = (1 - Sw)/(1 + R*(1 - Sw));
// the current effective residual non-wetting saturation must
// be smaller than the residual non-wetting saturation
assert(curSnr <= params.Snr());
return curSnr;
}
// returns the trapped effective non-wetting saturation for a
// given wetting phase saturation
template <class Evaluation>
static Evaluation trappedEffectiveSn_(const Params& params, const Evaluation& Sw)
{
const Evaluation& Swe = absoluteToEffectiveSw_(params, Sw);
Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
Scalar Snre = absoluteToEffectiveSw_(params, params.currentSnr());
return Snre*(Swe - SwePisc) / (1 - Snre - SwePisc);
}
// returns the apparent saturation of the wetting phase depending
// on the effective saturation
template <class Evaluation>
static Evaluation effectiveToApparentSw_(const Params& params, const Evaluation& Swe)
{
if (params.pisc() == nullptr ||
Swe <= absoluteToEffectiveSw_(params, params.pisc()->Sw()))
{
// we are on the main drainage curve, i.e. no non-wetting
// fluid is trapped -> apparent saturation == effective
// saturation
return Swe;
}
// we are on a imbibition or drainage curve which is not the
// main drainage curve -> apparent saturation == effective
// saturation + trapped effective saturation
return Swe + trappedEffectiveSn_(params, Swe);
}
// Returns the effective saturation to a given apparent one
template <class Evaluation>
static Evaluation apparentToEffectiveSw_(const Params& params, const Evaluation& Swapp)
{
Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
if (params.pisc() == nullptr || Swapp <= SwePisc) {
// we are on the main drainage curve, i.e.
// no non-wetting fluid is trapped
// -> apparent saturation == effective saturation
return Swapp;
}
Scalar Snre = absoluteToEffectiveSw_(params.currentSnr());
return
(Swapp*(1 - Snre - SwePisc) + Snre*SwePisc)
/(1 - SwePisc);
}
// find the loop on which the an effective
// saturation has to be
static ScanningCurve* findScanningCurve_(const Params& params, Scalar Sw)
{
if (params.pisc() == nullptr || Sw <= params.pisc()->Sw()) {
// we don't have a PISC yet, or the effective
// saturation is smaller than the saturation where the
// PISC begins. In this case are on the MDC
return params.mdc();
}
// if we have a primary imbibition curve, and our current
// effective saturation is higher than the beginning of
// the secondary drainage curve. this means we are on the
// PISC again.
if (params.pisc()->next() == nullptr ||
params.pisc()->next()->Sw() < Sw)
{
return params.pisc();
}
ScanningCurve* curve = params.csc()->next();
while (true) {
assert(curve != params.mdc()->prev());
if (curve->isValidAt_Sw(Sw)) {
return curve;
}
curve = curve->prev();
}
}
};
} // namespace Opm
#endif // PARKER_LENHARD_HPP
ParkerLenhardParams.hpp
Default parameter class for the Parker-Lenhard hysteresis model.
VanGenuchten.hpp
Implementation of the van Genuchten capillary pressure - saturation relation.
Opm::PLScanningCurve
Represents a scanning curve in the Parker-Lenhard hysteresis model.
Definition ParkerLenhard.hpp:48
Opm::PLScanningCurve::next
PLScanningCurve * next() const
Return the next scanning curve, i.e.
Definition ParkerLenhard.hpp:119
Opm::PLScanningCurve::~PLScanningCurve
~PLScanningCurve()
Destructor.
Definition ParkerLenhard.hpp:100
Opm::PLScanningCurve::PLScanningCurve
PLScanningCurve(Scalar Swr)
Constructs main imbibition curve.
Definition ParkerLenhard.hpp:58
Opm::PLScanningCurve::isDrain
bool isDrain()
Returns true iff the scanning curve is a drainage curve.
Definition ParkerLenhard.hpp:181
Opm::PLScanningCurve::setNext
void setNext(Scalar SwReversal, Scalar pcnwReversal, Scalar SwMiCurve, Scalar SwMdCurve)
Set the next scanning curve.
Definition ParkerLenhard.hpp:130
Opm::PLScanningCurve::isImbib
bool isImbib()
Returns true iff the scanning curve is a imbibition curve.
Definition ParkerLenhard.hpp:174
Opm::PLScanningCurve::Sw
Scalar Sw() const
Absolute wetting-phase saturation at the scanning curve's reversal point.
Definition ParkerLenhard.hpp:196
Opm::PLScanningCurve::isValidAt_Sw
bool isValidAt_Sw(Scalar SwReversal)
Returns true iff the given effective saturation Swei is within the scope of the curve,...
Definition ParkerLenhard.hpp:154
Opm::PLScanningCurve::SwMdc
Scalar SwMdc()
Apparent saturation of the last reversal point on the pressure MDC.
Definition ParkerLenhard.hpp:216
Opm::PLScanningCurve::pcnw
Scalar pcnw() const
Capillary pressure at the last reversal point.
Definition ParkerLenhard.hpp:202
Opm::PLScanningCurve::prev
PLScanningCurve * prev() const
Return the previous scanning curve, i.e.
Definition ParkerLenhard.hpp:112
Opm::PLScanningCurve::loopNum
int loopNum()
The loop number of the scanning curve.
Definition ParkerLenhard.hpp:189
Opm::PLScanningCurve::SwMic
Scalar SwMic()
Apparent saturation of the last reversal point on the pressure MIC.
Definition ParkerLenhard.hpp:209
Opm::ParkerLenhard
Implements the Parker-Lenhard twophase p_c-Sw hysteresis model.
Definition ParkerLenhard.hpp:240
Opm::ParkerLenhard::reset
static void reset(Params &params)
Resets the hysteresis model to the initial parameters on the main drainage curve.
Definition ParkerLenhard.hpp:289
Opm::ParkerLenhard::Sw
static Evaluation Sw(const Params &, const FluidState &)
Calculate the wetting phase saturations depending on the phase pressures.
Definition ParkerLenhard.hpp:424
Opm::ParkerLenhard::krn
static Evaluation krn(const Params &params, const FluidState &fs)
The relative permeability for the non-wetting phase of the params.
Definition ParkerLenhard.hpp:470
Opm::ParkerLenhard::implementsTwoPhaseSatApi
static const bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition ParkerLenhard.hpp:258
Opm::ParkerLenhard::krw
static Evaluation krw(const Params &params, const FluidState &fs)
The relative permeability for the wetting phase of the medium.
Definition ParkerLenhard.hpp:448
Opm::ParkerLenhard::update
static void update(Params &params, const FluidState &fs)
Set the current absolute saturation for the current timestep.
Definition ParkerLenhard.hpp:303
Opm::ParkerLenhard::relativePermeabilities
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
Returns the relative permeabilities of the phases dependening on the phase saturations.
Definition ParkerLenhard.hpp:364
Opm::ParkerLenhard::isPressureDependent
static const bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition ParkerLenhard.hpp:266
Opm::ParkerLenhard::isTemperatureDependent
static const bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition ParkerLenhard.hpp:270
Opm::ParkerLenhard::absoluteToApparentSw_
static Evaluation absoluteToApparentSw_(const Params &params, const Evaluation &Sw)
Convert an absolute wetting saturation to an apparent one.
Definition ParkerLenhard.hpp:491
Opm::ParkerLenhard::Sn
static Evaluation Sn(const Params &params, const FluidState &fs)
Calculate the non-wetting phase saturations depending on the phase pressures.
Definition ParkerLenhard.hpp:436
Opm::ParkerLenhard::pcnw
static Evaluation pcnw(const Params &params, const FluidState &fs)
Returns the capillary pressure dependend on the phase saturations.
Definition ParkerLenhard.hpp:377
Opm::ParkerLenhard::numPhases
static const int numPhases
The number of fluid phases.
Definition ParkerLenhard.hpp:247
Opm::ParkerLenhard::isSaturationDependent
static const bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition ParkerLenhard.hpp:262
Opm::ParkerLenhard::implementsTwoPhaseApi
static const bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition ParkerLenhard.hpp:254
Opm::ParkerLenhard::capillaryPressures
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
Returns the capillary pressure dependening on the phase saturations.
Definition ParkerLenhard.hpp:343
Opm::ParkerLenhard::saturations
static void saturations(Container &, const Params &, const FluidState &)
Returns the capillary pressure dependening on the phase saturations.
Definition ParkerLenhard.hpp:356
Opm::ParkerLenhard::isCompositionDependent
static const bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition ParkerLenhard.hpp:274
Opm::VanGenuchten::twoPhaseSatPcnw
static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation &Sw)
The saturation-capillary pressure curve according to van Genuchten using a material law specific API.
Definition VanGenuchten.hpp:194
Opm
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30