Implement a Mohr Coulomb plasticity model

doc/documentation/src/developer_guide/development_parts.rst

@@ -6,3 +6,4 @@ Development of specific code blocks
 .. toctree::
 
     materialdevelopment
+    mohr_coulomb_material

doc/documentation/src/developer_guide/mohr_coulomb_material.rst

@@ -0,0 +1,81 @@
+Mohr-Coulomb material
+=====================
+
+``MAT_Struct_MohrCoulomb`` is a finite-strain isotropic elastoplastic material. It uses the
+principal elastic logarithmic (Hencky) strains and their work-conjugate Kirchhoff stresses.
+The multiplicative history variable is the inverse plastic right Cauchy-Green tensor.
+
+Yield surface and flow rule
+---------------------------
+
+For ordered principal Kirchhoff stresses
+:math:`\tau_1 \geq \tau_2 \geq \tau_3`, the active Mohr-Coulomb plane is
+
+.. math::
+
+   f = k\tau_1-\tau_3-2\sqrt{k}\,c(\bar{\varepsilon}^{p}) \leq 0,
+   \qquad
+   k=\frac{1+\sin\phi}{1-\sin\phi}.
+
+The plastic potential has the same form with
+:math:`m=(1+\sin\psi)/(1-\sin\psi)`, where :math:`\psi` is the dilatancy angle.
+Setting ``DILATANCY_ANGLE`` equal to ``FRICTION_ANGLE`` gives associative plasticity.
+The return mapping works directly with active planes in principal stress space and treats a
+smooth face, both triaxial edges, and the hydrostatic apex separately.
+
+The cohesion hardening law combines linear and Voce terms:
+
+.. math::
+
+   c(\bar{\varepsilon}^{p}) =
+   c_0 + H\bar{\varepsilon}^{p}
+   + (c_\infty-c_0)\left(1-\exp(-b\bar{\varepsilon}^{p})\right).
+
+Set ``SATURATION_HARDENING`` to zero to disable the Voce term.
+
+Input
+-----
+
+.. code-block:: yaml
+
+   MAT_Struct_MohrCoulomb:
+     YOUNG_MODULUS: 1.0e5
+     POISSON_RATIO: 0.3
+     DENSITY: 2200
+     COHESION: 100
+     FRICTION_ANGLE: 0.5235987755982
+     DILATANCY_ANGLE: 0.4363323129986
+     LINEAR_HARDENING: 0
+     SATURATION_HARDENING: 500
+     HARDENING_EXP: 10
+     TOLERANCE: 1.0e-10
+     MAX_ITERATIONS: 100
+
+The material requires ``KINEM: nonlinear``. Angles are specified in radians. The dilatancy
+angle must be strictly positive and must not exceed the friction angle; a zero dilatancy angle
+does not provide an admissible hydrostatic apex flow direction.
+
+Output
+------
+
+The material registers ``plastic_strain``, ``accumulated_plastic_strain``,
+``accumulated_plastic_volumetric_strain``, and ``local_dissipated_energy`` as Gauss-point
+output quantities.
+
+Verification
+------------
+
+The calculated triaxial-compression and triaxial-tension meridians coincide with the analytical
+Mohr-Coulomb yield surface:
+
+.. image:: ../_assets/mohr_coulomb_yield_curve.png
+   :alt: Analytical and calculated Mohr-Coulomb yield meridians
+   :width: 80%
+
+The constitutive tests also cover elastic response, associative and non-associated flow,
+linear and Voce hardening, both edges, the apex, the consistent smooth-face tangent, output
+variables, and serialization.
+
+The return algorithm follows Clausen, Damkilde, and Andersen,
+*Efficient return algorithms for associated plasticity with multiple yield planes*,
+Computers & Structures 85 (2007), doi:10.1016/j.compstruc.2007.04.002.

src/core/legacy_enum_definitions/4C_legacy_enum_definitions_materials.cpp

@@ -101,6 +101,8 @@ std::string_view Core::Materials::to_string(Core::Materials::MaterialType materi
       return "MAT_Struct_ThermoStVenantK";
     case m_pldruckprag:
       return "MAT_Struct_DruckerPrager";
+    case m_plmohrcoulomb:
+      return "MAT_Struct_MohrCoulomb";
     case m_thermopllinelast:
       return "MAT_Struct_ThermoPlasticLinElast";
     case m_superelast:

src/core/legacy_enum_definitions/4C_legacy_enum_definitions_materials.hpp

@@ -133,6 +133,7 @@ namespace Core::Materials
     m_particle_wall_dem,      ///< particle wall material for DEM
     m_particle_pd,            ///< particle material for PD
     m_pldruckprag,       ///< Plastic linear elastic St.Venant Kirchhoff / Drucker Prager plasticity
+    m_plmohrcoulomb,     ///< Finite-strain logarithmic Mohr-Coulomb plasticity
     m_plgtn,             ///< Plastic linear elastic St.Venant Kirchhoff / GTN plasticity
     m_plelasthyper,      ///< general hyperelastic material for finite strain von-Mises plasticity
                          ///< using a semi-smooth Newton strategy (only in combination with such

src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp

@@ -20,12 +20,14 @@
 #include "4C_mat_inelastic_defgrad_factors_service.hpp"
 #include "4C_mat_micromaterial.hpp"
 #include "4C_mat_plasticdruckerprager.hpp"
+#include "4C_mat_plasticmohrcoulomb.hpp"
 #include "4C_mat_scatra_growth_remodel.hpp"
 #include "4C_mat_scatra_nonlocal_stimulus.hpp"
 #include "4C_porofluid_pressure_based_elast_scatra_input.hpp"
 #include "4C_structure_new_input.hpp"
 
 #include <filesystem>
+#include <numbers>
 #include <optional>
 #include <string>
 
@@ -58,6 +60,55 @@ std::unordered_map<Core::Materials::MaterialType, Core::IO::InputSpec> Global::v
         {.description = "Newtonian fluid"});
   }
 
+  /*----------------------------------------------------------------------*/
+  // Finite-strain logarithmic Mohr-Coulomb plasticity
+  {
+    using namespace Core::IO::InputSpecBuilders::Validators;
+    constexpr double half_pi = std::numbers::pi / 2.0;
+
+    known_materials[Core::Materials::m_plmohrcoulomb] = group("MAT_Struct_MohrCoulomb",
+        {
+            parameter<double>("YOUNG_MODULUS",
+                {.description = "Young's modulus", .validator = positive<double>()}),
+            parameter<double>("POISSON_RATIO",
+                {.description = "Poisson's ratio (must be in (-1, 0.5) for stability)",
+                    .validator = in_range(excl(-1.0), excl(0.5))}),
+            parameter<double>("DENSITY",
+                {.description = "Mass density", .validator = positive_or_zero<double>()}),
+            parameter<double>(
+                "COHESION", {.description = "Initial cohesion", .validator = positive<double>()}),
+            parameter<double>(
+                "FRICTION_ANGLE", {.description = "Friction angle in radians",
+                                      .validator = in_range(excl(0.0), excl(half_pi))}),
+            parameter<double>("DILATANCY_ANGLE",
+                {.description = "Dilatancy angle in radians "
+                                "(equal to friction angle for associated plasticity)",
+                    .validator = in_range(excl(0.0), excl(half_pi))}),
+            parameter<double>(
+                "LINEAR_HARDENING", {.description = "Linear isotropic cohesion-hardening modulus",
+                                        .default_value = 0.0,
+                                        .validator = positive_or_zero<double>()}),
+            parameter<double>("SATURATION_HARDENING",
+                {.description = "Asymptotic cohesion for Voce hardening; zero disables it",
+                    .default_value = 0.0,
+                    .validator = positive_or_zero<double>()}),
+            parameter<double>("HARDENING_EXP", {.description = "Voce hardening exponent",
+                                                   .default_value = 0.0,
+                                                   .validator = positive_or_zero<double>()}),
+            parameter<double>("TOLERANCE", {.description = "Local return-mapping tolerance",
+                                               .default_value = 1.0e-10,
+                                               .validator = positive<double>()}),
+            parameter<int>(
+                "MAX_ITERATIONS", {.description = "Maximum local return-mapping iterations",
+                                      .default_value = 50,
+                                      .validator = positive<int>()}),
+        },
+        {.description =
+                "Finite-strain Mohr-Coulomb plasticity based on elastic logarithmic strain and "
+                "Kirchhoff stress, with non-associated flow and linear plus Voce isotropic "
+                "cohesion hardening."});
+  }
+
   /*----------------------------------------------------------------------*/
   // Weakly compressible fluid according to Murnaghan-Tait
   {

src/mat/4C_mat_material_factory.cpp

@@ -114,6 +114,7 @@
 #include "4C_mat_plasticelasthyper.hpp"
 #include "4C_mat_plasticgtn.hpp"
 #include "4C_mat_plasticlinelast.hpp"
+#include "4C_mat_plasticmohrcoulomb.hpp"
 #include "4C_mat_plasticnlnlogneohooke.hpp"
 #include "4C_mat_poro_density_law.hpp"
 #include "4C_mat_poro_law.hpp"
@@ -249,6 +250,10 @@ std::unique_ptr<Core::Mat::PAR::Parameter> Mat::make_parameter(
     {
       return make_parameter_impl<Mat::PAR::PlasticDruckerPrager>(id, type, input_data);
     }
+    case Core::Materials::m_plmohrcoulomb:
+    {
+      return make_parameter_impl<Mat::PAR::PlasticMohrCoulomb>(id, type, input_data);
+    }
     case Core::Materials::m_plgtn:
     {
       return make_parameter_impl<Mat::PAR::PlasticGTN>(id, type, input_data);