Detailed Description

van Genuchten-Mualem relative permeability function for non-wetting phase in terms of effective wetting-phase saturation [25] :

\[k_{rel}^n= (1 - S_e)^{1/2} (1 - S_e^{1/m})^{2m}\]

with

\[S_e=\frac{S^L-S^L_r}{S^L_{\mbox{max}}-S^L_r}\]

where

\begin{eqnarray*} &S^L_r& \mbox{residual saturation of wetting phase,}\\ &S^L_{\mbox{max}}& \mbox{maximum saturation of wetting phase,}\\ &m\, \in (0, 1) & \mbox{ exponent.}\\ \end{eqnarray*}

The derivative of the relative permeability with respect to saturation is computed as

\[\frac{\mathrm{d} k_{rel}^n}{\mathrm{d}S^L}= -(\dfrac{[1-S_e^{1/m}]^{2m}}{2\sqrt{1-S_e}}+2\sqrt{1-S_e} {(1-S_e^{1/m})}^{2*m-1} S_e^{1/m-1})/(S^L_\mbox{max}-S^L_r) \]

As \(S^L \to S^L_\mbox{max}\), or \(S_e \to 1\), \(\dfrac{[1-S_e^{1/m}]^{2m}}{2\sqrt{1-S_e}}\) has a limit of zero.

Definition at line 47 of file RelPermNonWettingPhaseVanGenuchtenMualem.h.

#include <RelPermNonWettingPhaseVanGenuchtenMualem.h>

Inheritance diagram for MaterialPropertyLib::RelPermNonWettingPhaseVanGenuchtenMualem:

[legend]

Collaboration diagram for MaterialPropertyLib::RelPermNonWettingPhaseVanGenuchtenMualem:

[legend]

Public Member Functions
	RelPermNonWettingPhaseVanGenuchtenMualem (std::string name, const double S_L_r, const double S_n_r, const double m, const double krel_min)

void	checkScale () const override

PropertyDataType	value (VariableArray const &variable_array, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const override

PropertyDataType	dValue (VariableArray const &variable_array, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const override

double	computeSaturationForMinimumRelativePermeability () const

Public Member Functions inherited from MaterialPropertyLib::Property
virtual	~Property ()

virtual PropertyDataType	initialValue (ParameterLib::SpatialPosition const &pos, double const t) const

virtual PropertyDataType	value () const

virtual PropertyDataType	value (VariableArray const &variable_array, VariableArray const &variable_array_prev, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

virtual PropertyDataType	dValue (VariableArray const &variable_array, VariableArray const &variable_array_prev, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

virtual PropertyDataType	d2Value (VariableArray const &variable_array, Variable const variable1, Variable const variable2, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const
	Default implementation: 2nd derivative of any constant property is zero.

virtual void	setProperties (std::vector< std::unique_ptr< Phase > > const &phases)
	Default implementation:

void	setScale (std::variant< Medium , Phase , Component * > scale)

template<typename T >
T	initialValue (ParameterLib::SpatialPosition const &pos, double const t) const

template<typename T >
T	value () const

template<typename T >
T	value (VariableArray const &variable_array, VariableArray const &variable_array_prev, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

template<typename T >
T	value (VariableArray const &variable_array, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

template<typename T >
T	dValue (VariableArray const &variable_array, VariableArray const &variable_array_prev, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

template<typename T >
T	dValue (VariableArray const &variable_array, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

template<typename T >
T	d2Value (VariableArray const &variable_array, Variable const &variable1, Variable const &variable2, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

Private Attributes
const double	S_L_r_
	Residual saturation of wetting phase.

const double	S_L_max_
	Maximum saturation of wetting phase.

const double	m_
	Exponent \( m \).

const double	krel_min_
	Minimum relative permeability.

const double	S_L_for_krel_min_
	Liquid saturation that gives `krel_min_`.

Additional Inherited Members
Protected Attributes inherited from MaterialPropertyLib::Property
std::string	name_

PropertyDataType	value_
	The single value of a property.

PropertyDataType	dvalue_

std::variant< Medium , Phase , Component * >	scale_

Constructor & Destructor Documentation

◆ RelPermNonWettingPhaseVanGenuchtenMualem()

MaterialPropertyLib::RelPermNonWettingPhaseVanGenuchtenMualem::RelPermNonWettingPhaseVanGenuchtenMualem	(	std::string	name,
		const double	S_L_r,
		const double	S_n_r,
		const double	m,
		const double	krel_min )

Parameters

name	Name of the property,
S_L_r	Residual saturation of the wetting phase, \( S^L_r \)
S_n_r	Residual saturation of the non-wetting phase, \( S^n_{r} \)
m	Exponent, \( m \in [0,1]\)
krel_min	Minimum relative permeability, \( k_{rel}^n_{\mbox{min}}\)

Definition at line 32 of file RelPermNonWettingPhaseVanGenuchtenMualem.cpp.

    : S_L_r_(S_L_r),
      S_L_max_(1. - S_n_r),
      m_(m),
      krel_min_(krel_min),
      S_L_for_krel_min_(computeSaturationForMinimumRelativePermeability())
{
    name_ = std::move(name);
    checkVanGenuchtenExponentRange(m_);
}

References MaterialPropertyLib::checkVanGenuchtenExponentRange(), m_, MaterialPropertyLib::name, and MaterialPropertyLib::Property::name_.

Member Function Documentation

◆ checkScale()

void MaterialPropertyLib::RelPermNonWettingPhaseVanGenuchtenMualem::checkScale ( ) const

inlineoverridevirtual

Reimplemented from MaterialPropertyLib::Property.

Definition at line 66 of file RelPermNonWettingPhaseVanGenuchtenMualem.h.

    {
        if (!std::holds_alternative<Medium*>(scale_))
        {
            OGS_FATAL(
                "The property 'RelPermNonWettingPhaseVanGenuchtenMualem' is "
                "implemented on the 'media' scale only.");
        }
    }

References OGS_FATAL, and MaterialPropertyLib::Property::scale_.

◆ computeSaturationForMinimumRelativePermeability()

double MaterialPropertyLib::RelPermNonWettingPhaseVanGenuchtenMualem::computeSaturationForMinimumRelativePermeability ( ) const

Computes the saturation that gives the minimum relative permeability by using the Regula–Falsi Method.

Returns: \( S^L\) that gives the minimum relative permeability.

Definition at line 95 of file RelPermNonWettingPhaseVanGenuchtenMualem.cpp.

{
    auto f = [this](double S_L) -> double
    {
        return computeVanGenuchtenMualemValue(S_L, this->S_L_r_, this->S_L_max_,
                                              this->m_) -
               this->krel_min_;
    };
 
    auto root_finder =
        MathLib::Nonlinear::makeRegulaFalsi<MathLib::Nonlinear::Pegasus>(
            f, S_L_r_, S_L_max_);
 
    root_finder.step(1000);
 
    return root_finder.getResult();
}

References MaterialPropertyLib::computeVanGenuchtenMualemValue(), krel_min_, m_, S_L_max_, and S_L_r_.

◆ dValue()

PropertyDataType MaterialPropertyLib::RelPermNonWettingPhaseVanGenuchtenMualem::dValue	(	VariableArray const &	variable_array,
		Variable const	variable,
		ParameterLib::SpatialPosition const &	pos,
		double const	t,
		double const	dt ) const

overridevirtual

Returns: \( \frac{\mathrm{d} k_{rel}^n}{\mathrm{d}S^L} \).

Reimplemented from MaterialPropertyLib::Property.

Definition at line 61 of file RelPermNonWettingPhaseVanGenuchtenMualem.cpp.

{
    if (variable != Variable::liquid_saturation)
    {
        OGS_FATAL(
            "RelPermNonWettingPhaseVanGenuchtenMualem::dValue is implemented "
            "for the derivative with respect to liquid saturation only.");
    }
 
    const double S_L = variable_array.liquid_saturation;
    if (S_L < S_L_r_ || S_L > S_L_for_krel_min_)
    {
        return 0.0;
    }
 
    if (std::fabs(S_L - S_L_max_) < std::numeric_limits<double>::epsilon())
    {
        return 0.0;
    }
 
    const double Se = (S_L - S_L_r_) / (S_L_max_ - S_L_r_);
 
    const double val1 = std::sqrt(1.0 - Se);
    const double val2 = 1.0 - std::pow(Se, 1.0 / m_);
 
    return (-0.5 * std::pow(val2, 2.0 * m_) / val1 -
            2.0 * std::pow(Se, -1.0 + 1.0 / m_) * val1 *
                std::pow(val2, 2.0 * m_ - 1.0)) /
           (S_L_max_ - S_L_r_);
}

References MaterialPropertyLib::liquid_saturation, MaterialPropertyLib::VariableArray::liquid_saturation, m_, OGS_FATAL, S_L_for_krel_min_, S_L_max_, and S_L_r_.

◆ value()

PropertyDataType MaterialPropertyLib::RelPermNonWettingPhaseVanGenuchtenMualem::value	(	VariableArray const &	variable_array,
		ParameterLib::SpatialPosition const &	pos,
		double const	t,
		double const	dt ) const

overridevirtual

Returns: \(k_{rel}^n \).

Reimplemented from MaterialPropertyLib::Property.

Definition at line 48 of file RelPermNonWettingPhaseVanGenuchtenMualem.cpp.

{
    const double S_L =
        std::clamp(variable_array.liquid_saturation, S_L_r_, S_L_max_);
 
    const double krel =
        computeVanGenuchtenMualemValue(S_L, S_L_r_, S_L_max_, m_);
    return std::max(krel_min_, krel);
}