[GeoMechanicsApplication] Implement hardening process for Mohr-Coulomb, based on the modified class structure

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp

@@ -24,6 +24,16 @@ namespace
 
 using namespace Kratos;
 
+double CalculatePlasticMultiplier(const Vector& rSigmaTau,
+                                        const Vector& rDerivativeOfFlowFunction,
+                                        double        FrictionAngleInRadians,
+                                        double        Cohesion)
+{
+    const auto sin_phi   = std::sin(FrictionAngleInRadians);
+    const auto numerator = sin_phi * rDerivativeOfFlowFunction[0] + rDerivativeOfFlowFunction[1];
+    return (Cohesion * std::cos(FrictionAngleInRadians) - rSigmaTau[0] * sin_phi - rSigmaTau[1]) / numerator;
+}
+
 double CalculateApex(double FrictionAngleInRadians, double Cohesion)
 {
     return Cohesion / std::tan(FrictionAngleInRadians);
@@ -84,10 +94,8 @@ Vector ReturnStressAtRegularFailureZone(const Vector& rSigmaTau,
                                         double        FrictionAngleInRadians,
                                         double        Cohesion)
 {
-    const auto sin_phi   = std::sin(FrictionAngleInRadians);
-    const auto numerator = sin_phi * rDerivativeOfFlowFunction[0] + rDerivativeOfFlowFunction[1];
     const auto lambda =
-        (Cohesion * std::cos(FrictionAngleInRadians) - rSigmaTau[0] * sin_phi - rSigmaTau[1]) / numerator;
+        CalculatePlasticMultiplier(rSigmaTau, rDerivativeOfFlowFunction, FrictionAngleInRadians, Cohesion);
     return rSigmaTau + lambda * rDerivativeOfFlowFunction;
 }
 
@@ -113,34 +121,55 @@ bool CoulombWithTensionCutOffImpl::IsAdmissibleSigmaTau(const Vector& rTrialSigm
 
 Vector CoulombWithTensionCutOffImpl::DoReturnMapping(const Properties& rProperties,
                                                      const Vector&     rTrialSigmaTau,
-                                                     CoulombYieldSurface::CoulombAveragingType AveragingType) const
+                                                     CoulombYieldSurface::CoulombAveragingType AveragingType)
 {
-    const auto apex = CalculateApex(mCoulombYieldSurface.GetFrictionAngleInRadians(),
-                                    mCoulombYieldSurface.GetCohesion());
-
-    if (IsStressAtTensionApexReturnZone(rTrialSigmaTau, mTensionCutOff.GetTensileStrength(), apex)) {
-        return ReturnStressAtTensionApexReturnZone(mTensionCutOff.GetTensileStrength());
-    }
-
-    const auto corner_point =
-        CalculateCornerPoint(mCoulombYieldSurface.GetFrictionAngleInRadians(),
-                             mCoulombYieldSurface.GetCohesion(), mTensionCutOff.GetTensileStrength());
-    if (IsStressAtTensionCutoffReturnZone(rTrialSigmaTau, mTensionCutOff.GetTensileStrength(), apex, corner_point)) {
-        return ReturnStressAtTensionCutoffReturnZone(
-            rTrialSigmaTau, mTensionCutOff.DerivativeOfFlowFunction(rTrialSigmaTau),
-            mTensionCutOff.GetTensileStrength());
-    }
-
-    if (IsStressAtCornerReturnZone(
+    Vector result = ZeroVector(2);
+
+    const double kappa_start = mCoulombYieldSurface.GetKappa();
+    double       error       = 1.0;
+    int          counter     = 0;
+    while (error > 1.0e-8) {
+        const auto apex = CalculateApex(mCoulombYieldSurface.GetFrictionAngleInRadians(),
+                                        mCoulombYieldSurface.GetCohesion());
+
+        if (IsStressAtTensionApexReturnZone(rTrialSigmaTau, mTensionCutOff.GetTensileStrength(), apex)) {
+            return ReturnStressAtTensionApexReturnZone(mTensionCutOff.GetTensileStrength());
+        }
+
+        const auto corner_point = CalculateCornerPoint(mCoulombYieldSurface.GetFrictionAngleInRadians(),
+                                                       mCoulombYieldSurface.GetCohesion(),
+                                                       mTensionCutOff.GetTensileStrength());
+        if (IsStressAtTensionCutoffReturnZone(rTrialSigmaTau, mTensionCutOff.GetTensileStrength(),
+                                              apex, corner_point)) {
+            return ReturnStressAtTensionCutoffReturnZone(
+                rTrialSigmaTau, mTensionCutOff.DerivativeOfFlowFunction(rTrialSigmaTau),
+                mTensionCutOff.GetTensileStrength());
+        }
+
+        if (IsStressAtCornerReturnZone(
+                rTrialSigmaTau, mCoulombYieldSurface.DerivativeOfFlowFunction(rTrialSigmaTau, AveragingType),
+                corner_point)) {
+            result = corner_point;
+        }
+
+        else { // Regular failure region
+            result = ReturnStressAtRegularFailureZone(
+                rTrialSigmaTau, mCoulombYieldSurface.DerivativeOfFlowFunction(rTrialSigmaTau, AveragingType),
+                mCoulombYieldSurface.GetFrictionAngleInRadians(), mCoulombYieldSurface.GetCohesion());
+        }
+
+        const auto lambda = CalculatePlasticMultiplier(
             rTrialSigmaTau, mCoulombYieldSurface.DerivativeOfFlowFunction(rTrialSigmaTau, AveragingType),
-            corner_point)) {
-        return corner_point;
+            mCoulombYieldSurface.GetFrictionAngleInRadians(), mCoulombYieldSurface.GetCohesion());
+        const double kappa =
+            kappa_start + CalculateEquivalentPlasticStrain(rTrialSigmaTau, AveragingType, lambda);
+        mCoulombYieldSurface.SetKappa(kappa);
+
+        error = std::abs(mCoulombYieldSurface.YieldFunctionValue(result));
+        counter++;
+        if (counter > 50) break;
     }
-
-    // Regular failure region
-    return ReturnStressAtRegularFailureZone(
-        rTrialSigmaTau, mCoulombYieldSurface.DerivativeOfFlowFunction(rTrialSigmaTau, AveragingType),
-        mCoulombYieldSurface.GetFrictionAngleInRadians(), mCoulombYieldSurface.GetCohesion());
+    return result;
 }
 
 void CoulombWithTensionCutOffImpl::save(Serializer& rSerializer) const
@@ -155,4 +184,18 @@ void CoulombWithTensionCutOffImpl::load(Serializer& rSerializer)
     rSerializer.load("TensionCutOff", mTensionCutOff);
 }
 
+double CoulombWithTensionCutOffImpl::CalculateEquivalentPlasticStrain(const Vector& rSigmaTau,
+                                                                      CoulombYieldSurface::CoulombAveragingType AveragingType,
+                                                                      double lambda) const
+{
+    Vector     dGdsigma   = mCoulombYieldSurface.DerivativeOfFlowFunction(rSigmaTau, AveragingType);
+    const auto g1         = (dGdsigma[0] + dGdsigma[1]) * 0.5;
+    const auto g3         = (dGdsigma[0] - dGdsigma[1]) * 0.5;
+    const auto mean       = (g1 + g3) / 3.0;
+    const auto deviatoric = std::sqrt(std::pow(g1 - mean, 2) + std::pow(g3 - mean, 2));
+    const auto alpha      = std::sqrt(2.0 / 3.0) * deviatoric;
+    return -alpha * lambda;
+    ;
+}
+
 } // namespace Kratos

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.h

@@ -31,14 +31,18 @@ class CoulombWithTensionCutOffImpl
     explicit CoulombWithTensionCutOffImpl(const Properties& rMaterialProperties);
 
     [[nodiscard]] bool   IsAdmissibleSigmaTau(const Vector& rTrialSigmaTau) const;
-    [[nodiscard]] Vector DoReturnMapping(const Properties& rProperties,
-                                         const Vector&     rTrialSigmaTau,
-                                         CoulombYieldSurface::CoulombAveragingType AveragingType) const;
+    [[nodiscard]] Vector DoReturnMapping(const Properties&                         rProperties,
+                                         const Vector&                             rTrialSigmaTau,
+                                         CoulombYieldSurface::CoulombAveragingType AveragingType);
 
 private:
     CoulombYieldSurface mCoulombYieldSurface;
     TensionCutoff       mTensionCutOff;
 
+    double CalculateEquivalentPlasticStrain(const Vector&                             rSigmaTau,
+                                            CoulombYieldSurface::CoulombAveragingType AveragingType,
+                                            double                                    lambda) const;
+
     friend class Serializer;
     void save(Serializer& rSerializer) const;
     void load(Serializer& rSerializer);

applications/GeoMechanicsApplication/custom_constitutive/coulomb_yield_surface.cpp

@@ -29,6 +29,11 @@ CoulombYieldSurface::KappaDependentFunction MakeConstantFunction(double Value)
     return [Value](double /* unused kappa */) { return Value; };
 }
 
+CoulombYieldSurface::KappaDependentFunction MakeLinearFunction(double Value, double Coefficient)
+{
+    return [Value, Coefficient](double kappa) { return std::max(Value + Coefficient * kappa, 0.0); };
+}
+
 std::string GetCoulombHardeningTypeFrom(const Properties& rMaterialProperties)
 {
     auto result   = rMaterialProperties[GEO_COULOMB_HARDENING_TYPE];
@@ -40,35 +45,41 @@ std::string GetCoulombHardeningTypeFrom(const Properties& rMaterialProperties)
 CoulombYieldSurface::KappaDependentFunction MakeFrictionAngleCalculator(const Properties& rMaterialProperties)
 {
     const auto hardening_type = GetCoulombHardeningTypeFrom(rMaterialProperties);
-
     if (hardening_type == "none") {
         return MakeConstantFunction(ConstitutiveLawUtilities::GetFrictionAngleInRadians(rMaterialProperties));
     }
-
+    if (hardening_type == "linear") {
+        return MakeLinearFunction(ConstitutiveLawUtilities::GetFrictionAngleInRadians(rMaterialProperties),
+                                  rMaterialProperties[GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS](0));
+    }
     KRATOS_ERROR << "Cannot create a kappa-dependent function for the friction angle of material "
                  << rMaterialProperties.Id() << ": unknown hardening type '" << hardening_type << "'\n";
 }
 
 CoulombYieldSurface::KappaDependentFunction MakeCohesionCalculator(const Properties& rMaterialProperties)
 {
     const auto hardening_type = GetCoulombHardeningTypeFrom(rMaterialProperties);
-
     if (hardening_type == "none") {
         return MakeConstantFunction(ConstitutiveLawUtilities::GetCohesion(rMaterialProperties));
     }
-
+    if (hardening_type == "linear") {
+        return MakeLinearFunction(ConstitutiveLawUtilities::GetCohesion(rMaterialProperties),
+                                  rMaterialProperties[GEO_COHESION_FUNCTION_COEFFICIENTS](0));
+    }
     KRATOS_ERROR << "Cannot create a kappa-dependent function for the cohesion of material "
                  << rMaterialProperties.Id() << ": unknown hardening type '" << hardening_type << "'\n";
 }
 
 CoulombYieldSurface::KappaDependentFunction MakeDilatancyAngleCalculator(const Properties& rMaterialProperties)
 {
     const auto hardening_type = GetCoulombHardeningTypeFrom(rMaterialProperties);
-
     if (hardening_type == "none") {
         return MakeConstantFunction(MathUtils<>::DegreesToRadians(rMaterialProperties[GEO_DILATANCY_ANGLE]));
     }
-
+    if (hardening_type == "linear") {
+        return MakeLinearFunction(MathUtils<>::DegreesToRadians(rMaterialProperties[GEO_DILATANCY_ANGLE]),
+                                  rMaterialProperties[GEO_DILATANCY_ANGLE_FUNCTION_COEFFICIENTS](0));
+    }
     KRATOS_ERROR << "Cannot create a kappa-dependent function for the dilatancy angle of material "
                  << rMaterialProperties.Id() << ": unknown hardening type '" << hardening_type << "'\n";
 }
@@ -111,6 +122,10 @@ double CoulombYieldSurface::GetDilatancyAngleInRadians() const
     return mDilatancyAngleCalculator(mKappa);
 }
 
+double CoulombYieldSurface::GetKappa() const { return mKappa; }
+
+void CoulombYieldSurface::SetKappa(double kappa) { mKappa = kappa; }
+
 double CoulombYieldSurface::YieldFunctionValue(const Vector& rSigmaTau) const
 {
     return rSigmaTau[1] + rSigmaTau[0] * std::sin(GetFrictionAngleInRadians()) -

applications/GeoMechanicsApplication/custom_constitutive/coulomb_yield_surface.h

@@ -41,6 +41,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) CoulombYieldSurface : public YieldSu
     [[nodiscard]] double GetFrictionAngleInRadians() const;
     [[nodiscard]] double GetCohesion() const;
     [[nodiscard]] double GetDilatancyAngleInRadians() const;
+    [[nodiscard]] double GetKappa() const;
+
+    void SetKappa(double kappa);
 
     [[nodiscard]] double YieldFunctionValue(const Vector& rSigmaTau) const override;
     [[nodiscard]] Vector DerivativeOfFlowFunction(const Vector&) const override;

applications/GeoMechanicsApplication/geo_mechanics_application.cpp

@@ -378,6 +378,9 @@ void KratosGeoMechanicsApplication::Register()
     KRATOS_REGISTER_VARIABLE(GEO_DILATANCY_ANGLE)
     KRATOS_REGISTER_VARIABLE(GEO_TENSILE_STRENGTH)
     KRATOS_REGISTER_VARIABLE(GEO_COULOMB_HARDENING_TYPE)
+    KRATOS_REGISTER_VARIABLE(GEO_COHESION_FUNCTION_COEFFICIENTS)
+    KRATOS_REGISTER_VARIABLE(GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS)
+    KRATOS_REGISTER_VARIABLE(GEO_DILATANCY_ANGLE_FUNCTION_COEFFICIENTS)
 
     KRATOS_REGISTER_VARIABLE(SPECIFIC_HEAT_CAPACITY_WATER)
     KRATOS_REGISTER_VARIABLE(SPECIFIC_HEAT_CAPACITY_SOLID)

applications/GeoMechanicsApplication/geo_mechanics_application_variables.cpp

@@ -36,6 +36,9 @@ KRATOS_CREATE_VARIABLE(double, GEO_FRICTION_ANGLE)
 KRATOS_CREATE_VARIABLE(double, GEO_DILATANCY_ANGLE)
 KRATOS_CREATE_VARIABLE(double, GEO_TENSILE_STRENGTH)
 KRATOS_CREATE_VARIABLE(std::string, GEO_COULOMB_HARDENING_TYPE)
+KRATOS_CREATE_VARIABLE(Vector, GEO_COHESION_FUNCTION_COEFFICIENTS)
+KRATOS_CREATE_VARIABLE(Vector, GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS)
+KRATOS_CREATE_VARIABLE(Vector, GEO_DILATANCY_ANGLE_FUNCTION_COEFFICIENTS)
 
 KRATOS_CREATE_VARIABLE(double, SPECIFIC_HEAT_CAPACITY_WATER)
 KRATOS_CREATE_VARIABLE(double, SPECIFIC_HEAT_CAPACITY_SOLID)

applications/GeoMechanicsApplication/geo_mechanics_application_variables.h

@@ -62,6 +62,9 @@ KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, GEO_FRICTI
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, GEO_DILATANCY_ANGLE)
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, GEO_TENSILE_STRENGTH)
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, std::string, GEO_COULOMB_HARDENING_TYPE)
+KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, Vector, GEO_COHESION_FUNCTION_COEFFICIENTS)
+KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, Vector, GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS)
+KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, Vector, GEO_DILATANCY_ANGLE_FUNCTION_COEFFICIENTS)
 
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, SPECIFIC_HEAT_CAPACITY_WATER)
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, SPECIFIC_HEAT_CAPACITY_SOLID)

applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_mohr_coulomb_with_tension_cutoff.cpp

@@ -15,6 +15,7 @@
 #include "custom_constitutive/plane_strain.h"
 #include "custom_constitutive/three_dimensional.h"
 #include "custom_utilities/registration_utilities.h"
+#include "custom_utilities/ublas_utilities.h"
 #include "geo_mechanics_application_variables.h"
 #include "tests/cpp_tests/geo_mechanics_fast_suite.h"
 #include "tests/cpp_tests/test_utilities.h"
@@ -622,4 +623,56 @@ KRATOS_TEST_CASE_IN_SUITE(MohrCoulombWithTensionCutOff_CalculateMaterialResponse
                               expected_cauchy_stress_vector, tolerance);
 }
 
+KRATOS_TEST_CASE_IN_SUITE(MohrCoulombWithTensionCutOff_CalculateMaterialResponseCauchyAtRegularFailureZoneWithHardening,
+                          KratosGeoMechanicsFastSuiteWithoutKernel)
+{
+    // Arrange
+    auto       law = MohrCoulombWithTensionCutOff(std::make_unique<PlaneStrain>());
+    Properties properties;
+    properties.SetValue(GEO_COULOMB_HARDENING_TYPE, "Linear");
+    properties.SetValue(GEO_FRICTION_ANGLE, 30.0);
+    properties.SetValue(GEO_COHESION, 10.0);
+    properties.SetValue(GEO_DILATANCY_ANGLE, 0.0);
+    properties.SetValue(GEO_TENSILE_STRENGTH, 10.0);
+
+    properties.SetValue(GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS, UblasUtilities::CreateVector({0.0}));
+    properties.SetValue(GEO_COHESION_FUNCTION_COEFFICIENTS, UblasUtilities::CreateVector({1.3301270189}));
+    properties.SetValue(GEO_DILATANCY_ANGLE_FUNCTION_COEFFICIENTS, UblasUtilities::CreateVector({0.0}));
+
+    ConstitutiveLaw::Parameters parameters;
+    parameters.SetMaterialProperties(properties);
+    const auto dummy_element_geometry      = Geometry<Node>{};
+    const auto dummy_shape_function_values = Vector{};
+    law.InitializeMaterial(properties, dummy_element_geometry, dummy_shape_function_values);
+
+    // Act and Assert
+    auto cauchy_stress_vector = UblasUtilities::CreateVector({10.0, 0.0, -40.0, 0.0});
+    auto expected_cauchy_stress_vector =
+        UblasUtilities::CreateVector({4.69104576482, 0.0, -34.69104576482, 0.0});
+    constexpr double tolerance = 1.0e-8;
+    KRATOS_EXPECT_VECTOR_NEAR(CalculateMappedStressVector(cauchy_stress_vector, parameters, law),
+                              expected_cauchy_stress_vector, tolerance);
+
+    // Arrange
+    properties.SetValue(GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS, UblasUtilities::CreateVector({1.0}));
+    properties.SetValue(GEO_COHESION_FUNCTION_COEFFICIENTS, UblasUtilities::CreateVector({0.0}));
+
+    // Act and Assert
+    cauchy_stress_vector = UblasUtilities::CreateVector({10.0, 0.0, -40.0, 0.0});
+    expected_cauchy_stress_vector =
+        UblasUtilities::CreateVector({6.811560518996, 0.0, -36.811560518996, 0.0});
+    KRATOS_EXPECT_VECTOR_NEAR(CalculateMappedStressVector(cauchy_stress_vector, parameters, law),
+                              expected_cauchy_stress_vector, tolerance);
+
+    // Arrange
+    properties.SetValue(GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS, UblasUtilities::CreateVector({0.0}));
+    properties.SetValue(GEO_DILATANCY_ANGLE_FUNCTION_COEFFICIENTS, UblasUtilities::CreateVector({1.0}));
+
+    // Act and Assert
+    cauchy_stress_vector = UblasUtilities::CreateVector({10.0, 0.0, -40.0, 0.0});
+    expected_cauchy_stress_vector = UblasUtilities::CreateVector({8.08509437348, 0.0, -38.08509437348, 0.0});
+    KRATOS_EXPECT_VECTOR_NEAR(CalculateMappedStressVector(cauchy_stress_vector, parameters, law),
+                              expected_cauchy_stress_vector, tolerance);
+}
+
 } // namespace Kratos::Testing