diff --git a/opm/material/fluidmatrixinteractions/3p/3pParkerVanGenuchten.hpp b/opm/material/fluidmatrixinteractions/3p/3pParkerVanGenuchten.hpp deleted file mode 100644 index d9cf2f026..000000000 --- a/opm/material/fluidmatrixinteractions/3p/3pParkerVanGenuchten.hpp +++ /dev/null @@ -1,380 +0,0 @@ -/* - Copyright (C) 2008-2013 by Andreas Lauser - Copyright (C) 2012 by Holger Class - Copyright (C) 2012 by Vishal Jambhekar - - 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::ThreePParkerVanGenuchten - */ -#ifndef OPM_3P_PARKER_VAN_GENUCHTEN_HPP -#define OPM_3P_PARKER_VAN_GENUCHTEN_HPP - -#include "3pParkerVanGenuchtenParams.hpp" - -#include - -namespace Opm { -/*! - * \ingroup FluidMatrixInteractions - * - * \brief Implementation of van Genuchten's capillary pressure <-> - * saturation relation. - * - * \sa VanGenuchten, VanGenuchtenThreephase - */ -template > -class ThreePParkerVanGenuchten -{ - -public: - typedef ParamsT Params; - typedef typename Params::Scalar Scalar; - - /*! - * \brief The capillary pressure-saturation curve. - * - */ - static Scalar pC(const Params ¶ms, Scalar Sw) - { - OPM_THROW(std::logic_error, "Not implemented: Capillary pressures for three phases is not so simple! Use pCGN, pCNW, and pcGW"); - } - - static Scalar pCGW(const Params ¶ms, Scalar Sw) - { - /* - Sw = wetting phase saturation, or, - sum of wetting phase saturations - alpha : VanGenuchten-alpha - this function is copied from MUFTE/pml/constrel3p3cni.c */ - - Scalar r,Se,x,vg_m; - Scalar pc,pc_prime,Se_regu; - Scalar PC_VG_REG = 0.01; - - Se = (Sw-params.Swr())/(1.-params.Sgr()); - - /* Snr = 0.0; test version */ - - /* regularization */ - if (Se<0.0) Se=0.0; - if (Se>1.0) Se=1.0; - vg_m = 1.-1./params.vgN(); - - if (Se>PC_VG_REG && Se<1-PC_VG_REG) - { - r = std::pow(Se,-1/vg_m); - x = r-1; - vg_m = 1-vg_m; - x = std::pow(x,vg_m); - r = x/params.vgAlpha(); - return(r); - } - else - { - /* value and derivative at regularization point */ - if (Se<=PC_VG_REG) Se_regu = PC_VG_REG; else Se_regu = 1-PC_VG_REG; - pc = std::pow(std::pow(Se_regu,-1/vg_m)-1,1/params.vgN())/params.vgAlpha(); - pc_prime = std::pow(std::pow(Se_regu,-1/vg_m)-1,1/params.vgN()-1)*std::pow(Se_regu,-1/vg_m-1)*(-1/vg_m)/params.vgAlpha()/(1-params.Sgr()-params.Swr())/params.vgN(); - - /* evaluate tangential */ - r = (Se-Se_regu)*pc_prime+pc; - return(r/params.betaGW()); - } - } - - static Scalar pCNW(const Params ¶ms, Scalar Sw) - { - /* - Sw = wetting phase saturation, or, - sum of wetting phase saturations - alpha : VanGenuchten-alpha - this function is just copied from MUFTE/pml/constrel3p3cni.c */ - - Scalar r,Se,x,vg_m; - Scalar pc,pc_prime,Se_regu; - Scalar PC_VG_REG = 0.01; - - Se = (Sw-params.Swr())/(1.-params.Snr()); - - /* Snr = 0.0; test version */ - - /* regularization */ - if (Se<0.0) Se=0.0; - if (Se>1.0) Se=1.0; - vg_m = 1.-1./params.vgN(); - - if (Se>PC_VG_REG && Se<1-PC_VG_REG) - { - r = std::pow(Se,-1/vg_m); - x = r-1; - vg_m = 1-vg_m; - x = std::pow(x,vg_m); - r = x/params.vgAlpha(); - return(r); - } - else - { - /* value and derivative at regularization point */ - if (Se<=PC_VG_REG) Se_regu = PC_VG_REG; else Se_regu = 1-PC_VG_REG; - pc = std::pow(std::pow(Se_regu,-1/vg_m)-1,1/params.vgN())/params.vgAlpha(); - pc_prime = std::pow(std::pow(Se_regu,-1/vg_m)-1,1/params.vgN()-1)*std::pow(Se_regu,-1/vg_m-1)*(-1/vg_m)/params.vgAlpha()/(1-params.Snr()-params.Swr())/params.vgN(); - - /* evaluate tangential */ - r = (Se-Se_regu)*pc_prime+pc; - return(r/params.betaNW()); - } - } - - static Scalar pCGN(const Params ¶ms, Scalar St) - { - /* - St = sum of wetting (liquid) phase saturations - alpha : VanGenuchten-alpha - this function is just copied from MUFTE/pml/constrel3p3cni.c */ - - Scalar r,Se,x,vg_m; - Scalar pc,pc_prime,Se_regu; - Scalar PC_VG_REG = 0.01; - - Se = (St-params.Swrx())/(1.-params.Swrx()); - - /* Snr = 0.0; test version */ - - /* regularization */ - if (Se<0.0) Se=0.0; - if (Se>1.0) Se=1.0; - vg_m = 1.-1./params.vgN(); - - if (Se>PC_VG_REG && Se<1-PC_VG_REG) - { - r = std::pow(Se,-1/vg_m); - x = r-1; - vg_m = 1-vg_m; - x = std::pow(x,vg_m); - r = x/params.vgAlpha(); - return(r); - } - else - { - /* value and derivative at regularization point */ - if (Se<=PC_VG_REG) Se_regu = PC_VG_REG; else Se_regu = 1-PC_VG_REG; - pc = std::pow(std::pow(Se_regu,-1/vg_m)-1,1/params.vgN())/params.vgAlpha(); - pc_prime = std::pow(std::pow(Se_regu,-1/vg_m)-1,1/params.vgN()-1)*std::pow(Se_regu,-1/vg_m-1)*(-1/vg_m)/params.vgAlpha()/(1-params.Sgr()-params.Swrx())/params.vgN(); - - /* evaluate tangential */ - r = (Se-Se_regu)*pc_prime+pc; - return(r/params.betaGN()); - } - } - - static Scalar pCAlpha(const Params ¶ms, Scalar Sn) - { - /* continuous transition to zero */ - Scalar alpha,Sne; - - Sne=Sn; - /* regularization */ - if (Sne<=0.001) Sne=0.0; - if (Sne>=1.0) Sne=1.0; - - if (Sne>params.Snr()) alpha = 1.0; - else - { - if (params.Snr()>=0.001) alpha = Sne/params.Snr(); - else alpha = 0.0; - } - return(alpha); - } - - /*! - * \brief The saturation-capillary pressure curve. - * - */ - static Scalar Sw(const Params ¶ms, Scalar pC) - { - OPM_THROW(std::logic_error, "Not implemented: Sw(pc) for three phases not implemented! Do it yourself!"); - } - - /*! - * \brief Returns the partial derivative of the capillary - * pressure to the effective saturation. - * - */ - static Scalar dpC_dSw(const Params ¶ms, Scalar Sw) - { - OPM_THROW(std::logic_error, "Not implemented: dpC/dSw for three phases not implemented! Do it yourself!"); - } - - /*! - * \brief Returns the partial derivative of the effective - * saturation to the capillary pressure. - */ - static Scalar dSw_dpC(const Params ¶ms, Scalar pC) - { - OPM_THROW(std::logic_error, "Not implemented: dSw/dpC for three phases not implemented! Do it yourself!"); - } - - /*! - * \brief The relative permeability for the wetting phase of - * the medium implied by van Genuchten's - * parameterization. - * - * The permeability of water in a 3p system equals the standard 2p description. - * (see p61. in "Comparison of the Three-Phase Oil Relative Permeability Models" - * MOJDEH DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.) - * - * \param Sn saturation of the NAPL phase. - * \param Sg saturation of the gas phase. - * \param saturation saturation of the water phase. - * \param params Array of parameters. - */ - static Scalar krw(const Params ¶ms, Scalar saturation, Scalar Sn, Scalar Sg) - { - - //transformation to effective saturation - Scalar Se = (saturation - params.Swr()) / (1-params.Swr()); - - /* regularization */ - if(Se > 1.0) return 1.; - if(Se < 0.0) return 0.; - - Scalar r = 1. - std::pow(1 - std::pow(Se, 1/params.vgM()), params.vgM()); - return std::sqrt(Se)*r*r; - } - - /*! - * \brief The relative permeability for the non-wetting phase - * after the Model of Parker et al. (1987). - * - * See model 7 in "Comparison of the Three-Phase Oil Relative Permeability Models" - * MOJDEH DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83. - * or more comprehensive in - * "Estimation of primary drainage three-phase relative permeability for organic - * liquid transport in the vadose zone", Leonardo I. Oliveira, Avery H. Demond, - * Journal of Contaminant Hydrology 66 (2003), 261-285 - * - * - * \param Sw saturation of the water phase. - * \param Sg saturation of the gas phase. - * \param saturation saturation of the NAPL phase. - * \param params Array of parameters. - */ - static Scalar krn(const Params ¶ms, Scalar Sw, Scalar saturation, Scalar Sg) - { - - Scalar Swe = std::min((Sw - params.Swr()) / (1 - params.Swr()), 1.); - Scalar Ste = std::min((Sw + saturation - params.Swr()) / (1 - params.Swr()), 1.); - - // regularization - if(Swe <= 0.0) Swe = 0.; - if(Ste <= 0.0) Ste = 0.; - if(Ste - Swe <= 0.0) return 0.; - - Scalar krn_; - krn_ = std::pow(1 - std::pow(Swe, 1/params.vgM()), params.vgM()); - krn_ -= std::pow(1 - std::pow(Ste, 1/params.vgM()), params.vgM()); - krn_ *= krn_; - - if (params.krRegardsSnr()) - { - // regard Snr in the permeability of the n-phase, see Helmig1997 - Scalar resIncluded = std::max(std::min((saturation - params.Snr()/ (1-params.Swr())), 1.), 0.); - krn_ *= std::sqrt(resIncluded ); - } - else - krn_ *= std::sqrt(saturation / (1 - params.Swr())); // Hint: (Ste - Swe) = Sn / (1-Srw) - - - return krn_; - } - - - /*! - * \brief The relative permeability for the non-wetting phase - * of the medium implied by van Genuchten's - * parameterization. - * - * The permeability of gas in a 3p system equals the standard 2p description. - * (see p61. in "Comparison of the Three-Phase Oil Relative Permeability Models" - * MOJDEH DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.) - * - * \param Sw saturation of the water phase. - * \param Sn saturation of the NAPL phase. - * \param saturation saturation of the gas phase. - * \param params Array of parameters. - */ - static Scalar krg(const Params ¶ms, Scalar Sw, Scalar Sn, Scalar saturation) - { - - // Se = (Sw+Sn - Sgr)/(1-Sgr) - Scalar Se = std::min(((1-saturation) - params.Sgr()) / (1 - params.Sgr()), 1.); - - - /* regularization */ - if(Se > 1.0) return 0.0; - if(Se < 0.0) return 1.0; - Scalar scalFact = 1.; - if (saturation<=0.1) - { - scalFact = (saturation - params.Sgr())/(0.1 - params.Sgr()); - if (scalFact < 0.) scalFact = 0.; - } - - Scalar result = scalFact * std::pow(1 - Se, 1.0/3.) * std::pow(1 - std::pow(Se, 1/params.vgM()), 2*params.vgM()); - - return result; - } - - /*! - * \brief The relative permeability for a phase. - * \param Sw saturation of the water phase. - * \param Sg saturation of the gas phase. - * \param Sn saturation of the NAPL phase. - * \param params Array of parameters. - * \param phase indicator, The saturation of all phases. - */ - static Scalar kr(const Params ¶ms, const int phase, const Scalar Sw, const Scalar Sn, const Scalar Sg) - { - switch (phase) - { - case 0: - return krw(params, Sw, Sn, Sg); - break; - case 1: - return krn(params, Sw, Sn, Sg); - break; - case 2: - return krg(params, Sw, Sn, Sg); - break; - } - return 0; - } - - /* - * \brief the basis for calculating adsorbed NAPL in storage term - * \param bulk density of porous medium, adsorption coefficient - */ - static Scalar bulkDensTimesAdsorpCoeff (const Params ¶ms) - { - return params.rhoBulk() * params.KdNAPL(); - } -}; -} // namespace Opm - -#endif diff --git a/opm/material/fluidmatrixinteractions/3pAdapter.hpp b/opm/material/fluidmatrixinteractions/3pAdapter.hpp deleted file mode 100644 index 15bed6d6d..000000000 --- a/opm/material/fluidmatrixinteractions/3pAdapter.hpp +++ /dev/null @@ -1,105 +0,0 @@ -/* - Copyright (C) 2009-2013 by Andreas Lauser - Copyright (C) 2012 by Bernd Flemisch - - 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::ThreePAdapter - */ -#ifndef OPM_3P_ADAPTER_HPP -#define OPM_3P_ADAPTER_HPP - -#include -#include -#include - -#include - -namespace Opm { -/*! - * \ingroup FluidMatrixInteractions - * - * \brief Makes the three-phase capillary pressure-saturation relations - * available under the M-phase API for material laws. - */ -template -class ThreePAdapter -{ -public: - typedef typename ThreePLaw::Params Params; - typedef typename ThreePLaw::Scalar Scalar; - enum { numPhases = 3 }; - - /*! - * \brief The capillary pressure-saturation curve. - */ - template - static void capillaryPressures(ContainerT &values, - const Params ¶ms, - const FluidState &fluidState) - { - Scalar p_cgw = ThreePLaw::pCGW(params, fluidState.saturation(wPhaseIdx)); - Scalar p_cnw = ThreePLaw::pCNW(params, fluidState.saturation(wPhaseIdx)); - Scalar p_cgn = ThreePLaw::pCGN(params, - fluidState.saturation(wPhaseIdx) - + fluidState.saturation(nPhaseIdx)); - Scalar p_cAlpha = ThreePLaw::pCAlpha(params, - fluidState.saturation(nPhaseIdx)); - Scalar p_cnw1 = 0.0; - Valgrind::CheckDefined(p_cgw); - Valgrind::CheckDefined(p_cnw); - Valgrind::CheckDefined(p_cgn); - Valgrind::CheckDefined(p_cAlpha); - - values[gPhaseIdx] = 0; - values[nPhaseIdx] = - (p_cAlpha*p_cgn + (1 - p_cAlpha)*(p_cgw - p_cnw1)); - values[wPhaseIdx] = values[nPhaseIdx] - (p_cAlpha*p_cnw + (1 - p_cAlpha)*p_cnw1); - } - - static Scalar pCGW(const Params ¶ms, Scalar Sw) - { return ThreePLaw::pCGW(params, Sw); } - - /*! - * \brief The inverse of the capillary pressure-saturation curve. - */ - template - static void saturations(ContainerT &values, - const Params ¶ms, - const FluidState &fluidState) - { OPM_THROW(std::runtime_error, "Not implemented: Inverse capillary pressure curves"); } - - /*! - * \brief The relative permeability of all phases. - */ - template - static void relativePermeabilities(ContainerT &values, - const Params ¶ms, - const FluidState &fluidState) - { - Scalar Sw = fluidState.saturation(wPhaseIdx); - Scalar Sn = fluidState.saturation(nPhaseIdx); - Scalar Sg = fluidState.saturation(gPhaseIdx); - - values[wPhaseIdx] = ThreePLaw::krw(params, Sw, Sn, Sg); - values[nPhaseIdx] = ThreePLaw::krn(params, Sw, Sn, Sg); - values[gPhaseIdx] = ThreePLaw::krg(params, Sw, Sn, Sg); - } -}; -} // namespace Opm - -#endif diff --git a/opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp b/opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp new file mode 100644 index 000000000..ee9acf52e --- /dev/null +++ b/opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp @@ -0,0 +1,485 @@ +/* + Copyright (C) 2008-2013 by Andreas Lauser + Copyright (C) 2012 by Holger Class + Copyright (C) 2012 by Vishal Jambhekar + + 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::ThreePhaseParkerVanGenuchten + */ +#ifndef OPM_THREE_PHASE_PARKER_VAN_GENUCHTEN_HPP +#define OPM_THREE_PHASE_PARKER_VAN_GENUCHTEN_HPP + +#include "ThreePhaseParkerVanGenuchtenParams.hpp" + +#include + +namespace Opm { +/*! + * \ingroup FluidMatrixInteractions + * + * \brief Implementation of three-phase capillary pressure and + * relative permeability relations proposed by Parker and van + * Genuchten. + * + * Reference: J.B. Kool, J.C. Parker, M.Th. van Genuchten: Parameter + * Estimation for Unsaturated Flow and Transport Models -- A Review; + * Journal of Hydrology, 91 (1987) 255-293 + */ +template > +class ThreePhaseParkerVanGenuchten +{ +public: + static_assert(TraitsT::numPhases == 3, + "The number of phases considered by this capillary pressure " + "law is always three!"); + + typedef TraitsT Traits; + typedef ParamsT Params; + typedef typename Traits::Scalar Scalar; + + static const int numPhases = 3; + static const int wPhaseIdx = Traits::wPhaseIdx; + static const int nPhaseIdx = Traits::nPhaseIdx; + static const int gPhaseIdx = Traits::gPhaseIdx; + + //! 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 three phase capillary pressure law + * proposed by Parker and van Genuchten. + * + * This material law is valid for three fluid phases and only + * depends on the saturations. + * + * \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 &fluidState) + { + values[gPhaseIdx] = pcgn(params, fluidState); + values[nPhaseIdx] = 0; + values[wPhaseIdx] = - pcnw(params, fluidState); + } + + /*! + * \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 Scalar pcgn(const Params ¶ms, + const FluidState &fluidState) + { + Scalar PC_VG_REG = 0.01; + + // sum of liquid saturations + Scalar St = + fluidState.saturation(wPhaseIdx) + + fluidState.saturation(nPhaseIdx); + + Scalar Se = (St - params.Swrx())/(1. - params.Swrx()); + + // regularization + if (Se < 0.0) + Se=0.0; + if (Se > 1.0) + Se=1.0; + + if (Se>PC_VG_REG && Se<1-PC_VG_REG) + { + Scalar x = std::pow(Se,-1/params.vgM()) - 1; + return std::pow(x, 1 - params.vgM())/params.vgAlpha(); + } + + // value and derivative at regularization point + Scalar Se_regu; + if (Se<=PC_VG_REG) + Se_regu = PC_VG_REG; + else + Se_regu = 1-PC_VG_REG; + Scalar x = std::pow(Se_regu,-1/params.vgM())-1; + Scalar pc = std::pow(x, 1/params.vgN())/params.vgAlpha(); + Scalar pc_prime = + std::pow(x, 1/params.vgN()-1) + * std::pow(Se_regu,-1/params.vgM()-1) + / (-params.vgM()) + / params.vgAlpha() + / (1 - params.Sgr() - params.Swrx()) + / params.vgN(); + + // evaluate tangential + return ((Se-Se_regu)*pc_prime + pc)/params.betaGN(); + } + + /*! + * \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 Scalar pcnw(const Params ¶ms, + const FluidState &fluidState) + { + Scalar Sw = fluidState.saturation(wPhaseIdx); + Scalar Se = (Sw-params.Swr())/(1.-params.Snr()); + + Scalar PC_VG_REG = 0.01; + + // regularization + if (Se<0.0) + Se=0.0; + if (Se>1.0) + Se=1.0; + + if (Se>PC_VG_REG && Se<1-PC_VG_REG) { + Scalar x = std::pow(Se,-1/params.vgM()) - 1.0; + x = std::pow(x, 1 - params.vgM()); + return x/params.vgAlpha(); + } + + // value and derivative at regularization point + Scalar Se_regu; + if (Se<=PC_VG_REG) + Se_regu = PC_VG_REG; + else + Se_regu = 1.0 - PC_VG_REG; + + Scalar x = std::pow(Se_regu,-1/params.vgM())-1; + Scalar pc = std::pow(x, 1/params.vgN())/params.vgAlpha(); + Scalar pc_prime = + std::pow(x,1/params.vgN()-1) + * std::pow(Se_regu, -1.0/params.vgM() - 1) + / (-params.vgM()) + / params.vgAlpha() + / (1-params.Snr()-params.Swr()) + / params.vgN(); + + // evaluate tangential + return ((Se-Se_regu)*pc_prime + pc)/params.betaNW(); + } + + /*! + * \brief The saturation-capillary pressure curve. + * + */ + template + static void saturations(ContainerT &values, + const Params ¶ms, + const FluidState &fluidState) + { OPM_THROW(std::logic_error, "Not implemented: inverse capillary pressures"); } + + /*! + * \brief The saturation of the gas phase. + */ + template + static Scalar 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 Scalar 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 Scalar Sw(const Params ¶ms, + const FluidState &fluidState) + { OPM_THROW(std::logic_error, "Not implemented: Sw()"); } + + /*! + * \brief The relative permeability of all phases. + */ + template + static void relativePermeabilities(ContainerT &values, + const Params ¶ms, + const FluidState &fluidState) + { + values[wPhaseIdx] = krw(params, fluidState); + values[nPhaseIdx] = krn(params, fluidState); + values[gPhaseIdx] = krg(params, fluidState); + } + + /*! + * \brief The relative permeability for the wetting phase of the + * medium implied by van Genuchten's parameterization. + * + * The permeability of water in a 3p system equals the standard 2p description. + * (see p61. in "Comparison of the Three-Phase Oil Relative Permeability Models" + * MOJDEH DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.) + */ + template + static Scalar krw(const Params ¶ms, + const FluidState &fluidState) + { + + // transformation to effective saturation + Scalar Se = (fluidState.saturation(wPhaseIdx) - params.Swr()) / (1-params.Swr()); + + // regularization + if(Se > 1.0) return 1.; + if(Se < 0.0) return 0.; + + Scalar r = 1. - std::pow(1 - std::pow(Se, 1/params.vgM()), params.vgM()); + return std::sqrt(Se)*r*r; + } + + /*! + * \brief The relative permeability for the non-wetting phase + * due to the model of Parker et al. (1987). + * + * See model 7 of "Comparison of the Three-Phase Oil Relative + * Permeability Models" M. Delshad and G. A. Pope, Transport in + * Porous Media 4 (1989), 59-83; or -- more comprehensively -- + * "Estimation of primary drainage three-phase relative + * permeability for organic liquid transport in the vadose zone", + * L. I. Oliveira, A. H. Demond, Journal of Contaminant Hydrology + * 66 (2003), 261-285 + */ + template + static Scalar krn(const Params ¶ms, + const FluidState &fluidState) + { + Scalar Sn = fluidState.saturation(nPhaseIdx); + Scalar Sw = fluidState.saturation(wPhaseIdx); + Scalar Swe = std::min((Sw - params.Swr()) / (1 - params.Swr()), 1.); + Scalar Ste = std::min((Sw + Sn - params.Swr()) / (1 - params.Swr()), 1.); + + // regularization + if(Swe <= 0.0) Swe = 0.; + if(Ste <= 0.0) Ste = 0.; + if(Ste - Swe <= 0.0) return 0.; + + Scalar krn_; + krn_ = std::pow(1 - std::pow(Swe, 1/params.vgM()), params.vgM()); + krn_ -= std::pow(1 - std::pow(Ste, 1/params.vgM()), params.vgM()); + krn_ *= krn_; + + if (params.krRegardsSnr()) + { + // regard Snr in the permeability of the non-wetting + // phase, see Helmig1997 + Scalar resIncluded = + std::max(std::min(Sw - params.Snr() / (1-params.Swr()), 1.0), + 0.0); + krn_ *= std::sqrt(resIncluded ); + } + else + krn_ *= std::sqrt(Sn / (1 - params.Swr())); + + return krn_; + } + + + /*! + * \brief The relative permeability for the non-wetting phase + * of the medium implied by van Genuchten's + * parameterization. + * + * The permeability of gas in a three-phase system equals the + * standard two-phase description. (see p61. of "Comparison of the + * Three-Phase Oil Relative Permeability Models" M. Delshad and + * G. A. Pope, Transport in Porous Media 4 (1989), 59-83.) + */ + template + static Scalar krg(const Params ¶ms, + const FluidState &fluidState) + { + Scalar Sg = fluidState.saturation(gPhaseIdx); + Scalar Se = std::min(((1-Sg) - params.Sgr()) / (1 - params.Sgr()), 1.); + + // regularization + if(Se > 1.0) + return 0.0; + if(Se < 0.0) + return 1.0; + + Scalar scaleFactor = 1.; + if (Sg<=0.1) { + scaleFactor = (Sg - params.Sgr())/(0.1 - params.Sgr()); + if (scaleFactor < 0.) + scaleFactor = 0.; + } + + return scaleFactor + * std::pow(1 - Se, 1.0/3.) + * std::pow(1 - std::pow(Se, 1/params.vgM()), 2*params.vgM()); + } + + /*! + * \brief The derivative of all capillary pressures in regard to + * a given phase saturation. + */ + template + static void dCapillaryPressures_dSaturation(ContainerT &values, + const Params ¶ms, + const FluidState &state, + int satPhaseIdx) + { + OPM_THROW(std::logic_error, + "Not implemented: dCapillaryPressures_dSaturation()"); + } + + /*! + * \brief The derivative of all capillary pressures in regard to + * a given phase pressure. + */ + template + static void dCapillaryPressures_dPressure(ContainerT &values, + const Params ¶ms, + 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 + static void dCapillaryPressures_dTemperature(ContainerT &values, + const Params ¶ms, + 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 + static void dCapillaryPressures_dMoleFraction(ContainerT &values, + const Params ¶ms, + 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 + static void dRelativePermeabilities_dSaturation(ContainerT &values, + const Params ¶ms, + const FluidState &state, + int satPhaseIdx) + { + OPM_THROW(std::logic_error, + "Not implemented: dRelativePermeabilities_dSaturation()"); + } + + /*! + * \brief The derivative of all relative permeabilities in regard to + * a given phase pressure. + */ + template + static void dRelativePermeabilities_dPressure(ContainerT &values, + const Params ¶ms, + 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 + static void dRelativePermeabilities_dTemperature(ContainerT &values, + const Params ¶ms, + 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 + static void dRelativePermeabilities_dMoleFraction(ContainerT &values, + const Params ¶ms, + const FluidState &state, + int phaseIdx, + int compIdx) + { + // -> not composition dependent + for (int krPhaseIdx = 0; krPhaseIdx < numPhases; ++krPhaseIdx) + values[krPhaseIdx] = 0.0; + } +}; +} // namespace Opm + +#endif diff --git a/opm/material/fluidmatrixinteractions/3p/3pParkerVanGenuchtenParams.hpp b/opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchtenParams.hpp similarity index 63% rename from opm/material/fluidmatrixinteractions/3p/3pParkerVanGenuchtenParams.hpp rename to opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchtenParams.hpp index 4dcf0c14d..069b9757c 100644 --- a/opm/material/fluidmatrixinteractions/3p/3pParkerVanGenuchtenParams.hpp +++ b/opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchtenParams.hpp @@ -19,10 +19,10 @@ */ /*! * \file - * \copydoc Opm::ParkerVanGen3PParams + * \copydoc Opm::ThreePhaseParkerVanGenuchtenParams */ -#ifndef OPM_3P_PARKER_VAN_GENUCHTEN_PARAMS_HPP -#define OPM_3P_PARKER_VAN_GENUCHTEN_PARAMS_HPP +#ifndef OPM_THREE_PHASE_PARKER_VAN_GENUCHTEN_PARAMS_HPP +#define OPM_THREE_PHASE_PARKER_VAN_GENUCHTEN_PARAMS_HPP #include @@ -42,29 +42,31 @@ namespace Opm { * includes the residual saturations, as their handling is very * model-specific. */ -template -class ParkerVanGen3PParams +template +class ThreePhaseParkerVanGenuchtenParams { public: - typedef ScalarT Scalar; + typedef TraitsT Traits; + typedef typename Traits::Scalar Scalar; - ParkerVanGen3PParams() - {betaGW_ = betaNW_ = betaGN_ = 1.;} - - ParkerVanGen3PParams(Scalar vgAlpha, Scalar vgN, Scalar KdNAPL, Scalar rhoBulk, Dune::FieldVector residualSaturation, Scalar betaNW = 1., Scalar betaGN = 1., Scalar betaGW = 1., bool regardSnr=false) + ThreePhaseParkerVanGenuchtenParams() { - setVgAlpha(vgAlpha); - setVgN(vgN); - setSwr(residualSaturation[0]); - setSnr(residualSaturation[1]); - setSgr(residualSaturation[2]); - setSwrx(residualSaturation[3]); - setkrRegardsSnr(regardSnr); - setKdNAPL(KdNAPL); - setBetaNW(betaNW); - setBetaGN(betaGN); - setBetaGW(betaGW); - setRhoBulk(rhoBulk); + betaNW_ = 1.0; + betaGN_ = 1.0; + +#ifndef NDEBUG + finalized_ = false; +#endif + } + + /*! + * \brief Finish the initialization of the parameter object. + */ + void finalize() + { +#ifndef NDEBUG + finalized_ = true; +#endif } /*! @@ -72,7 +74,7 @@ public: * curve. */ Scalar vgAlpha() const - { return vgAlpha_; } + { assertFinalized_(); return vgAlpha_; } /*! * \brief Set the \f$\alpha\f$ shape parameter of van Genuchten's @@ -86,7 +88,7 @@ public: * curve. */ Scalar vgM() const - { return vgM_; } + { assertFinalized_(); return vgM_; } /*! * \brief Set the \f$m\f$ shape parameter of van Genuchten's @@ -102,7 +104,7 @@ public: * curve. */ Scalar vgN() const - { return vgN_; } + { assertFinalized_(); return vgN_; } /*! * \brief Set the \f$n\f$ shape parameter of van Genuchten's @@ -113,42 +115,11 @@ public: void setVgN(Scalar n) { vgN_ = n; vgM_ = 1 - 1/vgN_; } - /*! - * \brief Return the residual saturation. - */ - Scalar satResidual(int phaseIdx) const - { - switch (phaseIdx) - { - case 0: - return Swr_; - break; - case 1: - return Snr_; - break; - case 2: - return Sgr_; - break; - }; - OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx); - } - - /*! - * \brief Set all residual saturations. - */ - void setResiduals(Dune::FieldVector residualSaturation) - { - setSwr(residualSaturation[0]); - setSnr(residualSaturation[1]); - setSgr(residualSaturation[2]); - } - - /*! * \brief Return the residual wetting saturation. */ Scalar Swr() const - { return Swr_; } + { assertFinalized_(); return Swr_; } /*! * \brief Set the residual wetting saturation. @@ -160,7 +131,7 @@ public: * \brief Return the residual non-wetting saturation. */ Scalar Snr() const - { return Snr_; } + { assertFinalized_(); return Snr_; } /*! * \brief Set the residual non-wetting saturation. @@ -172,7 +143,7 @@ public: * \brief Return the residual gas saturation. */ Scalar Sgr() const - { return Sgr_; } + { assertFinalized_(); return Sgr_; } /*! * \brief Set the residual gas saturation. @@ -181,7 +152,7 @@ public: { Sgr_ = input; } Scalar Swrx() const - { return Swrx_; } + { assertFinalized_(); return Swrx_; } /*! * \brief Set the residual gas saturation. @@ -198,20 +169,14 @@ public: void setBetaGN(Scalar input) { betaGN_ = input; } - void setBetaGW(Scalar input) - { betaGW_ = input; } - /*! * \brief Return the values for the beta scaling parameters of capillary pressure between the phases */ Scalar betaNW() const - { return betaNW_; } + { assertFinalized_(); return betaNW_; } Scalar betaGN() const - { return betaGN_; } - - Scalar betaGW() const - { return betaGW_; } + { assertFinalized_(); return betaGN_; } /*! * \brief defines if residual n-phase saturation should be regarded in its relative permeability. @@ -222,32 +187,7 @@ public: * \brief Calls if residual n-phase saturation should be regarded in its relative permeability. */ bool krRegardsSnr() const - { return krRegardsSnr_; } - - - /*! - * \brief Return the bulk density of the porous medium - */ - Scalar rhoBulk() const - { return rhoBulk_; } - - /*! - * \brief Set the bulk density of the porous medium - */ - void setRhoBulk(Scalar input) - { rhoBulk_ = input; } - - /*! - * \brief Return the adsorption coefficient - */ - Scalar KdNAPL() const - { return KdNAPL_; } - - /*! - * \brief Set the adsorption coefficient - */ - void setKdNAPL(Scalar input) - { KdNAPL_ = input; } + { assertFinalized_(); return krRegardsSnr_; } void checkDefined() const { @@ -258,26 +198,32 @@ public: Valgrind::CheckDefined(Snr_); Valgrind::CheckDefined(Sgr_); Valgrind::CheckDefined(Swrx_); - Valgrind::CheckDefined(KdNAPL_); - Valgrind::CheckDefined(rhoBulk_); + Valgrind::CheckDefined(betaNW_); + Valgrind::CheckDefined(betaGN_); Valgrind::CheckDefined(krRegardsSnr_); } private: +#ifndef NDEBUG + void assertFinalized_() const + { assert(finalized_); } + + bool finalized_; +#else + void assertFinalized_() const + { } +#endif + Scalar vgAlpha_; Scalar vgM_; Scalar vgN_; Scalar Swr_; Scalar Snr_; Scalar Sgr_; - Scalar Swrx_; /* (Sw+Sn)_r */ - - Scalar KdNAPL_; - Scalar rhoBulk_; + Scalar Swrx_; // Swr + Snr Scalar betaNW_; Scalar betaGN_; - Scalar betaGW_; bool krRegardsSnr_ ; }; diff --git a/tests/test_fluidmatrixinteractions.cpp b/tests/test_fluidmatrixinteractions.cpp index 83cb8d105..6c46de2f9 100644 --- a/tests/test_fluidmatrixinteractions.cpp +++ b/tests/test_fluidmatrixinteractions.cpp @@ -36,6 +36,7 @@ #include #include #include +#include // include the helper classes to construct traits #include @@ -305,6 +306,12 @@ int main(int argc, char **argv) testThreePhaseApi(); //testThreePhaseSatApi(); } + { + typedef Opm::ThreePhaseParkerVanGenuchten MaterialLaw; + testGenericApi(); + testThreePhaseApi(); + //testThreePhaseSatApi(); + } { typedef Opm::NullMaterial MaterialLaw; testGenericApi();