Detailed Description

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

\[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*}

Definition at line 37 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

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. More...

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 Member Functions
double	computeValue (const double S_L) const

double	computeDerivative (const double S_L) const

double	computeSaturationForMinimumRelativePermeability () const

Private Attributes
const double	S_L_r_
	Residual saturation of wetting phase. More...

const double	S_L_max_
	Maximum saturation of wetting phase. More...

const double	m_
	Exponent \( m \). More...

const double	krel_min_
	Minimum relative permeability. More...

const double	S_L_for_krel_min_
	Liquid saturation that gives `krel_min_`. More...

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

PropertyDataType	value_
	The single value of a property. More...

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 24 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 56 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_.

◆ computeDerivative()

double MaterialPropertyLib::RelPermNonWettingPhaseVanGenuchtenMualem::computeDerivative ( const double S_L ) const

private

Computes

\[\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 limit zero.

Parameters

S_L	Liquid saturation.

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

Definition at line 84 of file RelPermNonWettingPhaseVanGenuchtenMualem.cpp.

 {
     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 m_, S_L_max_, and S_L_r_.

Referenced by computeSaturationForMinimumRelativePermeability(), and dValue().

◆ computeSaturationForMinimumRelativePermeability()

double MaterialPropertyLib::RelPermNonWettingPhaseVanGenuchtenMualem::computeSaturationForMinimumRelativePermeability ( ) const

private

Computes the saturation that gives the minimum relative permeability by using the Newton-Raphson method.

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

Definition at line 102 of file RelPermNonWettingPhaseVanGenuchtenMualem.cpp.

 {
     double S_for_k_rel_min = 0.5 * (S_L_r_ + S_L_max_);
     const double tolerance = 1.e-16;
     const double r0 = computeValue(S_for_k_rel_min) - krel_min_;
     for (int iterations = 0; iterations < 1000; ++iterations)
     {
         const double r = computeValue(S_for_k_rel_min) - krel_min_;
         S_for_k_rel_min = std::clamp(S_for_k_rel_min, S_L_r_, S_L_max_);
         S_for_k_rel_min -= r / computeDerivative(S_for_k_rel_min);
  
         if (std::fabs(r / r0) < tolerance)
         {
             return S_for_k_rel_min;
         }
     }
  
     OGS_FATAL(
         "The given minimum relative permeability, {:g}, is  not "
         "associated with the saturation  in the range of '('{:f}, {:f}')'. "
         "Please try another one in '(0, 1)', which should be close to zero",
         krel_min_, S_L_r_, S_L_max_);
 }

References computeDerivative(), computeValue(), krel_min_, OGS_FATAL, MathLib::r, S_L_max_, S_L_r_, and ParameterLib::tolerance.

◆ computeValue()

double MaterialPropertyLib::RelPermNonWettingPhaseVanGenuchtenMualem::computeValue ( const double S_L ) const

private

Parameters

S_L	Liquid saturation.

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

Definition at line 76 of file RelPermNonWettingPhaseVanGenuchtenMualem.cpp.

 {
     const double Se = (S_L - S_L_r_) / (S_L_max_ - S_L_r_);
     return std::sqrt(1.0 - Se) *
            std::pow(1.0 - std::pow(Se, 1.0 / m_), 2.0 * m_);
 }

References m_, S_L_max_, and S_L_r_.

Referenced by computeSaturationForMinimumRelativePermeability(), and value().

◆ 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 54 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 = std::get<double>(
         variable_array[static_cast<int>(Variable::liquid_saturation)]);
     if (S_L < S_L_r_ || S_L > S_L_for_krel_min_)
     {
         return 0.0;
     }
  
     return computeDerivative(S_L);
 }

References computeDerivative(), MaterialPropertyLib::liquid_saturation, OGS_FATAL, and S_L_for_krel_min_.

◆ 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 40 of file RelPermNonWettingPhaseVanGenuchtenMualem.cpp.

 {
     const double S_L = std::clamp(
         std::get<double>(
             variable_array[static_cast<int>(Variable::liquid_saturation)]),
         S_L_r_, S_L_max_);
  
     const double krel = computeValue(S_L);
     return std::max(krel_min_, krel);
 }