[MPL/EDZ_permeability] Added a failure criterion for the tensile stress status

Documentation/ProjectFile/properties/property/PermeabilityMohrCoulombFailureIndexModel/t_tensile_strength_parameter.md

+\copydoc MaterialPropertyLib::PermeabilityMohrCoulombFailureIndexModel::t_sigma_max_

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

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)
     {

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_;
 };
 

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);
 }