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/MaterialLib/MPL/Properties/CreatePermeabilityMohrCoulombFailureIndexModel.cpp b/MaterialLib/MPL/Properties/CreatePermeabilityMohrCoulombFailureIndexModel.cpp
index 1b2037ed5dee755274a99d001867966c42eb66ce..257376af9a8ec7562c550eb816b0519747fd9200 100644
--- a/MaterialLib/MPL/Properties/CreatePermeabilityMohrCoulombFailureIndexModel.cpp
+++ b/MaterialLib/MPL/Properties/CreatePermeabilityMohrCoulombFailureIndexModel.cpp
@@ -68,16 +68,19 @@ std::unique_ptr<Property> createPermeabilityMohrCoulombFailureIndexModel(
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,
- local_coordinate_system);
+ t_sigma_max, local_coordinate_system);
}
return std::make_unique<PermeabilityMohrCoulombFailureIndexModel<3>>(
std::move(property_name), parameter_k0, kr, b, c, phi, max_k,
- local_coordinate_system);
+ t_sigma_max, local_coordinate_system);
}
} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp
index f04b9217993b42aabbc70477c35c3e5e5be4760e..fc17f83128f7abd767a670060c7a1af9267fdb8f 100644
--- a/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp
+++ b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp
@@ -31,7 +31,7 @@ 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 k_max, double const t_sigma_max,
ParameterLib::CoordinateSystem const* const local_coordinate_system)
: k0_(k0),
kr_(kr),
@@ -39,8 +39,20 @@ PermeabilityMohrCoulombFailureIndexModel<DisplacementDim>::
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);
}
@@ -85,15 +97,24 @@ PermeabilityMohrCoulombFailureIndexModel<DisplacementDim>::value(
}
double const sigma_m = 0.5 * (sigma[2] + sigma[0]);
- double tau_tt = c_ * std::cos(phi_) - sigma_m * std::sin(phi_);
- const double apex_cut_offset = 0.001;
- if (std::fabs(tau_tt) < apex_cut_offset)
+
+ double const tau_m = 0.5 * std::fabs(sigma[2] - sigma[0]);
+ double f = 0.0;
+ if (sigma_m > t_sigma_max_)
{
- tau_tt = apex_cut_offset * c_ * std::cos(phi_);
- }
+ // tensile failure criterion
+ f = sigma_m / t_sigma_max_;
+
+ double const tau_tt =
+ c_ * std::cos(phi_) - t_sigma_max_ * std::sin(phi_);
- // +- tau_t = tau_tt
- double const f = 0.5 * std::fabs(sigma[2] - sigma[0]) / tau_tt;
+ f = std::max(f, tau_m / tau_tt);
+ }
+ else
+ {
+ // von Mohr Coulomb failure criterion
+ f = tau_m / (c_ * std::cos(phi_) - sigma_m * std::sin(phi_));
+ }
if (f >= 1.0)
{
diff --git a/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.h b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.h
index 2a03b63660416d5ee2b0cc003a0206515e919f89..80c0212a47c2ca8a7fadf3bef9502081cf86c971 100644
--- a/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.h
+++ b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.h
@@ -36,17 +36,22 @@ namespace MaterialPropertyLib
* \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. 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
* tensile stress, and \f$\phi\f$ the internal friction angle.
*
- * The failure index is calculated by
+ * The failure index of the Mohr Coulomb model is calculated by
* \f[
- * f=\frac{\tau_m }{\cos(\phi)\tau(\sigma_m)}
+ * f_{MC}=\frac{\tau_m }{\cos(\phi)\tau(\sigma_m)}
* \f]
* with
* \f$\tau_m=(\sigma_3-\sigma_1)/2\f$
@@ -54,6 +59,23 @@ namespace MaterialPropertyLib
* 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$\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=f_{MC}, & \sigma_{m} <\sigma^t_{max}\\
+ * f=max(f_{MC}, f_t), & \sigma_{m} \geq \sigma^t_{max}\\
+ * \end{cases}
+ * \f]
+ *
* Note: the conventional mechanics notations are used, which mean that tensile
* stress is positive.
*
@@ -65,7 +87,7 @@ 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 k_max, double const t_sigma_max,
ParameterLib::CoordinateSystem const* const local_coordinate_system);
void checkScale() const override;
@@ -92,6 +114,9 @@ private:
/// Maximum permeability.
double const k_max_;
+
+ /// Tensile strength parameter.
+ double const t_sigma_max_;
ParameterLib::CoordinateSystem const* const local_coordinate_system_;
};
diff --git a/Tests/MaterialLib/TestMPLPermeabilityMohrCoulombFailureIndexModel.cpp b/Tests/MaterialLib/TestMPLPermeabilityMohrCoulombFailureIndexModel.cpp
index 12aa77dea390656427284b4de45d877d1711b385..5cb62e30945fe140b37cfac11dcb333cc9179eb0 100644
--- a/Tests/MaterialLib/TestMPLPermeabilityMohrCoulombFailureIndexModel.cpp
+++ b/Tests/MaterialLib/TestMPLPermeabilityMohrCoulombFailureIndexModel.cpp
@@ -9,8 +9,6 @@
* Created on June 8, 2020, 9:17 AM
*/
-#pragma once
-
#include <gtest/gtest.h>
#include <Eigen/Eigen>
@@ -27,16 +25,16 @@
TEST(MaterialPropertyLib, PermeabilityMohrCoulombFailureIndexModel)
{
- ParameterLib::ConstantParameter<double> k0 =
- ParameterLib::ConstantParameter<double>("k0", 1.e-20);
+ 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, nullptr);
+ "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>;
@@ -62,7 +60,7 @@ TEST(MaterialPropertyLib, PermeabilityMohrCoulombFailureIndexModel)
double const k_expected = 1.1398264890628033e-15;
- ASSERT_LE(std::fabs((k_expected - k(0, 0))) / k_expected, 1e-10)
+ 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);
@@ -97,7 +95,9 @@ TEST(MaterialPropertyLib, PermeabilityMohrCoulombFailureIndexModel)
auto const k_apex_f =
MPL::formEigenTensor<3>(k_model.value(vars, pos, t, dt));
- ASSERT_LE(std::fabs(k_max - k_apex_f(0, 0)) / k_max, 1e-19)
+ 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);
}