Detailed Description

template<int DisplacementDim>
class MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel< DisplacementDim >

A failure index dependent permeability model [41].

\[ \mathbf{k} = \mathbf{k}_0+ H(f-1) k_\text{r} \mathrm{e}^{b f}\mathbf{I}\]

where \(\mathbf{k}_0\) is the intrinsic permeability of the undamaged material, \(H\) is the Heaviside step function, \(f\) is the failure index, \(k_\text{r}\) is a reference permeability, \(b\) is a fitting parameter. \(k_\text{r}\) and \(b\) can be calibrated by experimental data.

The failure index \(f\) is calculated from the Mohr Coulomb failure criterion comparing an acting shear stress for the shear dominated failure. The tensile failure is governed by an input parameter of tensile_strength_parameter .

The Mohr Coulomb failure criterion [24] takes the form

\[\tau(\sigma)=c-\sigma \mathrm{tan} \phi\]

with \(\tau\) the shear stress, \(c\) the cohesion, \(\sigma\) the normal stress, and \(\phi\) the internal friction angle.

The failure index of the Mohr Coulomb model is calculated by

\[ f_{MC}=\frac{|\tau_m| }{\cos(\phi)\tau(\sigma_m)} \]

with \(\tau_m=(\sigma_3-\sigma_1)/2\) and \(\sigma_m=(\sigma_1+\sigma_3)/2\), where \(\sigma_1\) and \(\sigma_3\) are the minimum and maximum shear stress, respectively.

The tensile failure index is calculated by

\[ f_{t} = \sigma_m / \sigma^t_{max} \]

with, \(0 < \sigma^t_{max} < c \tan(\phi) \), a parameter of tensile strength for the cutting of the apex of the Mohr Coulomb model.

The tensile stress status is determined by a condition of \(\sigma_m> \sigma^t_{max}\). The failure index is then calculated by

\[ f = \begin{cases} f_{MC}, & \sigma_{m} \leq \sigma^t_{max}\\ max(f_{MC}, f_t), & \sigma_{m} > \sigma^t_{max}\\ \end{cases} \]

The computed permeability components are restricted with an upper bound, i.e. \(\mathbf{k}:=k_{ij} < k_{max}\).

If \(\mathbf{k}_0\) is orthogonal, i.e input two or three numbers for its diagonal entries, a coordinate system rotation of \(\mathbf{k}\) is possible if it is needed.

Note: the conventional mechanics notations are used, which mean that tensile stress is positive.

Definition at line 91 of file PermeabilityMohrCoulombFailureIndexModel.h.

#include <PermeabilityMohrCoulombFailureIndexModel.h>

Inheritance diagram for MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel< DisplacementDim >:

[legend]

Collaboration diagram for MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel< DisplacementDim >:

[legend]

Public Member Functions
	PermeabilityMohrCoulombFailureIndexModel (std::string name, ParameterLib::Parameter< double > const &k0, double const kr, double const b, double const c, double const phi, double const k_max, double const t_sigma_max, ParameterLib::CoordinateSystem const *const local_coordinate_system)

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.

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
ParameterLib::Parameter< double > const &	k0_

double const	kr_
	Reference permeability.

double const	b_
	Fitting parameter.

double const	c_
	Cohesion.

double const	phi_
	Angle of internal friction.

double const	k_max_
	Maximum permeability.

double const	t_sigma_max_
	Tensile strength parameter.

ParameterLib::CoordinateSystem const *const	local_coordinate_system_

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

◆ PermeabilityMohrCoulombFailureIndexModel()

template<int DisplacementDim>

MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel< DisplacementDim >::PermeabilityMohrCoulombFailureIndexModel	(	std::string	name,
		ParameterLib::Parameter< double > const &	k0,
		double const	kr,
		double const	b,
		double const	c,
		double const	phi,
		double const	k_max,
		double const	t_sigma_max,
		ParameterLib::CoordinateSystem const *const	local_coordinate_system )

Definition at line 31 of file PermeabilityMohrCoulombFailureIndexModel.cpp.

    : k0_(k0),
      kr_(kr),
      b_(b),
      c_(c),
      phi_(boost::math::constants::degree<double>() * phi),
      k_max_(k_max),
      t_sigma_max_(t_sigma_max),
      local_coordinate_system_(local_coordinate_system)
{
    const double t_sigma_upper = c_ / std::tan(phi_);
    if (t_sigma_max_ <= 0.0 || t_sigma_max_ > t_sigma_upper ||
        std::fabs(t_sigma_max_ - t_sigma_upper) <
            std::numeric_limits<double>::epsilon())
    {
        OGS_FATAL(
            "Tensile strength parameter of {:e} is out of the range (0, "
            "c/tan(phi)) = (0, {:e})",
            t_sigma_max_, t_sigma_upper);
    }
 
    name_ = std::move(name);
}

References MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel< DisplacementDim >::c_, MaterialPropertyLib::name, MaterialPropertyLib::Property::name_, OGS_FATAL, MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel< DisplacementDim >::phi_, and MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel< DisplacementDim >::t_sigma_max_.

Member Function Documentation

◆ checkScale()

template<int DisplacementDim>

void MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel< DisplacementDim >::checkScale ( ) const

overridevirtual

Reimplemented from MaterialPropertyLib::Property.

Definition at line 61 of file PermeabilityMohrCoulombFailureIndexModel.cpp.

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

References OGS_FATAL.

◆ dValue()

template<int DisplacementDim>

PropertyDataType MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel< DisplacementDim >::dValue	(	VariableArray const &	variable_array,
		Variable const	variable,
		ParameterLib::SpatialPosition const &	pos,
		double const	t,
		double const	dt ) const

overridevirtual

This virtual method will compute the property derivative value based on the variables that are passed as arguments with the default implementation using empty variables array for the previous time step.

The default implementation of this method only returns the property value derivative without altering it.

Reimplemented from MaterialPropertyLib::Property.

Definition at line 152 of file PermeabilityMohrCoulombFailureIndexModel.cpp.

{
    if (variable == Variable::mechanical_strain)
    {
        return 0.;
    }
 
    OGS_FATAL(
        "The derivative of the intrinsic permeability k(sigma, ...) with "
        "respect to stress tensor (sigma) is not implemented because that "
        "dk/du is normally omitted.");
}

References MaterialPropertyLib::mechanical_strain, and OGS_FATAL.

◆ value()

template<int DisplacementDim>

PropertyDataType MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel< DisplacementDim >::value	(	VariableArray const &	variable_array,
		ParameterLib::SpatialPosition const &	pos,
		double const	t,
		double const	dt ) const

overridevirtual

This virtual method will compute the property value based on the variables that are passed as arguments with the default implementation using empty variables array for the previous time step.

Reimplemented from MaterialPropertyLib::Property.

Definition at line 74 of file PermeabilityMohrCoulombFailureIndexModel.cpp.

{
    auto const& stress_vector =
        std::get<SymmetricTensor<DisplacementDim>>(variable_array.total_stress);
 
    auto const& stress_tensor =
        formEigenTensor<3>(static_cast<PropertyDataType>(stress_vector));
 
    Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>>
        eigenvalue_solver(stress_tensor);
 
    // Principle stress
    auto const sigma = eigenvalue_solver.eigenvalues();
 
    auto k_data = k0_(t, pos);
 
    double const max_sigma = std::max(std::fabs(sigma[0]), std::fabs(sigma[2]));
 
    if (max_sigma < std::numeric_limits<double>::epsilon())
    {
        return fromVector(k_data);
    }
 
    double const sigma_m = 0.5 * (sigma[2] + sigma[0]);
 
    double const tau_m = 0.5 * std::fabs(sigma[2] - sigma[0]);
    double f = 0.0;
    if (sigma_m > t_sigma_max_)
    {
        // tensile failure criterion
        f = sigma_m / t_sigma_max_;
 
        double const tau_tt =
            c_ * std::cos(phi_) - t_sigma_max_ * std::sin(phi_);
 
        f = std::max(f, tau_m / tau_tt);
    }
    else
    {
        // Mohr Coulomb failure criterion
        f = tau_m / (c_ * std::cos(phi_) - sigma_m * std::sin(phi_));
    }
 
    if (f >= 1.0)
    {
        const double exp_value = std::exp(b_ * f);
        std::transform(begin(k_data), end(k_data), begin(k_data),
                       [&](double const k_i)
                       { return std::min(k_i + kr_ * exp_value, k_max_); });
    }
 
    // Local coordinate transformation is only applied for the case that the
    // initial intrinsic permeability is given with orthotropic assumption.
    if (local_coordinate_system_ && (k_data.size() == DisplacementDim))
    {
        Eigen::Matrix<double, DisplacementDim, DisplacementDim> const e =
            local_coordinate_system_->transformation<DisplacementDim>(pos);
        Eigen::Matrix<double, DisplacementDim, DisplacementDim> k =
            Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
 
        for (int i = 0; i < DisplacementDim; ++i)
        {
            Eigen::Matrix<double, DisplacementDim, DisplacementDim> const
                ei_otimes_ei = e.col(i) * e.col(i).transpose();
 
            k += k_data[i] * ei_otimes_ei;
        }
        return k;
    }
 
    return fromVector(k_data);
}