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

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp

@@ -15,131 +15,135 @@
 
 #include "custom_constitutive/coulomb_with_tension_cut_off_impl.h"
 #include "custom_utilities/constitutive_law_utilities.h"
+#include "custom_utilities/ublas_utilities.h"
 #include "geo_mechanics_application_variables.h"
 #include "includes/properties.h"
 #include "includes/serializer.h"
 
-namespace
+namespace Kratos
 {
 
-using namespace Kratos;
-
-double CalculateApex(double FrictionAngleInRadians, double Cohesion)
+CoulombWithTensionCutOffImpl::CoulombWithTensionCutOffImpl(const Properties& rMaterialProperties)
+    : mCoulombYieldSurface{rMaterialProperties}, mTensionCutOff{rMaterialProperties[GEO_TENSILE_STRENGTH]}
 {
-    return Cohesion / std::tan(FrictionAngleInRadians);
-}
+    if (rMaterialProperties.Has(GEO_ABS_YIELD_FUNCTION_TOLERANCE)) {
+        mAbsoluteYieldFunctionValueTolerance = rMaterialProperties[GEO_ABS_YIELD_FUNCTION_TOLERANCE];
+    }
 
-Vector CalculateCornerPoint(double FrictionAngleInRadians, double Cohesion, double TensileStrength)
-{
-    // Check whether the tension cut-off lies beyond the apex
-    auto result = Vector{ZeroVector(2)};
-    result[0]   = CalculateApex(FrictionAngleInRadians, Cohesion);
-    if (TensileStrength > result[0]) return result;
-
-    result[0] = (TensileStrength - Cohesion * std::cos(FrictionAngleInRadians)) /
-                (1.0 - std::sin(FrictionAngleInRadians));
-    result[1] = (Cohesion * std::cos(FrictionAngleInRadians) - TensileStrength * std::sin(FrictionAngleInRadians)) /
-                (1.0 - std::sin(FrictionAngleInRadians));
-    return result;
+    if (rMaterialProperties.Has(GEO_MAX_PLASTIC_ITERATIONS)) {
+        mMaxNumberOfPlasticIterations = rMaterialProperties[GEO_MAX_PLASTIC_ITERATIONS];
+    }
 }
 
-bool IsStressAtTensionApexReturnZone(const Vector& rTrialSigmaTau, double TensileStrength, double Apex)
+bool CoulombWithTensionCutOffImpl::IsAdmissibleSigmaTau(const Vector& rTrialSigmaTau) const
 {
-    return TensileStrength < Apex && rTrialSigmaTau[0] - rTrialSigmaTau[1] - TensileStrength > 0.0;
+    const auto coulomb_yield_function_value = mCoulombYieldSurface.YieldFunctionValue(rTrialSigmaTau);
+    const auto     tension_yield_function_value = mTensionCutOff.YieldFunctionValue(rTrialSigmaTau);
+    constexpr auto tolerance                    = 1.0e-10;
+    const auto     coulomb_tolerance = tolerance * (1.0 + std::abs(coulomb_yield_function_value));
+    const auto     tension_tolerance = tolerance * (1.0 + std::abs(tension_yield_function_value));
+    return coulomb_yield_function_value < coulomb_tolerance && tension_yield_function_value < tension_tolerance;
 }
 
-bool IsStressAtTensionCutoffReturnZone(const Vector& rTrialSigmaTau, double TensileStrength, double Apex, const Vector& rCornerPoint)
+Vector CoulombWithTensionCutOffImpl::DoReturnMapping(const Vector& rTrialSigmaTau,
+                                                     CoulombYieldSurface::CoulombAveragingType AveragingType,
+                                                     double& kappa_start)
 {
-    return TensileStrength < Apex &&
-           rCornerPoint[1] - rTrialSigmaTau[1] - rCornerPoint[0] + rTrialSigmaTau[0] > 0.0;
-}
+    Vector result = ZeroVector(2);
 
-bool IsStressAtCornerReturnZone(const Vector& rTrialSigmaTau,
-                                const Vector& rDerivativeOfFlowFunction,
-                                const Vector& rCornerPoint)
-{
-    return (rTrialSigmaTau[0] - rCornerPoint[0]) * rDerivativeOfFlowFunction[1] -
-               (rTrialSigmaTau[1] - rCornerPoint[1]) * rDerivativeOfFlowFunction[0] >=
-           0.0;
+    if (kappa_start < 0.0) {
+        kappa_start = mCoulombYieldSurface.GetKappa();
+    } else {
+        mCoulombYieldSurface.SetKappa(kappa_start);
+    }
+
+    for (auto counter = std::size_t{0}; counter < mMaxNumberOfPlasticIterations; ++counter) {
+        if (IsStressAtTensionApexReturnZone(rTrialSigmaTau)) {
+            return ReturnStressAtTensionApexReturnZone();
+        }
+
+        if (IsStressAtTensionCutoffReturnZone(rTrialSigmaTau)) {
+            return ReturnStressAtTensionCutoffReturnZone(rTrialSigmaTau);
+        }
+
+        if (IsStressAtCornerReturnZone(rTrialSigmaTau, AveragingType)) {
+            result = CalculateCornerPoint();
+        } else { // Regular failure region
+            result = ReturnStressAtRegularFailureZone(rTrialSigmaTau, AveragingType);
+        }
+
+        const auto kappa = kappa_start + mCoulombYieldSurface.CalculateEquivalentPlasticStrainIncrement(
+                                             rTrialSigmaTau, AveragingType);
+        mCoulombYieldSurface.SetKappa(kappa);
+
+        if (std::abs(mCoulombYieldSurface.YieldFunctionValue(result)) < mAbsoluteYieldFunctionValueTolerance) {
+            break;
+        }
+    }
+    return result;
 }
 
-Vector ReturnStressAtTensionApexReturnZone(double TensileStrength)
+Vector CoulombWithTensionCutOffImpl::CalculateCornerPoint() const
 {
-    auto result = Vector{ZeroVector{2}};
-    result[0]   = TensileStrength;
+    // Check whether the tension cut-off lies beyond the apex
+    auto result                 = Vector{ZeroVector(2)};
+    result[0]                   = mCoulombYieldSurface.CalculateApex();
+    const auto tensile_strength = mTensionCutOff.GetTensileStrength();
+    if (tensile_strength > result[0]) return result;
+
+    const auto cohesion = mCoulombYieldSurface.GetCohesion();
+    const auto sin_phi  = std::sin(mCoulombYieldSurface.GetFrictionAngleInRadians());
+    const auto cos_phi  = std::cos(mCoulombYieldSurface.GetFrictionAngleInRadians());
+    result[0]           = (tensile_strength - cohesion * cos_phi) / (1.0 - sin_phi);
+    result[1]           = (cohesion * cos_phi - tensile_strength * sin_phi) / (1.0 - sin_phi);
     return result;
 }
 
-Vector ReturnStressAtTensionCutoffReturnZone(const Vector& rSigmaTau,
-                                             const Vector& rDerivativeOfFlowFunction,
-                                             double        TensileStrength)
+bool CoulombWithTensionCutOffImpl::IsStressAtTensionApexReturnZone(const Vector& rTrialSigmaTau) const
 {
-    const auto lambda = (TensileStrength - rSigmaTau[0] - rSigmaTau[1]) /
-                        (rDerivativeOfFlowFunction[0] + rDerivativeOfFlowFunction[1]);
-    return rSigmaTau + lambda * rDerivativeOfFlowFunction;
+    const auto tensile_strength = mTensionCutOff.GetTensileStrength();
+    return tensile_strength < mCoulombYieldSurface.CalculateApex() &&
+           rTrialSigmaTau[0] - rTrialSigmaTau[1] - tensile_strength > 0.0;
 }
 
-Vector ReturnStressAtRegularFailureZone(const Vector& rSigmaTau,
-                                        const Vector& rDerivativeOfFlowFunction,
-                                        double        FrictionAngleInRadians,
-                                        double        Cohesion)
+bool CoulombWithTensionCutOffImpl::IsStressAtTensionCutoffReturnZone(const Vector& rTrialSigmaTau) const
 {
-    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;
-    return rSigmaTau + lambda * rDerivativeOfFlowFunction;
+    const auto corner_point = CalculateCornerPoint();
+    return mTensionCutOff.GetTensileStrength() < mCoulombYieldSurface.CalculateApex() &&
+           corner_point[1] - rTrialSigmaTau[1] - corner_point[0] + rTrialSigmaTau[0] > 0.0;
 }
 
-} // namespace
-
-namespace Kratos
+bool CoulombWithTensionCutOffImpl::IsStressAtCornerReturnZone(const Vector& rTrialSigmaTau,
+                                                              CoulombYieldSurface::CoulombAveragingType AveragingType) const
 {
+    const auto corner_point = CalculateCornerPoint();
+    const auto derivative_of_flow_function =
+        mCoulombYieldSurface.DerivativeOfFlowFunction(rTrialSigmaTau, AveragingType);
+    return (rTrialSigmaTau[0] - corner_point[0]) * derivative_of_flow_function[1] -
+               (rTrialSigmaTau[1] - corner_point[1]) * derivative_of_flow_function[0] >=
+           0.0;
+}
 
-CoulombWithTensionCutOffImpl::CoulombWithTensionCutOffImpl(const Properties& rMaterialProperties)
-    : mCoulombYieldSurface{rMaterialProperties}, mTensionCutOff{rMaterialProperties[GEO_TENSILE_STRENGTH]}
+Vector CoulombWithTensionCutOffImpl::ReturnStressAtTensionApexReturnZone() const
 {
+    return UblasUtilities::CreateVector({mTensionCutOff.GetTensileStrength(), 0.0});
 }
 
-bool CoulombWithTensionCutOffImpl::IsAdmissibleSigmaTau(const Vector& rTrialSigmaTau) const
+Vector CoulombWithTensionCutOffImpl::ReturnStressAtTensionCutoffReturnZone(const Vector& rSigmaTau) const
 {
-    const auto coulomb_yield_function_value = mCoulombYieldSurface.YieldFunctionValue(rTrialSigmaTau);
-    const auto     tension_yield_function_value = mTensionCutOff.YieldFunctionValue(rTrialSigmaTau);
-    constexpr auto tolerance                    = 1.0e-10;
-    const auto     coulomb_tolerance = tolerance * (1.0 + std::abs(coulomb_yield_function_value));
-    const auto     tension_tolerance = tolerance * (1.0 + std::abs(tension_yield_function_value));
-    return coulomb_yield_function_value < coulomb_tolerance && tension_yield_function_value < tension_tolerance;
+    const auto derivative_of_flow_function = mTensionCutOff.DerivativeOfFlowFunction(rSigmaTau);
+    const auto lambda_tc = (mTensionCutOff.GetTensileStrength() - rSigmaTau[0] - rSigmaTau[1]) /
+                           (derivative_of_flow_function[0] + derivative_of_flow_function[1]);
+    return rSigmaTau + lambda_tc * derivative_of_flow_function;
 }
 
-Vector CoulombWithTensionCutOffImpl::DoReturnMapping(const Vector& rTrialSigmaTau,
-                                                     CoulombYieldSurface::CoulombAveragingType AveragingType) const
+Vector CoulombWithTensionCutOffImpl::ReturnStressAtRegularFailureZone(const Vector& rSigmaTau,
+                                                                      CoulombYieldSurface::CoulombAveragingType AveragingType) const
 {
-    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)) {
-        return corner_point;
-    }
-
-    // Regular failure region
-    return ReturnStressAtRegularFailureZone(
-        rTrialSigmaTau, mCoulombYieldSurface.DerivativeOfFlowFunction(rTrialSigmaTau, AveragingType),
-        mCoulombYieldSurface.GetFrictionAngleInRadians(), mCoulombYieldSurface.GetCohesion());
+    const auto derivative_of_flow_function =
+        mCoulombYieldSurface.DerivativeOfFlowFunction(rSigmaTau, AveragingType);
+    const auto lambda = mCoulombYieldSurface.CalculatePlasticMultiplier(rSigmaTau, derivative_of_flow_function);
+    return rSigmaTau + lambda * derivative_of_flow_function;
 }
 
 void CoulombWithTensionCutOffImpl::save(Serializer& rSerializer) const

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.h

@@ -31,12 +31,25 @@ class CoulombWithTensionCutOffImpl
     explicit CoulombWithTensionCutOffImpl(const Properties& rMaterialProperties);
 
     [[nodiscard]] bool   IsAdmissibleSigmaTau(const Vector& rTrialSigmaTau) const;
-    [[nodiscard]] Vector DoReturnMapping(const Vector& rTrialSigmaTau,
-                                         CoulombYieldSurface::CoulombAveragingType AveragingType) const;
+    [[nodiscard]] Vector DoReturnMapping(const Vector&                             rTrialSigmaTau,
+                                         CoulombYieldSurface::CoulombAveragingType AveragingType,
+                                         double&                                   kappa_start);
 
 private:
     CoulombYieldSurface mCoulombYieldSurface;
     TensionCutoff       mTensionCutOff;
+    double              mAbsoluteYieldFunctionValueTolerance{1.0e-8};
+    std::size_t         mMaxNumberOfPlasticIterations{100};
+
+    [[nodiscard]] Vector CalculateCornerPoint() const;
+    [[nodiscard]] bool   IsStressAtTensionApexReturnZone(const Vector& rTrialSigmaTau) const;
+    [[nodiscard]] bool   IsStressAtTensionCutoffReturnZone(const Vector& rTrialSigmaTau) const;
+    [[nodiscard]] bool   IsStressAtCornerReturnZone(const Vector& rTrialSigmaTau,
+                                                    CoulombYieldSurface::CoulombAveragingType AveragingType) const;
+    [[nodiscard]] Vector ReturnStressAtTensionApexReturnZone() const;
+    [[nodiscard]] Vector ReturnStressAtTensionCutoffReturnZone(const Vector& rSigmaTau) const;
+    [[nodiscard]] Vector ReturnStressAtRegularFailureZone(const Vector& rSigmaTau,
+                                                          CoulombYieldSurface::CoulombAveragingType AveragingType) const;
 
     friend class Serializer;
     void save(Serializer& rSerializer) const;