[GeoMechanicsApplication] Tests for Mohr Coulomb (#13302)

* Registered 2D and 3D laws for Mohr-Coulomb

* Fixes

* Added an extra region

* Added unit tests for the new zone

* clang format

* fixed documentation

* Added unit tests

* Changes based or review comments

* Changes based on review requests

* Fixed a bug + changes based on review

* fix in tests

* fixes based on review

* changes based on review

* fixes

* fixes

* fixes in readme

* used SetStrainSize/SetSpaceDimension and ransformPrincipalStressesToSigmaAndTau is taked out of three functions

* cleaned unit tests

* made the tolerance relative

* response to the review comments

* reverted properties settings

* Changes based on review

* typos

* Cleanup tests

* Added unit tests for 3D

* Fixes based on review

---------

Co-authored-by: markelov <gennady.markelov@deltares.nl>

Author: mnabideltares

applications/GeoMechanicsApplication/custom_constitutive/README.md

@@ -95,7 +95,8 @@ To incorporate the Mohr-Coulomb model with tensile cutoff in numerical simulatio
 
 4. Evaluate the condition and mapping
   - If the trial stress falls in the elastic zone, it stays unchanged. No mapping is applied.
-  - If the trial stress falls in the axial zone. The trial stress then needs to be mapped back to the cutoff curve.
+  - If the trial stress falls in the tensile apex return zone. The trial stress then needs to be mapped back to the apex.
+  - If the trial stress falls in the tensile cutoff zone. The trial stress then needs to be mapped back to the tension cutoff line.
   - If it falls in the tensile corner return zone, then it needs to be mapped to the corner point.
   - In the case of regular failure zone, then it is mapped back to the Mohr-Coulomb curve along the normal direction of flow function. The flow function is defined by
 

applications/GeoMechanicsApplication/custom_constitutive/documentation_data/Mohr-Coulomb-with-tension_cutoff-formulations.pdf

@@ -65,19 +65,27 @@ Vector CalculateCornerPoint(double FrictionAngle, double Cohesion, double Tensil
     return result;
 }
 
-Vector ReturnStressAtAxialZone(const Vector& rPrincipalTrialStressVector, double TensileStrength)
+Vector ReturnStressAtTensionApexReturnZone(const Vector& rPrincipalTrialStressVector, double TensileStrength)
 {
     auto result = rPrincipalTrialStressVector;
     result[0]   = TensileStrength;
-    result[2]   = TensileStrength - rPrincipalTrialStressVector[0] + rPrincipalTrialStressVector[2];
+    result[2]   = result[0];
     return result;
 }
 
+Vector ReturnStressAtTensionCutoffReturnZone(const Vector& rPrincipalTrialStressVector,
+                                             const Vector& rDerivativeOfFlowFunction,
+                                             double        TensileStrength)
+{
+    const auto lambda = (TensileStrength - rPrincipalTrialStressVector[0]) / rDerivativeOfFlowFunction[0];
+    return rPrincipalTrialStressVector + lambda * rDerivativeOfFlowFunction;
+}
+
 Vector ReturnStressAtCornerReturnZone(const Vector& rPrincipalTrialStressVector, const Vector& rCornerPoint)
 {
-    Vector result = rPrincipalTrialStressVector;
-    result[0]     = rCornerPoint[0] + rCornerPoint[1];
-    result[2]     = rCornerPoint[0] - rCornerPoint[1];
+    auto result = rPrincipalTrialStressVector;
+    result[0]   = rCornerPoint[0] + rCornerPoint[1];
+    result[2]   = rCornerPoint[0] - rCornerPoint[1];
     return result;
 }
 
@@ -122,7 +130,7 @@ Vector& MohrCoulombWithTensionCutOff::GetValue(const Variable<Vector>& rThisVari
     if (rThisVariable == CAUCHY_STRESS_VECTOR) {
         rValue = mStressVector;
     } else {
-        KRATOS_ERROR << "Can't get value of " << rThisVariable.Name() << ": unsupported variable\n";
+        rValue = ConstitutiveLaw::GetValue(rThisVariable, rValue);
     }
     return rValue;
 }
@@ -190,6 +198,27 @@ void MohrCoulombWithTensionCutOff::InitializeMaterial(const Properties& rMateria
     mTensionCutOff = TensionCutoff(rMaterialProperties[GEO_TENSILE_STRENGTH]);
 }
 
+void MohrCoulombWithTensionCutOff::InitializeMaterialResponseCauchy(Parameters& rValues)
+{
+    if (!mIsModelInitialized) {
+        mStressVectorFinalized = rValues.GetStressVector();
+        mStrainVectorFinalized = rValues.GetStrainVector();
+        mIsModelInitialized    = true;
+    }
+}
+
+void MohrCoulombWithTensionCutOff::GetLawFeatures(Features& rFeatures)
+{
+    rFeatures.SetOptions(mpConstitutiveDimension->GetSpatialType());
+    rFeatures.SetOptions(INFINITESIMAL_STRAINS);
+    rFeatures.SetOptions(ISOTROPIC);
+
+    rFeatures.SetStrainMeasure(StrainMeasure_Infinitesimal);
+
+    rFeatures.SetStrainSize(GetStrainSize());
+    rFeatures.SetSpaceDimension(WorkingSpaceDimension());
+}
+
 void MohrCoulombWithTensionCutOff::CalculateMaterialResponseCauchy(ConstitutiveLaw::Parameters& rParameters)
 {
     const auto& r_prop = rParameters.GetMaterialProperties();
@@ -209,12 +238,18 @@ void MohrCoulombWithTensionCutOff::CalculateMaterialResponseCauchy(ConstitutiveL
             CalculateCornerPoint(MathUtils<>::DegreesToRadians(r_prop[GEO_FRICTION_ANGLE]),
                                  r_prop[GEO_COHESION], r_prop[GEO_TENSILE_STRENGTH]);
 
-        if (IsStressAtAxialZone(principal_trial_stress_vector, r_prop[GEO_TENSILE_STRENGTH], apex, corner_point)) {
-            principal_trial_stress_vector =
-                ReturnStressAtAxialZone(principal_trial_stress_vector, r_prop[GEO_TENSILE_STRENGTH]);
-        } else if (IsStressAtCornerReturnZone(principal_trial_stress_vector,
-                                              MathUtils<>::DegreesToRadians(r_prop[GEO_DILATANCY_ANGLE]),
-                                              corner_point)) {
+        if (const auto trial_sigma_tau = TransformPrincipalStressesToSigmaAndTau(principal_trial_stress_vector);
+            IsStressAtTensionApexReturnZone(trial_sigma_tau, r_prop[GEO_TENSILE_STRENGTH], apex)) {
+            principal_trial_stress_vector = ReturnStressAtTensionApexReturnZone(
+                principal_trial_stress_vector, r_prop[GEO_TENSILE_STRENGTH]);
+        } else if (IsStressAtTensionCutoffReturnZone(trial_sigma_tau, r_prop[GEO_TENSILE_STRENGTH],
+                                                     apex, corner_point)) {
+            principal_trial_stress_vector = ReturnStressAtTensionCutoffReturnZone(
+                principal_trial_stress_vector,
+                mTensionCutOff.DerivativeOfFlowFunction(principal_trial_stress_vector),
+                r_prop[GEO_TENSILE_STRENGTH]);
+        } else if (IsStressAtCornerReturnZone(
+                       trial_sigma_tau, MathUtils<>::DegreesToRadians(r_prop[GEO_DILATANCY_ANGLE]), corner_point)) {
             principal_trial_stress_vector =
                 ReturnStressAtCornerReturnZone(principal_trial_stress_vector, corner_point);
         } else {
@@ -230,30 +265,41 @@ void MohrCoulombWithTensionCutOff::CalculateMaterialResponseCauchy(ConstitutiveL
 
     mStressVector = StressStrainUtilities::RotatePrincipalStresses(
         principal_trial_stress_vector, rotation_matrix, mpConstitutiveDimension->GetStrainSize());
+
+    rParameters.GetStressVector() = mStressVector;
 }
 
 bool MohrCoulombWithTensionCutOff::IsAdmissiblePrincipalStressState(const Vector& rPrincipalStresses) const
 {
-    return mCoulombYieldSurface.YieldFunctionValue(rPrincipalStresses) <= 0.0 &&
-           mTensionCutOff.YieldFunctionValue(rPrincipalStresses) <= 0.0;
+    constexpr auto tolerance          = 1.0e-10;
+    const auto coulomb_yield_function = mCoulombYieldSurface.YieldFunctionValue(rPrincipalStresses);
+    const auto tension_yield_function = mTensionCutOff.YieldFunctionValue(rPrincipalStresses);
+    const auto coulomb_tolerance      = tolerance * (1.0 + std::abs(coulomb_yield_function));
+    const auto tension_tolerance      = tolerance * (1.0 + std::abs(tension_yield_function));
+    return coulomb_yield_function < coulomb_tolerance && tension_yield_function < tension_tolerance;
+}
+
+bool MohrCoulombWithTensionCutOff::IsStressAtTensionApexReturnZone(const Vector& rTrialSigmaTau,
+                                                                   double        TensileStrength,
+                                                                   double        Apex) const
+{
+    return TensileStrength < Apex && rTrialSigmaTau[0] - rTrialSigmaTau[1] - TensileStrength > 0.0;
 }
 
-bool MohrCoulombWithTensionCutOff::IsStressAtAxialZone(const Vector& rPrincipalTrialStresses,
-                                                       double        TensileStrength,
-                                                       double        Apex,
-                                                       const Vector& rCornerPoint) const
+bool MohrCoulombWithTensionCutOff::IsStressAtTensionCutoffReturnZone(const Vector& rTrialSigmaTau,
+                                                                     double        TensileStrength,
+                                                                     double        Apex,
+                                                                     const Vector& rCornerPoint) const
 {
-    const auto trial_tau = TransformPrincipalStressesToSigmaAndTau(rPrincipalTrialStresses)[1];
-    return TensileStrength < Apex && trial_tau <= rCornerPoint[1] &&
-           mTensionCutOff.YieldFunctionValue(rPrincipalTrialStresses) >= 0.0;
+    return TensileStrength < Apex &&
+           rCornerPoint[1] - rTrialSigmaTau[1] - rCornerPoint[0] + rTrialSigmaTau[0] > 0.0;
 }
 
-bool MohrCoulombWithTensionCutOff::IsStressAtCornerReturnZone(const Vector& rPrincipalTrialStresses,
+bool MohrCoulombWithTensionCutOff::IsStressAtCornerReturnZone(const Vector& rTrialSigmaTau,
                                                               double        DilatancyAngle,
                                                               const Vector& rCornerPoint)
 {
-    const auto trial_sigma_tau = TransformPrincipalStressesToSigmaAndTau(rPrincipalTrialStresses);
-    return trial_sigma_tau[0] - rCornerPoint[0] - (trial_sigma_tau[1] - rCornerPoint[1]) * std::sin(DilatancyAngle) >= 0.0;
+    return rTrialSigmaTau[0] - rCornerPoint[0] - (rTrialSigmaTau[1] - rCornerPoint[1]) * std::sin(DilatancyAngle) >= 0.0;
 }
 
 Vector MohrCoulombWithTensionCutOff::CalculateTrialStressVector(const Vector& rStrainVector,

applications/GeoMechanicsApplication/custom_constitutive/documentation_data/mohr-coulomb-with-tension-cutoff-zones.svg

@@ -32,6 +32,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) MohrCoulombWithTensionCutOff : publi
 public:
     KRATOS_CLASS_POINTER_DEFINITION(MohrCoulombWithTensionCutOff);
 
+    MohrCoulombWithTensionCutOff() = default;
     explicit MohrCoulombWithTensionCutOff(std::unique_ptr<ConstitutiveLawDimension> pConstitutiveDimension);
 
     // Copying is not allowed. Use member `Clone` instead.
@@ -52,6 +53,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) MohrCoulombWithTensionCutOff : publi
     void                                   InitializeMaterial(const Properties&     rMaterialProperties,
                                                               const Geometry<Node>& rElementGeometry,
                                                               const Vector&         rShapeFunctionsValues) override;
+    void    InitializeMaterialResponseCauchy(Parameters& rValues) override;
+    void    GetLawFeatures(Features& rFeatures) override;
     Vector& GetValue(const Variable<Vector>& rThisVariable, Vector& rValue) override;
     using ConstitutiveLaw::GetValue;
     void SetValue(const Variable<Vector>& rVariable, const Vector& rValue, const ProcessInfo& rCurrentProcessInfo) override;
@@ -69,16 +72,20 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) MohrCoulombWithTensionCutOff : publi
     Vector                                    mStrainVectorFinalized;
     CoulombYieldSurface                       mCoulombYieldSurface;
     TensionCutoff                             mTensionCutOff;
+    bool                                      mIsModelInitialized = false;
 
     [[nodiscard]] Vector CalculateTrialStressVector(const Vector& rStrainVector,
                                                     double        YoungsModulus,
                                                     double        PoissonsRatio) const;
     [[nodiscard]] bool   IsAdmissiblePrincipalStressState(const Vector& rPrincipalStresses) const;
-    [[nodiscard]] bool   IsStressAtAxialZone(const Vector& rPrincipalTrialStresses,
-                                             double        TensileStrength,
-                                             double        Apex,
-                                             const Vector& rCornerPoint) const;
-    [[nodiscard]] static bool IsStressAtCornerReturnZone(const Vector& rPrincipalTrialStresses,
+    [[nodiscard]] bool   IsStressAtTensionApexReturnZone(const Vector& rTrialSigmaTau,
+                                                         double        TensileStrength,
+                                                         double        Apex) const;
+    [[nodiscard]] bool   IsStressAtTensionCutoffReturnZone(const Vector& rTrialSigmaTau,
+                                                           double        TensileStrength,
+                                                           double        Apex,
+                                                           const Vector& rCornerPoint) const;
+    [[nodiscard]] static bool IsStressAtCornerReturnZone(const Vector& rTrialSigmaTau,
                                                          double        DilatancyAngle,
                                                          const Vector& rCornerPoint);
 }; // Class MohrCoulombWithTensionCutOff

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp

@@ -366,6 +366,9 @@ void KratosGeoMechanicsApplication::Register()
 
     KRATOS_REGISTER_CONSTITUTIVE_LAW("Geo_IncrementalLinearElasticInterfaceLaw", mIncrementalLinearElasticInterfaceLaw)
 
+    KRATOS_REGISTER_CONSTITUTIVE_LAW("GeoMohrCoulombWithTensionCutOff2D", mMohrCoulombWithTensionCutOff2D)
+    KRATOS_REGISTER_CONSTITUTIVE_LAW("GeoMohrCoulombWithTensionCutOff3D", mMohrCoulombWithTensionCutOff3D)
+
     // Register Variables
     KRATOS_REGISTER_VARIABLE(VELOCITY_COEFFICIENT)
     KRATOS_REGISTER_VARIABLE(DT_PRESSURE_COEFFICIENT)

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.h

@@ -122,6 +122,7 @@
 #include "custom_constitutive/linear_elastic_2D_interface_law.h"
 #include "custom_constitutive/linear_elastic_3D_interface_law.h"
 #include "custom_constitutive/linear_elastic_plane_stress_2D_law.h"
+#include "custom_constitutive/mohr_coulomb_with_tension_cutoff.h"
 #include "custom_constitutive/plane_strain.h"
 #include "custom_constitutive/small_strain_udsm_2D_interface_law.hpp"
 #include "custom_constitutive/small_strain_udsm_2D_plane_strain_law.hpp"
@@ -970,6 +971,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) KratosGeoMechanicsApplication : publ
 
     const GeoIncrementalLinearElasticInterfaceLaw mIncrementalLinearElasticInterfaceLaw;
 
+    const MohrCoulombWithTensionCutOff mMohrCoulombWithTensionCutOff2D{std::make_unique<PlaneStrain>()};
+    const MohrCoulombWithTensionCutOff mMohrCoulombWithTensionCutOff3D{std::make_unique<ThreeDimensional>()};
     ///@}
 
 }; // Class KratosGeoMechanicsApplication