fluid-matrix interactions: add a tabulated two-phase pc/kr law which uses splines

compared to the variant which uses linear interpolation, this makes
the resulting curves smoother and thus hopefully improves the
convergence rates of the simulations.
This commit is contained in:
Andreas Lauser 2014-08-06 13:49:25 +02:00
parent 1108f24c33
commit 34f50e6613
5 changed files with 512 additions and 10 deletions

View File

@ -87,7 +87,6 @@ public:
static const int numPhases = 3;
static const int waterPhaseIdx = Traits::wettingPhaseIdx;
static const int nonWettingPhaseIdx = Traits::nonWettingPhaseIdx;
static const int oilPhaseIdx = Traits::nonWettingPhaseIdx;
static const int gasPhaseIdx = Traits::gasPhaseIdx;
@ -235,7 +234,7 @@ public:
const FluidState &fluidState)
{
values[waterPhaseIdx] = krw(params, fluidState);
values[nonWettingPhaseIdx] = krn(params, fluidState);
values[oilPhaseIdx] = krn(params, fluidState);
values[gasPhaseIdx] = krg(params, fluidState);
}

View File

@ -23,19 +23,14 @@
#ifndef OPM_PIECEWISE_LINEAR_TWO_PHASE_MATERIAL_PARAMS_HPP
#define OPM_PIECEWISE_LINEAR_TWO_PHASE_MATERIAL_PARAMS_HPP
#include <opm/parser/eclipse/Utility/SwofTable.hpp>
#include <opm/parser/eclipse/Utility/SgofTable.hpp>
#include <vector>
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief Specification of the material parameters for the van
* Genuchten constitutive relations.
*
* In this implementation setting either the \f$n\f$ or \f$m\f$ shape
* parameter automatically calculates the other. I.e. they cannot be
* set independently.
* \brief Specification of the material parameters for a two-phase material law which
* uses a table and piecewise constant interpolation.
*/
template<class TraitsT>
class PiecewiseLinearTwoPhaseMaterialParams

View File

@ -0,0 +1,349 @@
/*
Copyright (C) 2009-2013 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::SplineTwoPhaseMaterial
*/
#ifndef OPM_SPLINE_TWO_PHASE_MATERIAL_HPP
#define OPM_SPLINE_TWO_PHASE_MATERIAL_HPP
#include "SplineTwoPhaseMaterialParams.hpp"
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Exceptions.hpp>
#include <algorithm>
#include <cmath>
#include <cassert>
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief Implementation of a tabulated capillary pressure and relperm law which uses
* spline curves as interpolation functions.
*/
template <class TraitsT, class ParamsT = SplineTwoPhaseMaterialParams<TraitsT> >
class SplineTwoPhaseMaterial : public TraitsT
{
typedef typename ParamsT::SamplePoints SamplePoints;
public:
//! The traits class for this material law
typedef TraitsT Traits;
//! The type of the parameter objects for this law
typedef ParamsT Params;
//! The type of the scalar values for this law
typedef typename Traits::Scalar Scalar;
//! The number of fluid phases
static const int numPhases = Traits::numPhases;
static_assert(numPhases == 2,
"The piecewise linear two-phase capillary pressure law only"
"applies to the case of two fluid phases");
//! Specify whether this material law implements the two-phase
//! convenience API
static const bool implementsTwoPhaseApi = true;
//! Specify whether this material law implements the two-phase
//! convenience API which only depends on the phase saturations
static const bool implementsTwoPhaseSatApi = true;
//! 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 curve.
*/
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
}
/*!
* \brief The saturations of the fluid phases starting from their
* pressure differences.
*/
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{ OPM_THROW(std::logic_error, "Not implemented: saturations()"); }
/*!
* \brief The relative permeabilities
*/
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
}
/*!
* \brief The derivative of all capillary pressures in regard to
* a given phase saturation.
*/
template <class ContainerT, class FluidState>
static void dCapillaryPressures_dSaturation(ContainerT &values,
const Params &params,
const FluidState &state,
int satPhaseIdx)
{
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx) {
values[Traits::nonWettingPhaseIdx] =
params.pcnwSpline().evalDerivative(state.saturation(Traits::wettingPhaseIdx),
/*extrapolate=*/true);
}
}
/*!
* \brief The derivative of all capillary pressures in regard to
* a given phase pressure.
*/
template <class ContainerT, class FluidState>
static void dCapillaryPressures_dPressure(ContainerT &values,
const Params &params,
const FluidState &state,
int pPhaseIdx)
{
// -> not pressure dependent
for (int pcPhaseIdx = 0; pcPhaseIdx < numPhases; ++pcPhaseIdx)
values[pcPhaseIdx] = 0.0;
}
/*!
* \brief The derivative of all capillary pressures in regard to
* temperature.
*/
template <class ContainerT, class FluidState>
static void dCapillaryPressures_dTemperature(ContainerT &values,
const Params &params,
const FluidState &state)
{
// -> not temperature dependent
for (int pcPhaseIdx = 0; pcPhaseIdx < numPhases; ++pcPhaseIdx)
values[pcPhaseIdx] = 0.0;
}
/*!
* \brief The derivative of all capillary pressures in regard to
* a given mole fraction of a component in a phase.
*/
template <class ContainerT, class FluidState>
static void dCapillaryPressures_dMoleFraction(ContainerT &values,
const Params &params,
const FluidState &state,
int phaseIdx,
int compIdx)
{
// -> not composition dependent
for (int pcPhaseIdx = 0; pcPhaseIdx < numPhases; ++pcPhaseIdx)
values[pcPhaseIdx] = 0.0;
}
/*!
* \brief The derivative of all relative permeabilities in regard to
* a given phase saturation.
*/
template <class ContainerT, class FluidState>
static void dRelativePermeabilities_dSaturation(ContainerT &values,
const Params &params,
const FluidState &state,
int satPhaseIdx)
{
if (satPhaseIdx == Traits::wettingPhaseIdx) {
values[Traits::wettingPhaseIdx] =
params.pcnwSpline().evalDerivative(state.saturation(Traits::wettingPhaseIdx),
/*extrapolate=*/true);
values[Traits::nonWettingPhaseIdx] = 0;
}
else {
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] =
- params.pcnwSpline().evalDerivative(1 - state.saturation(Traits::nonWettingPhaseIdx),
/*extrapolate=*/true);
}
}
/*!
* \brief The derivative of all relative permeabilities in regard to
* a given phase pressure.
*/
template <class ContainerT, class FluidState>
static void dRelativePermeabilities_dPressure(ContainerT &values,
const Params &params,
const FluidState &state,
int pPhaseIdx)
{
// -> not pressure dependent
for (int krPhaseIdx = 0; krPhaseIdx < numPhases; ++krPhaseIdx)
values[krPhaseIdx] = 0.0;
}
/*!
* \brief The derivative of all relative permeabilities in regard to
* temperature.
*/
template <class ContainerT, class FluidState>
static void dRelativePermeabilities_dTemperature(ContainerT &values,
const Params &params,
const FluidState &state)
{
// -> not temperature dependent
for (int krPhaseIdx = 0; krPhaseIdx < numPhases; ++krPhaseIdx)
values[krPhaseIdx] = 0.0;
}
/*!
* \brief The derivative of all relative permeabilities in regard to
* a given mole fraction of a component in a phase.
*/
template <class ContainerT, class FluidState>
static void dRelativePermeabilities_dMoleFraction(ContainerT &values,
const Params &params,
const FluidState &state,
int phaseIdx,
int compIdx)
{
// -> not composition dependent
for (int krPhaseIdx = 0; krPhaseIdx < numPhases; ++krPhaseIdx)
values[krPhaseIdx] = 0.0;
}
/*!
* \brief The capillary pressure-saturation curve
*/
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
return twoPhaseSatPcnw(params, Sw);
}
/*!
* \brief The saturation-capillary pressure curve
*/
static Scalar twoPhaseSatPcnw(const Params &params, Scalar Sw)
{ return params.pcnwSpline().eval(Sw, /*extrapolate=*/true); }
/*!
* \brief The saturation-capillary pressure curve
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{ OPM_THROW(std::logic_error, "Not implemented: Sw()"); }
static Scalar twoPhaseSatSw(const Params &params, Scalar pC)
{ OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatSw()"); }
/*!
* \brief Calculate the non-wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw(params, fs); }
static Scalar twoPhaseSatSn(const Params &params, Scalar pC)
{ return 1 - twoPhaseSatSw(params, pC); }
/*!
* \brief The partial derivative of the capillary pressure with
* regard to the saturation
*/
template <class FluidState>
static Scalar dPcnw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDPcnw_dSw(const Params &params, Scalar Sw)
{
assert(0 < Sw && Sw < 1);
return params.pcnwSpline().evalDerivative(Sw, /*extrapolate=*/true);
}
/*!
* \brief The relative permeability for the wetting phase of the
* porous medium
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{ return params.krwSpline().eval(Sw, /*extrapolate=*/true); }
/*!
* \brief The derivative of the relative permeability of the
* wetting phase in regard to the wetting saturation of the
* porous medium
*/
template <class FluidState>
static Scalar dKrw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDkrw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDKrw_dSw(const Params &params, Scalar Sw)
{ return params.krwSpline().evalDerivative(Sw, /*extrapolate=*/true); }
/*!
* \brief The relative permeability for the non-wetting phase
* of the porous medium
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{ return params.krnSpline().eval(Sw, /*extrapolate=*/true); }
/*!
* \brief The derivative of the relative permeability for the
* non-wetting phase in regard to the wetting saturation of
* the porous medium
*/
template <class FluidState>
static Scalar dKrn_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDkrn_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDKrn_dSw(const Params &params, Scalar Sw)
{ return params.krnSpline().evalDerivative(Sw, /*extrapolate=*/true); }
};
} // namespace Opm
#endif

View File

@ -0,0 +1,152 @@
/*
Copyright (C) 2009-2013 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::SplineTwoPhaseMaterialParams
*/
#ifndef OPM_SPLINE_TWO_PHASE_MATERIAL_PARAMS_HPP
#define OPM_SPLINE_TWO_PHASE_MATERIAL_PARAMS_HPP
#include <opm/core/utility/Spline.hpp>
#include <vector>
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief Specification of the material parameters for a two-phase material law which
* uses a table and spline-based interpolation.
*/
template<class TraitsT>
class SplineTwoPhaseMaterialParams
{
typedef typename TraitsT::Scalar Scalar;
public:
typedef std::vector<Scalar> SamplePoints;
typedef Opm::Spline<Scalar> Spline;
typedef typename Spline::SplineType SplineType;
typedef TraitsT Traits;
SplineTwoPhaseMaterialParams()
{
#ifndef NDEBUG
finalized_ = false;
#endif
}
/*!
* \brief Calculate all dependent quantities once the independent
* quantities of the parameter object have been set.
*/
void finalize()
{
#ifndef NDEBUG
finalized_ = true;
#endif
}
/*!
* \brief Return the sampling points for the capillary pressure curve.
*
* This curve is assumed to depend on the wetting phase saturation
*/
const Spline& pcnwSpline() const
{ assertFinalized_(); return pcwnSpline_; }
/*!
* \brief Set the sampling points for the capillary pressure curve.
*
* This curve is assumed to depend on the wetting phase saturation
*/
void setPcnwSamples(const SamplePoints& SwSamplePoints,
const SamplePoints& pcnwSamplePoints,
SplineType splineType = Spline::Monotonic)
{
assert(SwSamplePoints.size() == pcnwSamplePoints.size());
pcwnSpline_.setXYContainers(SwSamplePoints, pcnwSamplePoints, splineType);
}
/*!
* \brief Return the sampling points for the relative permeability
* curve of the wetting phase.
*
* This curve is assumed to depend on the wetting phase saturation
*/
const Spline& krwSpline() const
{ assertFinalized_(); return krwSpline_; }
/*!
* \brief Set the sampling points for the relative permeability
* curve of the wetting phase.
*
* This curve is assumed to depend on the wetting phase saturation
*/
void setKrwSamples(const SamplePoints& SwSamplePoints,
const SamplePoints& krwSamplePoints,
SplineType splineType = Spline::Monotonic)
{
assert(SwSamplePoints.size() == krwSamplePoints.size());
krwSpline_.setXYContainers(SwSamplePoints, krwSamplePoints, splineType);
}
/*!
* \brief Return the sampling points for the relative permeability
* curve of the non-wetting phase.
*
* This curve is assumed to depend on the wetting phase saturation
*/
const Spline& krnSpline() const
{ assertFinalized_(); return krnSpline_; }
/*!
* \brief Set the sampling points for the relative permeability
* curve of the non-wetting phase.
*
* This curve is assumed to depend on the wetting phase saturation
*/
void setKrnSamples(const SamplePoints& SwSamplePoints,
const SamplePoints& krnSamplePoints,
SplineType splineType = Spline::Monotonic)
{
assert(SwSamplePoints.size() == krnSamplePoints.size());
krnSpline_.setXYContainers(SwSamplePoints, krnSamplePoints, splineType);
}
private:
#ifndef NDEBUG
void assertFinalized_() const
{ assert(finalized_); }
bool finalized_;
#else
void assertFinalized_() const
{ }
#endif
Spline SwSpline_;
Spline pcwnSpline_;
Spline krwSpline_;
Spline krnSpline_;
};
} // namespace Opm
#endif

View File

@ -36,6 +36,7 @@
#include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
#include <opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp>
// include the helper classes to construct traits
@ -335,6 +336,12 @@ int main(int argc, char **argv)
testTwoPhaseApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseSatApi<MaterialLaw, TwoPhaseFluidState>();
}
{
typedef Opm::SplineTwoPhaseMaterial<TwoPhaseTraits> MaterialLaw;
testGenericApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseSatApi<MaterialLaw, TwoPhaseFluidState>();
}
{
typedef Opm::VanGenuchten<TwoPhaseTraits> MaterialLaw;
testGenericApi<MaterialLaw, TwoPhaseFluidState>();