diff --git a/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/c_PermeabilityMohrCoulombFailureIndexModel.md b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/c_PermeabilityMohrCoulombFailureIndexModel.md
new file mode 100644
index 0000000000000000000000000000000000000000..b671162f260fc5a7223f73337a12ccc0541e10c3
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/c_PermeabilityMohrCoulombFailureIndexModel.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel
diff --git a/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_cohesion.md b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_cohesion.md
new file mode 100644
index 0000000000000000000000000000000000000000..7a1c464913b0fdb253c25e529b126627568e68ad
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_cohesion.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel::c_
diff --git a/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_fitting_factor.md b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_fitting_factor.md
new file mode 100644
index 0000000000000000000000000000000000000000..3a291b79aca744cb66bb5119af04ab64b4b37559
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_fitting_factor.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel::b_
diff --git a/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_friction_angle.md b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_friction_angle.md
new file mode 100644
index 0000000000000000000000000000000000000000..c74a04f77aa76905d2900d2c5e3b1c3d806c5295
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_friction_angle.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel::phi_
diff --git a/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_initial_permeability.md b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_initial_permeability.md
new file mode 100644
index 0000000000000000000000000000000000000000..b33969884fd89b8a0e0f481700b38a5dd51c68b1
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_initial_permeability.md
@@ -0,0 +1,2 @@
+\copydoc MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel::k0_
+Parameter name to be used as initial permeability.
diff --git a/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_maximum_permeability.md b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_maximum_permeability.md
new file mode 100644
index 0000000000000000000000000000000000000000..a0c095b15ffb81a737c17bc2da2fb764823437ac
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_maximum_permeability.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel::k_max_
diff --git a/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_reference_permeability.md b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_reference_permeability.md
new file mode 100644
index 0000000000000000000000000000000000000000..d3609db42d6bd2b5e22c3d43f5ab968ee303b13f
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_reference_permeability.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel::kr_
diff --git a/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_tensile_strength_parameter.md b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_tensile_strength_parameter.md
new file mode 100644
index 0000000000000000000000000000000000000000..ba0f3bdd3b56fa9fdd122a3f407daf84b214a3c4
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_tensile_strength_parameter.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel::t_sigma_max_
diff --git a/Documentation/bibliography/other.bib b/Documentation/bibliography/other.bib
index 3b8c3bd2de6af76225c2727e760cf153630c6602..f62df77bd36ee13f0866ce974acd5f6dde2abf8a 100644
--- a/Documentation/bibliography/other.bib
+++ b/Documentation/bibliography/other.bib
@@ -140,3 +140,28 @@
doi = "10.1007/s10596-013-9341-7",
publisher={Springer}
}
+
+@article{wangetalPerm2020,
+ title={Analysis of coupled thermal-hydro-mechanical processes during small
+ scale in situ heater experiment in Callovo-Oxfordian clay rock introducing
+ a failure-index permeability mode},
+ author={Wang, Wenqing and Shao, Hua and Nagel, Thomas and Fischer,
+ Thomas and Kolditz, Olaf},
+ journal={international journal of rock mechanics and mining sciences},
+ volume={under review},
+ number={},
+ pages={},
+ year={2020},
+}
+
+@article{labuz2012mohr,
+ title={Mohr--Coulomb failure criterion},
+ author={Labuz, Joseph F and Zang, Arno},
+ journal={Rock mechanics and rock engineering},
+ volume={45},
+ number={6},
+ pages={975--979},
+ year={2012},
+ publisher={Springer}
+}
+
diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp
index c91ee91a2c8d45c20dd46e4a2caade0cde551a34..bea2a9af96711165fd635adc814488f0acb04f0d 100644
--- a/MaterialLib/MPL/CreateProperty.cpp
+++ b/MaterialLib/MPL/CreateProperty.cpp
@@ -27,7 +27,7 @@
namespace
{
std::unique_ptr<MaterialPropertyLib::Property> createProperty(
- int const /*geometry_dimension*/,
+ int const geometry_dimension,
BaseLib::ConfigTree const& config,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
ParameterLib::CoordinateSystem const* const local_coordinate_system,
@@ -73,6 +73,13 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
return createIdealGasLaw(config);
}
+ if (boost::iequals(property_type,
+ "PermeabilityMohrCoulombFailureIndexModel"))
+ {
+ return createPermeabilityMohrCoulombFailureIndexModel(
+ geometry_dimension, config, parameters, local_coordinate_system);
+ }
+
if (boost::iequals(property_type, "PermeabilityOrthotropicPowerLaw"))
{
return createPermeabilityOrthotropicPowerLaw(config,
diff --git a/MaterialLib/MPL/Properties/CreatePermeabilityMohrCoulombFailureIndexModel.cpp b/MaterialLib/MPL/Properties/CreatePermeabilityMohrCoulombFailureIndexModel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..257376af9a8ec7562c550eb816b0519747fd9200
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreatePermeabilityMohrCoulombFailureIndexModel.cpp
@@ -0,0 +1,86 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on June 5, 2020, 9:07 AM
+ */
+
+#include "CreatePermeabilityMohrCoulombFailureIndexModel.h"
+
+#include <string>
+
+#include "BaseLib/ConfigTree.h"
+#include "Parameter.h"
+#include "ParameterLib/CoordinateSystem.h"
+#include "ParameterLib/Parameter.h"
+#include "ParameterLib/Utils.h"
+#include "PermeabilityMohrCoulombFailureIndexModel.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Property> createPermeabilityMohrCoulombFailureIndexModel(
+ int const geometry_dimension,
+ BaseLib::ConfigTree const& config,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+ ParameterLib::CoordinateSystem const* const local_coordinate_system)
+{
+ if ((geometry_dimension != 2) && (geometry_dimension != 3))
+ {
+ OGS_FATAL(
+ "The PermeabilityMohrCoulombFailureIndexModel is implemented only "
+ "for 2D or 3D problems");
+ }
+
+ //! \ogs_file_param{properties__property__type}
+ config.checkConfigParameter("type",
+ "PermeabilityMohrCoulombFailureIndexModel");
+
+ // Second access for storage.
+ //! \ogs_file_param{properties__property__name}
+ auto property_name = config.peekConfigParameter<std::string>("name");
+
+ DBUG("Create PermeabilityMohrCoulombFailureIndexModel property {:s}.",
+ property_name);
+
+ std::string const& parameter_name =
+ //! \ogs_file_param{properties__property__PermeabilityMohrCoulombFailureIndexModel__initial_permeability}
+ config.getConfigParameter<std::string>(
+ "initial_permeability");
+ auto const& parameter_k0 = ParameterLib::findParameter<double>(
+ parameter_name, parameters, 0, nullptr);
+
+ auto const kr =
+ //! \ogs_file_param{properties__property__PermeabilityMohrCoulombFailureIndexModel__reference_permeability}
+ config.getConfigParameter<double>("reference_permeability");
+ auto const b =
+ //! \ogs_file_param{properties__property__PermeabilityMohrCoulombFailureIndexModel__fitting_factor}
+ config.getConfigParameter<double>("fitting_factor");
+ auto const c =
+ //! \ogs_file_param{properties__property__PermeabilityMohrCoulombFailureIndexModel__cohesion}
+ config.getConfigParameter<double>("cohesion");
+ auto const phi =
+ //! \ogs_file_param{properties__property__PermeabilityMohrCoulombFailureIndexModel__friction_angle}
+ config.getConfigParameter<double>("friction_angle");
+ auto const max_k =
+ //! \ogs_file_param{properties__property__PermeabilityMohrCoulombFailureIndexModel__maximum_permeability}
+ config.getConfigParameter<double>("maximum_permeability");
+ auto const t_sigma_max =
+ //! \ogs_file_param{properties__property__PermeabilityMohrCoulombFailureIndexModel__tensile_strength_parameter}
+ config.getConfigParameter<double>("tensile_strength_parameter");
+
+ if (geometry_dimension == 2)
+ {
+ return std::make_unique<PermeabilityMohrCoulombFailureIndexModel<2>>(
+ std::move(property_name), parameter_k0, kr, b, c, phi, max_k,
+ t_sigma_max, local_coordinate_system);
+ }
+
+ return std::make_unique<PermeabilityMohrCoulombFailureIndexModel<3>>(
+ std::move(property_name), parameter_k0, kr, b, c, phi, max_k,
+ t_sigma_max, local_coordinate_system);
+}
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreatePermeabilityMohrCoulombFailureIndexModel.h b/MaterialLib/MPL/Properties/CreatePermeabilityMohrCoulombFailureIndexModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..9350b0a9ae907271d41c1b2d7145bf044e75d20c
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreatePermeabilityMohrCoulombFailureIndexModel.h
@@ -0,0 +1,40 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on June 5, 2020, 9:07 AM
+ */
+
+#pragma once
+
+#include <memory>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace ParameterLib
+{
+struct CoordinateSystem;
+struct ParameterBase;
+} // namespace ParameterLib
+
+namespace MaterialPropertyLib
+{
+class Property;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Property> createPermeabilityMohrCoulombFailureIndexModel(
+ int const geometry_dimension,
+ BaseLib::ConfigTree const& config,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+ ParameterLib::CoordinateSystem const* const local_coordinate_system);
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateProperties.h b/MaterialLib/MPL/Properties/CreateProperties.h
index d57b0c292ea7b838c2766f2da4f8c076b79561d9..bc85705aabaa88312e77b35791712c061bd2ed29 100644
--- a/MaterialLib/MPL/Properties/CreateProperties.h
+++ b/MaterialLib/MPL/Properties/CreateProperties.h
@@ -26,6 +26,7 @@
#include "CreateIdealGasLaw.h"
#include "CreateLinear.h"
#include "CreateParameter.h"
+#include "CreatePermeabilityMohrCoulombFailureIndexModel.h"
#include "CreatePermeabilityOrthotropicPowerLaw.h"
#include "CreatePorosityFromMassBalance.h"
#include "CreateSaturationDependentSwelling.h"
diff --git a/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..48ccccaf6b9421bfea13a808565d76f29e80feba
--- /dev/null
+++ b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp
@@ -0,0 +1,164 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on June 4, 2020, 10:13 AM
+ */
+
+#include "PermeabilityMohrCoulombFailureIndexModel.h"
+
+#include <boost/math/constants/constants.hpp>
+#include <cmath>
+#include <limits>
+
+#include "BaseLib/Error.h"
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
+#include "MaterialLib/MPL/Utils/GetSymmetricTensor.h"
+#include "MathLib/KelvinVector.h"
+#include "MathLib/MathTools.h"
+#include "ParameterLib/CoordinateSystem.h"
+#include "ParameterLib/Parameter.h"
+
+namespace MaterialPropertyLib
+{
+template <int DisplacementDim>
+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)
+ : 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);
+}
+
+template <int DisplacementDim>
+void PermeabilityMohrCoulombFailureIndexModel<DisplacementDim>::checkScale()
+ const
+{
+ if (!std::holds_alternative<Medium*>(scale_))
+ {
+ OGS_FATAL(
+ "The property 'PermeabilityMohrCoulombFailureIndexModel' is "
+ "implemented on the 'medium' scale only.");
+ }
+}
+
+template <int DisplacementDim>
+PropertyDataType
+PermeabilityMohrCoulombFailureIndexModel<DisplacementDim>::value(
+ VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& pos, double const t,
+ double const /*dt*/) const
+{
+ auto const& stress_vector = std::get<SymmetricTensor<DisplacementDim>>(
+ variable_array[static_cast<int>(Variable::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);
+ for (auto& k_i : k_data)
+ {
+ k_i = 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);
+}
+
+template <int DisplacementDim>
+PropertyDataType
+PermeabilityMohrCoulombFailureIndexModel<DisplacementDim>::dValue(
+ VariableArray const& /*variable_array*/, Variable const /*variable*/,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ 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.");
+}
+template class PermeabilityMohrCoulombFailureIndexModel<2>;
+template class PermeabilityMohrCoulombFailureIndexModel<3>;
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.h b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..0813a064ac28d9111554a14a3fa7039ec2aae44f
--- /dev/null
+++ b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.h
@@ -0,0 +1,131 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on June 4, 2020, 10:13 AM
+ */
+
+#pragma once
+
+#include "MaterialLib/MPL/Property.h"
+#include "MaterialLib/MPL/VariableType.h"
+
+namespace ParameterLib
+{
+struct CoordinateSystem;
+template <typename T>
+struct Parameter;
+} // namespace ParameterLib
+
+namespace MaterialPropertyLib
+{
+/**
+ * \brief A failure index dependent permeability model \cite wangetalPerm2020
+ *
+ * \f[ \mathbf{k} =
+ * \mathbf{k}_0+ H(f-1) k_\text{r} \mathrm{e}^{b f}\mathbf{I}\f]
+ *
+ * where
+ * \f$\mathbf{k}_0\f$ is the intrinsic permeability
+ * of the undamaged material,
+ * \f$H\f$ is the Heaviside step function, \f$f\f$ is the failure index,
+ * \f$k_\text{r}\f$ is a reference permeability,
+ * \f$b\f$ is a fitting parameter.
+ * \f$k_\text{r}\f$ and \f$b\f$ can be calibrated by experimental data.
+ *
+ * The failure index \f$f\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 \cite labuz2012mohr takes
+ * the form
+ * \f[\tau(\sigma)=c-\sigma \mathrm{tan} \phi\f]
+ * with \f$\tau\f$ the shear stress, \f$c\f$ the cohesion, \f$\sigma\f$ the
+ * normal stress, and \f$\phi\f$ the internal friction angle.
+ *
+ * The failure index of the Mohr Coulomb model is calculated by
+ * \f[
+ * f_{MC}=\frac{|\tau_m| }{\cos(\phi)\tau(\sigma_m)}
+ * \f]
+ * with
+ * \f$\tau_m=(\sigma_3-\sigma_1)/2\f$
+ * and \f$\sigma_m=(\sigma_1+\sigma_3)/2\f$,
+ * where \f$\sigma_1\f$ and \f$\sigma_3\f$ are the minimum and maximum shear
+ * stress, respectively.
+ *
+ * The tensile failure index is calculated by
+ * \f[
+ * f_{t} = \sigma_m / \sigma^t_{max}
+ * \f]
+ * with, \f$0 < \sigma^t_{max} < c \tan(\phi) \f$, 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 \f$\sigma_m>
+ * \sigma^t_{max}\f$. The failure index is then calculated by
+ * \f[
+ * f =
+ * \begin{cases}
+ * f_{MC}, & \sigma_{m} \leq \sigma^t_{max}\\
+ * max(f_{MC}, f_t), & \sigma_{m} > \sigma^t_{max}\\
+ * \end{cases}
+ * \f]
+ *
+ * The computed permeability components are restricted with an upper bound,
+ * i.e. \f$\mathbf{k}:=k_{ij} < k_{max}\f$.
+ *
+ * If \f$\mathbf{k}_0\f$ is orthogonal, i.e input two or three numbers
+ * for its diagonal entries, a coordinate system rotation of \f$\mathbf{k}\f$
+ * is possible if it is needed.
+ *
+ * Note: the conventional mechanics notations are used, which mean that tensile
+ * stress is positive.
+ *
+ */
+template <int DisplacementDim>
+class PermeabilityMohrCoulombFailureIndexModel final : public Property
+{
+public:
+ 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;
+
+private:
+ /// Intrinsic permeability for undamaged material. It can be a scalar or
+ /// tensor for anisotropic material.
+ ParameterLib::Parameter<double> const& k0_;
+ /// Reference permeability.
+ double const kr_;
+ /// Fitting parameter.
+ double const b_;
+ /// Cohesion.
+ double const c_;
+ /// Angle of internal friction.
+ double const phi_;
+
+ /// Maximum permeability.
+ double const k_max_;
+
+ /// Tensile strength parameter.
+ double const t_sigma_max_;
+ ParameterLib::CoordinateSystem const* const local_coordinate_system_;
+};
+
+} // namespace MaterialPropertyLib
diff --git a/Tests/MaterialLib/TestMPLPermeabilityMohrCoulombFailureIndexModel.cpp b/Tests/MaterialLib/TestMPLPermeabilityMohrCoulombFailureIndexModel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5cb62e30945fe140b37cfac11dcb333cc9179eb0
--- /dev/null
+++ b/Tests/MaterialLib/TestMPLPermeabilityMohrCoulombFailureIndexModel.cpp
@@ -0,0 +1,103 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on June 8, 2020, 9:17 AM
+ */
+
+#include <gtest/gtest.h>
+
+#include <Eigen/Eigen>
+#include <boost/math/constants/constants.hpp>
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Properties/CreatePermeabilityMohrCoulombFailureIndexModel.h"
+#include "MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
+#include "MathLib/KelvinVector.h"
+#include "ParameterLib/ConstantParameter.h"
+#include "TestMPL.h"
+#include "Tests/TestTools.h"
+
+TEST(MaterialPropertyLib, PermeabilityMohrCoulombFailureIndexModel)
+{
+ ParameterLib::ConstantParameter<double> const k0("k0", 1.e-20);
+ double const kr = 1.0e-19;
+ double const b = 3.0;
+ double const c = 1.0e+6;
+ double const phi = 15.0;
+ double const k_max = 1.e-10;
+ double const t_sigma_max = c;
+
+ auto const k_model = MPL::PermeabilityMohrCoulombFailureIndexModel<3>(
+ "k_f", k0, kr, b, c, phi, k_max, t_sigma_max, nullptr);
+
+ const int symmetric_tensor_size = 6;
+ using SymmetricTensor = Eigen::Matrix<double, symmetric_tensor_size, 1>;
+
+ // Under failure, i,e stress beyond the yield.
+ SymmetricTensor stress;
+ stress[0] = -1.36040e+7;
+ stress[1] = -1.78344e+7;
+ stress[2] = -1.36627e+7;
+ stress[3] = -105408;
+ stress[4] = -25377.2;
+ stress[5] = -1.39944e+7;
+
+ MPL::VariableArray vars;
+ vars[static_cast<int>(MPL::Variable::stress)].emplace<SymmetricTensor>(
+ stress);
+
+ ParameterLib::SpatialPosition const pos;
+ double const t = std::numeric_limits<double>::quiet_NaN();
+ double const dt = std::numeric_limits<double>::quiet_NaN();
+
+ auto const k = MPL::formEigenTensor<3>(k_model.value(vars, pos, t, dt));
+
+ double const k_expected = 1.1398264890628033e-15;
+
+ ASSERT_LE(std::fabs(k_expected - k(0, 0)) / k_expected, 1e-10)
+ << "for expected changed permeability " << k_expected
+ << " and for computed changed permeability." << k(0, 0);
+
+ // Stress outside of the yield surface. No change in permeability
+ stress[0] = -1.2e+7;
+ stress[1] = -1.5e+7;
+ stress[2] = -1.2e+7;
+ stress[3] = -1e+5;
+ stress[4] = -2200.2;
+ stress[5] = -8e+5;
+ vars[static_cast<int>(MPL::Variable::stress)] = stress;
+ auto const k_non_f =
+ MPL::formEigenTensor<3>(k_model.value(vars, pos, t, dt));
+
+ auto const k_non_f_expected = k0(t, pos)[0];
+
+ ASSERT_LE(std::fabs(k_non_f_expected - k_non_f(0, 0)) / k_non_f_expected,
+ 1e-19)
+ << "for expected untouched permeability" << k_non_f_expected
+ << " and for computed untouched permeability." << k_non_f(0, 0);
+
+ // Stress at the apex. No change in permeability.
+ const double val =
+ 2.0 * c / std::tan(phi * boost::math::constants::degree<double>());
+ stress[0] = 0.7 * val;
+ stress[1] = 0.3 * val;
+ stress[2] = 0.3 * val;
+ stress[3] = 0.0;
+ stress[4] = 0.0;
+ stress[5] = 0.0;
+ vars[static_cast<int>(MPL::Variable::stress)] = stress;
+ auto const k_apex_f =
+ MPL::formEigenTensor<3>(k_model.value(vars, pos, t, dt));
+
+ double const k_apex_expacted = 7.2849707418474819e-15;
+ ASSERT_LE(std::fabs(k_apex_expacted - k_apex_f(0, 0)) / k_apex_expacted,
+ 1e-19)
+ << "for expected untouched permeability" << k_non_f_expected
+ << " and for computed untouched permeability." << k_apex_f(0, 0);
+}