diff --git a/opm/material/fluidmatrixinteractions/EclStone2Material.hpp b/opm/material/fluidmatrixinteractions/EclStone2Material.hpp new file mode 100644 index 000000000..fbd910efb --- /dev/null +++ b/opm/material/fluidmatrixinteractions/EclStone2Material.hpp @@ -0,0 +1,326 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright (C) 2008-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 . +*/ +/*! + * \file + * \copydoc Opm::EclStone2Material + */ +#ifndef OPM_ECL_STONE2_MATERIAL_HPP +#define OPM_ECL_STONE2_MATERIAL_HPP + +#include "EclStone2MaterialParams.hpp" + +#include +#include + +#include +#include + +#include + +namespace Opm { + +/*! + * \ingroup material + * + * \brief Implements the second phase capillary pressure/relperm law suggested by Stone + * as used by the ECLipse simulator. + * + * This material law is valid for three fluid phases and only depends + * on the saturations. + * + * The required two-phase relations are supplied by means of template + * arguments and can be an arbitrary other material laws. (Provided + * that these only depend on saturation.) + */ +template > +class EclStone2Material : public TraitsT +{ +public: + typedef GasOilMaterialLawT GasOilMaterialLaw; + typedef OilWaterMaterialLawT OilWaterMaterialLaw; + + // some safety checks + static_assert(TraitsT::numPhases == 3, + "The number of phases considered by this capillary pressure " + "law is always three!"); + static_assert(GasOilMaterialLaw::numPhases == 2, + "The number of phases considered by the gas-oil capillary " + "pressure law must be two!"); + static_assert(OilWaterMaterialLaw::numPhases == 2, + "The number of phases considered by the oil-water capillary " + "pressure law must be two!"); + static_assert(std::is_same::value, + "The two two-phase capillary pressure laws must use the same " + "type of floating point values."); + + static_assert(GasOilMaterialLaw::implementsTwoPhaseSatApi, + "The gas-oil material law must implement the two-phase saturation " + "only API to for the default Ecl capillary pressure law!"); + static_assert(OilWaterMaterialLaw::implementsTwoPhaseSatApi, + "The oil-water material law must implement the two-phase saturation " + "only API to for the default Ecl capillary pressure law!"); + + typedef TraitsT Traits; + typedef ParamsT Params; + typedef typename Traits::Scalar Scalar; + + static const int numPhases = 3; + static const int waterPhaseIdx = Traits::wettingPhaseIdx; + static const int oilPhaseIdx = Traits::nonWettingPhaseIdx; + static const int gasPhaseIdx = Traits::gasPhaseIdx; + + //! Specify whether this material law implements the two-phase + //! convenience API + static const bool implementsTwoPhaseApi = false; + + //! Specify whether this material law implements the two-phase + //! convenience API which only depends on the phase saturations + static const bool implementsTwoPhaseSatApi = false; + + //! 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 Implements the default three phase capillary pressure law + * used by the ECLipse simulator. + * + * This material law is valid for three fluid phases and only + * depends on the saturations. + * + * The required two-phase relations are supplied by means of template + * arguments and can be an arbitrary other material laws. + * + * \param values Container for the return values + * \param params Parameters + * \param state The fluid state + */ + template + static void capillaryPressures(ContainerT &values, + const Params ¶ms, + const FluidState &state) + { + typedef typename std::remove_reference::type Evaluation; + values[gasPhaseIdx] = pcgn(params, state); + values[oilPhaseIdx] = 0; + values[waterPhaseIdx] = - pcnw(params, state); + Valgrind::CheckDefined(values[gasPhaseIdx]); + Valgrind::CheckDefined(values[oilPhaseIdx]); + Valgrind::CheckDefined(values[waterPhaseIdx]); + } + + /*! + * \brief Capillary pressure between the gas and the non-wetting + * liquid (i.e., oil) phase. + * + * This is defined as + * \f[ + * p_{c,gn} = p_g - p_n + * \f] + */ + template + static Evaluation pcgn(const Params ¶ms, + const FluidState &fs) + { + typedef MathToolbox FsToolbox; + + const auto& Sw = 1.0 - FsToolbox::template toLhs(fs.saturation(gasPhaseIdx)); + return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw); + } + + /*! + * \brief Capillary pressure between the non-wetting liquid (i.e., + * oil) and the wetting liquid (i.e., water) phase. + * + * This is defined as + * \f[ + * p_{c,nw} = p_n - p_w + * \f] + */ + template + static Evaluation pcnw(const Params ¶ms, + const FluidState &fs) + { + typedef MathToolbox FsToolbox; + + const auto& Sw = FsToolbox::template toLhs(fs.saturation(waterPhaseIdx)); + Valgrind::CheckDefined(Sw); + const auto& result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw); + Valgrind::CheckDefined(result); + return result; + } + + /*! + * \brief The inverse of the capillary pressure + */ + template + static void saturations(ContainerT &values, + const Params ¶ms, + const FluidState &fs) + { + OPM_THROW(std::logic_error, "Not implemented: saturations()"); + } + + /*! + * \brief The saturation of the gas phase. + */ + template + static Evaluation Sg(const Params ¶ms, + const FluidState &fluidState) + { + OPM_THROW(std::logic_error, "Not implemented: Sg()"); + } + + /*! + * \brief The saturation of the non-wetting (i.e., oil) phase. + */ + template + static Evaluation Sn(const Params ¶ms, + const FluidState &fluidState) + { + OPM_THROW(std::logic_error, "Not implemented: Sn()"); + } + + /*! + * \brief The saturation of the wetting (i.e., water) phase. + */ + template + static Evaluation Sw(const Params ¶ms, + const FluidState &fluidState) + { + OPM_THROW(std::logic_error, "Not implemented: Sw()"); + } + + /*! + * \brief The relative permeability of all phases. + * + * The relative permeability of the water phase it uses the same + * value as the relative permeability for water in the water-oil + * law with \f$S_o = 1 - S_w\f$. The gas relative permebility is + * taken from the gas-oil material law, but with \f$S_o = 1 - + * S_g\f$. The relative permeability of the oil phase is + * calculated using the relative permeabilities of the oil phase + * in the two two-phase systems. + * + * A more detailed description can be found in the "Three phase + * oil relative permeability models" section of the ECLipse + * technical description. + */ + template + static void relativePermeabilities(ContainerT &values, + const Params ¶ms, + const FluidState &fluidState) + { + typedef typename std::remove_reference::type Evaluation; + + values[waterPhaseIdx] = krw(params, fluidState); + values[oilPhaseIdx] = krn(params, fluidState); + values[gasPhaseIdx] = krg(params, fluidState); + } + + /*! + * \brief The relative permeability of the gas phase. + */ + template + static Evaluation krg(const Params ¶ms, + const FluidState &fluidState) + { + typedef MathToolbox FsToolbox; + + const Evaluation& Sw = 1 - FsToolbox::template toLhs(fluidState.saturation(gasPhaseIdx)); + return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw); + } + + /*! + * \brief The relative permeability of the wetting phase. + */ + template + static Evaluation krw(const Params ¶ms, + const FluidState &fluidState) + { + typedef MathToolbox FsToolbox; + + const Evaluation& Sw = FsToolbox::template toLhs(fluidState.saturation(waterPhaseIdx)); + return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw); + } + + /*! + * \brief The relative permeability of the non-wetting (i.e., oil) phase. + */ + template + static Evaluation krn(const Params ¶ms, + const FluidState &fluidState) + { + typedef MathToolbox FsToolbox; + + Scalar Swco = params.Swl(); + const Evaluation& Sw = FsToolbox::template toLhs(fluidState.saturation(waterPhaseIdx)); + const Evaluation& Sg = FsToolbox::template toLhs(fluidState.saturation(gasPhaseIdx)); + + Scalar krocw = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Swco); + Evaluation krow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw); + Evaluation krw = OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw); + Evaluation krg = GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), 1 - Sg); + Evaluation krog = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg); + + return krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg); + } + + /*! + * \brief Update the hysteresis parameters after a time step. + * + * This assumes that the nested two-phase material laws are parameters for + * EclHysteresisLaw. If they are not, calling this methid will cause a compiler + * error. (But not calling it will still work.) + */ + template + static void updateHysteresis(Params ¶ms, const FluidState &fluidState) + { + typedef MathToolbox FsToolbox; + + Scalar Sw = FsToolbox::value(fluidState.saturation(waterPhaseIdx)); + Scalar Sg = FsToolbox::value(fluidState.saturation(gasPhaseIdx)); + + params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/Sw); + params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krnSw=*/1 - Sg); + } +}; +} // namespace Opm + +#endif diff --git a/opm/material/fluidmatrixinteractions/EclStone2MaterialParams.hpp b/opm/material/fluidmatrixinteractions/EclStone2MaterialParams.hpp new file mode 100644 index 000000000..18d764494 --- /dev/null +++ b/opm/material/fluidmatrixinteractions/EclStone2MaterialParams.hpp @@ -0,0 +1,145 @@ +// -*- 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 . +*/ +/*! + * \file + * \copydoc Opm::EclStone2MaterialParams + */ +#ifndef OPM_ECL_STONE2_MATERIAL_PARAMS_HPP +#define OPM_ECL_STONE2_MATERIAL_PARAMS_HPP + +#include +#include +#include + +namespace Opm { + +/*! + * \brief Default implementation for the parameters required by the + * three-phase capillary pressure/relperm Stone 2 model used by + * Eclipse. + * + * Essentially, this class just stores the two parameter objects for + * the twophase capillary pressure laws. + */ +template +class EclStone2MaterialParams +{ + typedef typename Traits::Scalar Scalar; + enum { numPhases = 3 }; +public: + typedef GasOilParamsT GasOilParams; + typedef OilWaterParamsT OilWaterParams; + + /*! + * \brief The default constructor. + */ + EclStone2MaterialParams() + { +#ifndef NDEBUG + finalized_ = false; +#endif + } + + /*! + * \brief Finish the initialization of the parameter object. + */ + void finalize() + { +#ifndef NDEBUG + finalized_ = true; +#endif + } + + /*! + * \brief The parameter object for the gas-oil twophase law. + */ + const GasOilParams& gasOilParams() const + { assertFinalized_(); return *gasOilParams_; } + + /*! + * \brief The parameter object for the gas-oil twophase law. + */ + GasOilParams& gasOilParams() + { assertFinalized_(); return *gasOilParams_; } + + /*! + * \brief Set the parameter object for the gas-oil twophase law. + */ + void setGasOilParams(std::shared_ptr val) + { gasOilParams_ = val; } + + /*! + * \brief The parameter object for the oil-water twophase law. + */ + const OilWaterParams& oilWaterParams() const + { assertFinalized_(); return *oilWaterParams_; } + + /*! + * \brief The parameter object for the oil-water twophase law. + */ + OilWaterParams& oilWaterParams() + { assertFinalized_(); return *oilWaterParams_; } + + /*! + * \brief Set the parameter object for the oil-water twophase law. + */ + void setOilWaterParams(std::shared_ptr val) + { oilWaterParams_ = val; } + + /*! + * \brief Set the saturation of "connate" water. + * + * According to + * + * http://www.glossary.oilfield.slb.com/en/Terms/c/connate_water.aspx + * + * the connate water is the water which is trapped in the pores of the rock during + * the rock's formation. For our application, this is basically a reduction of the + * rock's porosity... + */ + void setSwl(Scalar val) + { Swl_ = val; } + + /*! + * \brief Return the saturation of "connate" water. + */ + Scalar Swl() const + { assertFinalized_(); return Swl_; } + +private: +#ifndef NDEBUG + void assertFinalized_() const + { assert(finalized_); } + + bool finalized_; +#else + void assertFinalized_() const + { } +#endif + + std::shared_ptr gasOilParams_; + std::shared_ptr oilWaterParams_; + + Scalar Swl_; +}; +} // namespace Opm + +#endif diff --git a/tests/test_fluidmatrixinteractions.cpp b/tests/test_fluidmatrixinteractions.cpp index 003f8745d..c06785063 100644 --- a/tests/test_fluidmatrixinteractions.cpp +++ b/tests/test_fluidmatrixinteractions.cpp @@ -46,6 +46,7 @@ #include #include #include +#include // include the helper classes to construct traits #include @@ -333,6 +334,15 @@ int main(int argc, char **argv) testThreePhaseApi(); //testThreePhaseSatApi(); } + { + typedef Opm::BrooksCorey TwoPhaseMaterial; + typedef Opm::EclStone2Material MaterialLaw; + testGenericApi(); + testThreePhaseApi(); + //testThreePhaseSatApi(); + } { typedef Opm::ThreePhaseParkerVanGenuchten MaterialLaw; testGenericApi();