ThreePhaseParkerVanGenuchten.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:
/*
  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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_THREE_PHASE_PARKER_VAN_GENUCHTEN_HPP
#define OPM_THREE_PHASE_PARKER_VAN_GENUCHTEN_HPP

#include "ThreePhaseParkerVanGenuchtenParams.hpp"

#include <opm/material/common/MathToolbox.hpp>

#include <algorithm>

namespace Opm {
template <class TraitsT,
          class ParamsT = ThreePhaseParkerVanGenuchtenParams<TraitsT> >
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 wettingPhaseIdx = Traits::wettingPhaseIdx;
    static const int nonWettingPhaseIdx = Traits::nonWettingPhaseIdx;
    static const int gasPhaseIdx = Traits::gasPhaseIdx;

    static const bool implementsTwoPhaseApi = false;

    static const bool implementsTwoPhaseSatApi = false;

    static const bool isSaturationDependent = true;

    static const bool isPressureDependent = false;

    static const bool isTemperatureDependent = false;

    static const bool isCompositionDependent = false;

    template <class ContainerT, class FluidState>
    static void capillaryPressures(ContainerT &values,
                                   const Params &params,
                                   const FluidState &fluidState)
    {
        typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;

        values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, fluidState);
        values[nonWettingPhaseIdx] = 0;
        values[wettingPhaseIdx] = - pcnw<FluidState, Evaluation>(params, fluidState);
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation pcgn(const Params &params, const FluidState &fluidState)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
        typedef MathToolbox<Evaluation> Toolbox;

        Scalar PC_VG_REG = 0.01;

        // sum of liquid saturations
        const auto& St =
            FsToolbox::template toLhs<Evaluation>(fluidState.saturation(wettingPhaseIdx))
            + FsToolbox::template toLhs<Evaluation>(fluidState.saturation(nonWettingPhaseIdx));

        Evaluation 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)
        {
            const Evaluation& x = Toolbox::pow(Se,-1/params.vgM()) - 1;
            return Toolbox::pow(x, 1.0 - 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;
        const Evaluation& x = std::pow(Se_regu,-1/params.vgM())-1;
        const Evaluation& pc = Toolbox::pow(x, 1.0/params.vgN())/params.vgAlpha();
        const Evaluation& pc_prime =
            Toolbox::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();
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation pcnw(const Params &params, const FluidState &fluidState)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
        typedef MathToolbox<Evaluation> Toolbox;

        const Evaluation& Sw =
            FsToolbox::template toLhs<Evaluation>(fluidState.saturation(wettingPhaseIdx));
        Evaluation 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) {
            Evaluation x = Toolbox::pow(Se,-1/params.vgM()) - 1.0;
            x = Toolbox::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;

        const Evaluation& x = std::pow(Se_regu,-1/params.vgM())-1;
        const Evaluation& pc = Toolbox::pow(x, 1/params.vgN())/params.vgAlpha();
        const Evaluation& pc_prime =
            Toolbox::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();
    }

    template <class ContainerT, class FluidState>
    static void saturations(ContainerT &/*values*/,
                            const Params &/*params*/,
                            const FluidState &/*fluidState*/)
    { OPM_THROW(std::logic_error, "Not implemented: inverse capillary pressures"); }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation Sg(const Params &/*params*/, const FluidState &/*fluidState*/)
    { OPM_THROW(std::logic_error, "Not implemented: Sg()"); }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation Sn(const Params &/*params*/, const FluidState &/*fluidState*/)
    { OPM_THROW(std::logic_error, "Not implemented: Sn()"); }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation Sw(const Params &/*params*/, const FluidState &/*fluidState*/)
    { OPM_THROW(std::logic_error, "Not implemented: Sw()"); }

    template <class ContainerT, class FluidState>
    static void relativePermeabilities(ContainerT &values,
                                       const Params &params,
                                       const FluidState &fluidState)
    {
        typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;

        values[wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
        values[nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
        values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation krw(const Params &params, const FluidState &fluidState)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
        typedef MathToolbox<Evaluation> Toolbox;

        const Evaluation& Sw =
            FsToolbox::template toLhs<Evaluation>(fluidState.saturation(wettingPhaseIdx));
        // transformation to effective saturation
        const Evaluation& Se = (Sw - params.Swr()) / (1-params.Swr());

        // regularization
        if(Se > 1.0) return 1.;
        if(Se < 0.0) return 0.;

        const Evaluation& r = 1. - Toolbox::pow(1 - Toolbox::pow(Se, 1/params.vgM()), params.vgM());
        return Toolbox::sqrt(Se)*r*r;
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation krn(const Params &params, const FluidState &fluidState)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
        typedef MathToolbox<Evaluation> Toolbox;

        const Evaluation& Sn =
            FsToolbox::template toLhs<Evaluation>(fluidState.saturation(nonWettingPhaseIdx));
        const Evaluation& Sw =
            FsToolbox::template toLhs<Evaluation>(fluidState.saturation(wettingPhaseIdx));
        Evaluation Swe = Toolbox::min((Sw - params.Swr()) / (1 - params.Swr()), 1.);
        Evaluation Ste = Toolbox::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.;

        Evaluation krn_;
        krn_ = Toolbox::pow(1 - Toolbox::pow(Swe, 1/params.vgM()), params.vgM());
        krn_ -= Toolbox::pow(1 - Toolbox::pow(Ste, 1/params.vgM()), params.vgM());
        krn_ *= krn_;

        if (params.krRegardsSnr())
        {
            // regard Snr in the permeability of the non-wetting
            // phase, see Helmig1997
            const Evaluation& resIncluded =
                Toolbox::max(Toolbox::min(Sw - params.Snr() / (1-params.Swr()), 1.0), 0.0);
            krn_ *= Toolbox::sqrt(resIncluded );
        }
        else
            krn_ *= Toolbox::sqrt(Sn / (1 - params.Swr()));

        return krn_;
    }


    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation krg(const Params &params, const FluidState &fluidState)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
        typedef MathToolbox<Evaluation> Toolbox;

        const Evaluation& Sg =
            FsToolbox::template toLhs<Evaluation>(fluidState.saturation(gasPhaseIdx));
        const Evaluation& Se = Toolbox::min(((1-Sg) - params.Sgr()) / (1 - params.Sgr()), 1.);

        // regularization
        if(Se > 1.0)
            return 0.0;
        if(Se < 0.0)
            return 1.0;

        Evaluation scaleFactor = 1.;
        if (Sg<=0.1) {
            scaleFactor = (Sg - params.Sgr())/(0.1 - params.Sgr());
            if (scaleFactor < 0.)
                return 0.0;
        }

        return scaleFactor
            * Toolbox::pow(1 - Se, 1.0/3.)
            * Toolbox::pow(1 - Toolbox::pow(Se, 1/params.vgM()), 2*params.vgM());
    }
};
} // namespace Opm

#endif