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

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp

@@ -14,83 +14,25 @@
 //
 
 #include "custom_constitutive/coulomb_with_tension_cut_off_impl.h"
-#include "custom_processes/set_automated_initial_variable_process.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
-{
-
-using namespace Kratos;
-
-Vector CalculateCornerPoint(double FrictionAngleInRadians, double Cohesion, double TensileStrength, double Apex)
-{
-    // Check whether the tension cut-off lies beyond the apex
-    auto result = Vector{ZeroVector(2)};
-    result[0]   = Apex;
-    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;
-}
-
-bool IsStressAtTensionApexReturnZone(const Vector& rTrialSigmaTau, double TensileStrength, double Apex)
-{
-    return TensileStrength < Apex && rTrialSigmaTau[0] - rTrialSigmaTau[1] - TensileStrength > 0.0;
-}
-
-bool IsStressAtTensionCutoffReturnZone(const Vector& rTrialSigmaTau, double TensileStrength, double Apex, const Vector& rCornerPoint)
-{
-    return TensileStrength < Apex &&
-           rCornerPoint[1] - rTrialSigmaTau[1] - rCornerPoint[0] + rTrialSigmaTau[0] > 0.0;
-}
-
-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;
-}
-
-Vector ReturnStressAtTensionApexReturnZone(double TensileStrength)
-{
-    auto result = Vector{ZeroVector{2}};
-    result[0]   = TensileStrength;
-    return result;
-}
-
-Vector ReturnStressAtTensionCutoffReturnZone(const Vector& rSigmaTau,
-                                             const Vector& rDerivativeOfFlowFunction,
-                                             double        TensileStrength)
-{
-    const auto lambda = (TensileStrength - rSigmaTau[0] - rSigmaTau[1]) /
-                        (rDerivativeOfFlowFunction[0] + rDerivativeOfFlowFunction[1]);
-    return rSigmaTau + lambda * rDerivativeOfFlowFunction;
-}
-
-Vector ReturnStressAtRegularFailureZone(const Vector&        rSigmaTau,
-                                        CoulombYieldSurface& rCoulombYieldSurface,
-                                        const Vector&        rDerivativeOfFlowFunction)
-{
-    const auto lambda = rCoulombYieldSurface.CalculatePlasticMultiplier(rSigmaTau, rDerivativeOfFlowFunction);
-    return rSigmaTau + lambda * rDerivativeOfFlowFunction;
-}
-
-} // namespace
-
 namespace Kratos
 {
 
 CoulombWithTensionCutOffImpl::CoulombWithTensionCutOffImpl(const Properties& rMaterialProperties)
     : mCoulombYieldSurface{rMaterialProperties}, mTensionCutOff{rMaterialProperties[GEO_TENSILE_STRENGTH]}
 {
+    if (rMaterialProperties.Has(GEO_ABS_YIELD_FUNCTION_TOLERANCE)) {
+        mAbsoluteYieldFunctionValueTolerance = rMaterialProperties[GEO_ABS_YIELD_FUNCTION_TOLERANCE];
+    }
+
+    if (rMaterialProperties.Has(GEO_MAX_PLASTIC_ITERATIONS)) {
+        mMaxNumberOfPlasticIterations = rMaterialProperties[GEO_MAX_PLASTIC_ITERATIONS];
+    }
 }
 
 bool CoulombWithTensionCutOffImpl::IsAdmissibleSigmaTau(const Vector& rTrialSigmaTau) const
@@ -104,51 +46,106 @@ bool CoulombWithTensionCutOffImpl::IsAdmissibleSigmaTau(const Vector& rTrialSigm
 }
 
 Vector CoulombWithTensionCutOffImpl::DoReturnMapping(const Vector& rTrialSigmaTau,
-                                                     CoulombYieldSurface::CoulombAveragingType AveragingType)
+                                                     CoulombYieldSurface::CoulombAveragingType AveragingType,
+                                                     double& kappa_start)
 {
-    Vector     result      = ZeroVector(2);
-    const auto kappa_start = mCoulombYieldSurface.GetKappa();
+    Vector result = ZeroVector(2);
 
-    for (unsigned int counter = 0; counter < mCoulombYieldSurface.GetMaxIterations(); ++counter) {
-        const auto apex = mCoulombYieldSurface.CalculateApex();
+    if (kappa_start < 0.0) {
+        kappa_start = mCoulombYieldSurface.GetKappa();
+    } else {
+        mCoulombYieldSurface.SetKappa(kappa_start);
+    }
 
-        if (IsStressAtTensionApexReturnZone(rTrialSigmaTau, mTensionCutOff.GetTensileStrength(), apex)) {
-            return ReturnStressAtTensionApexReturnZone(mTensionCutOff.GetTensileStrength());
+    for (auto counter = std::size_t{0}; counter < mMaxNumberOfPlasticIterations; ++counter) {
+        if (IsStressAtTensionApexReturnZone(rTrialSigmaTau)) {
+            return ReturnStressAtTensionApexReturnZone();
         }
 
-        const auto corner_point = CalculateCornerPoint(mCoulombYieldSurface.GetFrictionAngleInRadians(),
-                                                       mCoulombYieldSurface.GetCohesion(),
-                                                       mTensionCutOff.GetTensileStrength(), apex);
-
-        if (IsStressAtTensionCutoffReturnZone(rTrialSigmaTau, mTensionCutOff.GetTensileStrength(),
-                                              apex, corner_point)) {
-            return ReturnStressAtTensionCutoffReturnZone(
-                rTrialSigmaTau, mTensionCutOff.DerivativeOfFlowFunction(rTrialSigmaTau),
-                mTensionCutOff.GetTensileStrength());
+        if (IsStressAtTensionCutoffReturnZone(rTrialSigmaTau)) {
+            return ReturnStressAtTensionCutoffReturnZone(rTrialSigmaTau);
         }
 
-        if (IsStressAtCornerReturnZone(
-                rTrialSigmaTau, mCoulombYieldSurface.DerivativeOfFlowFunction(rTrialSigmaTau, AveragingType),
-                corner_point)) {
-            result = corner_point;
+        if (IsStressAtCornerReturnZone(rTrialSigmaTau, AveragingType)) {
+            result = CalculateCornerPoint();
         } else { // Regular failure region
-            result = ReturnStressAtRegularFailureZone(
-                rTrialSigmaTau, mCoulombYieldSurface,
-                mCoulombYieldSurface.DerivativeOfFlowFunction(rTrialSigmaTau, AveragingType));
+            result = ReturnStressAtRegularFailureZone(rTrialSigmaTau, AveragingType);
         }
 
-        const auto lambda = mCoulombYieldSurface.CalculatePlasticMultiplier(
-            rTrialSigmaTau, mCoulombYieldSurface.DerivativeOfFlowFunction(rTrialSigmaTau, AveragingType));
-        const double kappa = kappa_start + mCoulombYieldSurface.CalculateEquivalentPlasticStrain(
-                                               rTrialSigmaTau, AveragingType, lambda);
+        const auto kappa = kappa_start + mCoulombYieldSurface.CalculateEquivalentPlasticStrainIncrement(
+                                             rTrialSigmaTau, AveragingType);
         mCoulombYieldSurface.SetKappa(kappa);
 
-        double error = std::abs(mCoulombYieldSurface.YieldFunctionValue(result));
-        if (error < mCoulombYieldSurface.GetConvergenceTolerance()) break;
+        if (std::abs(mCoulombYieldSurface.YieldFunctionValue(result)) < mAbsoluteYieldFunctionValueTolerance) {
+            break;
+        }
     }
     return result;
 }
 
+Vector CoulombWithTensionCutOffImpl::CalculateCornerPoint() const
+{
+    // 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;
+}
+
+bool CoulombWithTensionCutOffImpl::IsStressAtTensionApexReturnZone(const Vector& rTrialSigmaTau) const
+{
+    const auto tensile_strength = mTensionCutOff.GetTensileStrength();
+    return tensile_strength < mCoulombYieldSurface.CalculateApex() &&
+           rTrialSigmaTau[0] - rTrialSigmaTau[1] - tensile_strength > 0.0;
+}
+
+bool CoulombWithTensionCutOffImpl::IsStressAtTensionCutoffReturnZone(const Vector& rTrialSigmaTau) const
+{
+    const auto corner_point = CalculateCornerPoint();
+    return mTensionCutOff.GetTensileStrength() < mCoulombYieldSurface.CalculateApex() &&
+           corner_point[1] - rTrialSigmaTau[1] - corner_point[0] + rTrialSigmaTau[0] > 0.0;
+}
+
+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;
+}
+
+Vector CoulombWithTensionCutOffImpl::ReturnStressAtTensionApexReturnZone() const
+{
+    return UblasUtilities::CreateVector({mTensionCutOff.GetTensileStrength(), 0.0});
+}
+
+Vector CoulombWithTensionCutOffImpl::ReturnStressAtTensionCutoffReturnZone(const Vector& rSigmaTau) const
+{
+    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::ReturnStressAtRegularFailureZone(const Vector& rSigmaTau,
+                                                                      CoulombYieldSurface::CoulombAveragingType AveragingType) const
+{
+    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
 {
     rSerializer.save("CoulombYieldSurface", mCoulombYieldSurface);

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.h

@@ -32,11 +32,24 @@ class CoulombWithTensionCutOffImpl
 
     [[nodiscard]] bool   IsAdmissibleSigmaTau(const Vector& rTrialSigmaTau) const;
     [[nodiscard]] Vector DoReturnMapping(const Vector&                             rTrialSigmaTau,
-                                         CoulombYieldSurface::CoulombAveragingType AveragingType);
+                                         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;