Merge pull request #12414 from KratosMultiphysics/ImprovementsCLs

[ContitutiveLawsApp] Minor PR on iso/kin hardening plasticity integrators
KratosMultiphysics · Jun 3, 2024 · 3da6b5b · 3da6b5b
2 parents 4dabe00 + 888a094
commit 3da6b5b
Show file tree

Hide file tree

Showing 2 changed files with 28 additions and 131 deletions.
diff --git a/..._constitutive/auxiliary_files/cl_integrators/generic_cl_integrator_kinematic_plasticity.h b/..._constitutive/auxiliary_files/cl_integrators/generic_cl_integrator_kinematic_plasticity.h
@@ -108,7 +108,9 @@ class GenericConstitutiveLawIntegratorKinematicPlasticity
         ExponentialSoftening = 1,
         InitialHardeningExponentialSoftening = 2,
         PerfectPlasticity = 3,
-        CurveFittingHardening = 4
+        CurveFittingHardening = 4,
+        LinearExponentialSoftening = 5,
+        CurveDefinedByPoints = 6
     };
 
     enum class KinematicHardeningType
@@ -283,13 +285,13 @@ class GenericConstitutiveLawIntegratorKinematicPlasticity
         BoundedArrayType h_capa = ZeroVector(VoigtSize);
         double J2, tensile_indicator_factor, compression_indicator_factor, slope, hardening_parameter, equivalent_plastic_strain, I1;
 
-        YieldSurfaceType::CalculateEquivalentStress( rPredictiveStressVector, rStrainVector, rUniaxialStress, rValues);
+        YieldSurfaceType::CalculateEquivalentStress(rPredictiveStressVector, rStrainVector, rUniaxialStress, rValues);
         AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateI1Invariant(rPredictiveStressVector, I1);
         AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateJ2Invariant(rPredictiveStressVector, I1, deviator, J2);
         CalculateDerivativeYieldSurface(rPredictiveStressVector, deviator, J2, rYieldSurfaceDerivative, rValues);
         CalculateDerivativePlasticPotential(rPredictiveStressVector, deviator, J2, rDerivativePlasticPotential, rValues);
-        CalculateIndicatorsFactors(rPredictiveStressVector, tensile_indicator_factor,compression_indicator_factor);
-        CalculatePlasticDissipation(rPredictiveStressVector, tensile_indicator_factor,compression_indicator_factor, rPlasticStrainIncrement,rPlasticDissipation, h_capa, rValues, CharacteristicLength);
+        CalculateIndicatorsFactors(rPredictiveStressVector + rBackStressVector, tensile_indicator_factor,compression_indicator_factor); // The real stress state is used.
+        CalculatePlasticDissipation(rPredictiveStressVector + rBackStressVector, tensile_indicator_factor,compression_indicator_factor, rPlasticStrainIncrement,rPlasticDissipation, h_capa, rValues, CharacteristicLength); // The real stress state is used.
         CalculateEquivalentPlasticStrain(rPredictiveStressVector, rUniaxialStress, rPlasticStrain, tensile_indicator_factor, rValues, equivalent_plastic_strain);
         CalculateEquivalentStressThreshold(rPlasticDissipation, tensile_indicator_factor,compression_indicator_factor, rThreshold, slope, rValues, equivalent_plastic_strain, CharacteristicLength);
         CalculateHardeningParameter(rDerivativePlasticPotential, slope, h_capa, hardening_parameter);
@@ -437,8 +439,8 @@ class GenericConstitutiveLawIntegratorKinematicPlasticity
         )
     {
         GenericConstitutiveLawIntegratorPlasticity<YieldSurfaceType>::CalculatePlasticDissipation(
-            rPredictiveStressVector,TensileIndicatorFactor,CompressionIndicatorFactor,
-            PlasticStrainInc,rPlasticDissipation,rHCapa,rValues,CharacteristicLength);
+            rPredictiveStressVector, TensileIndicatorFactor, CompressionIndicatorFactor,
+            PlasticStrainInc, rPlasticDissipation, rHCapa, rValues, CharacteristicLength);
     }
 
     /**
@@ -464,124 +466,8 @@ class GenericConstitutiveLawIntegratorKinematicPlasticity
         )
     {
         GenericConstitutiveLawIntegratorPlasticity<YieldSurfaceType>::CalculateEquivalentStressThreshold(
-            PlasticDissipation,TensileIndicatorFactor,CompressionIndicatorFactor,rEquivalentStressThreshold,
-            rSlope,rValues,EquivalentPlasticStrain,CharacteristicLength);
-    }
-
-    /**
-     * @brief This method computes the uniaxial threshold using a linear softening
-     * @param PlasticDissipation The internal variable of energy dissipation due to plasticity
-     * @param TensileIndicatorFactor The tensile indicator
-     * @param CompressionIndicatorFactor The compressive indicator
-     * @param rEquivalentStressThreshold The maximum uniaxial stress of the linear behaviour
-     * @param rSlope The slope of the PlasticDiss-Threshold curve
-     * @param rValues Parameters of the constitutive law
-     */
-    static void CalculateEquivalentStressThresholdHardeningCurveLinearSoftening(
-        const double PlasticDissipation,
-        const double TensileIndicatorFactor,
-        const double CompressionIndicatorFactor,
-        double& rEquivalentStressThreshold,
-        double& rSlope,
-        ConstitutiveLaw::Parameters& rValues
-        )
-    {
-        GenericConstitutiveLawIntegratorPlasticity<YieldSurfaceType>::CalculateEquivalentStressThresholdHardeningCurveLinearSoftening(
-            PlasticDissipation,TensileIndicatorFactor,CompressionIndicatorFactor,rEquivalentStressThreshold,rSlope,rValues);
-    }
-
-    /**
-     * @brief This method computes the uniaxial threshold using a exponential softening
-     * @param PlasticDissipation The internal variable of energy dissipation due to plasticity
-     * @param TensileIndicatorFactor The tensile indicator
-     * @param CompressionIndicatorFactor The compressive indicator
-     * @param rEquivalentStressThreshold The maximum uniaxial stress of the linear behaviour
-     * @param rSlope The slope of the PlasticDiss-Threshold curve
-     * @param rValues Parameters of the constitutive law
-     */
-    static void CalculateEquivalentStressThresholdHardeningCurveExponentialSoftening(
-        const double PlasticDissipation,
-        const double TensileIndicatorFactor,
-        const double CompressionIndicatorFactor,
-        double& rEquivalentStressThreshold,
-        double& rSlope,
-        ConstitutiveLaw::Parameters& rValues,
-        const double CharacteristicLength
-        )
-    {
-        GenericConstitutiveLawIntegratorPlasticity<YieldSurfaceType>::CalculateEquivalentStressThresholdHardeningCurveExponentialSoftening(
-            PlasticDissipation,TensileIndicatorFactor,CompressionIndicatorFactor,rEquivalentStressThreshold,rSlope,rValues,CharacteristicLength);
-    }
-
-    /**
-     * @brief This method computes the uniaxial threshold using a hardening-softening law
-     * @param PlasticDissipation The internal variable of energy dissipation due to plasticity
-     * @param TensileIndicatorFactor The tensile indicator
-     * @param CompressionIndicatorFactor The compressive indicator
-     * @param rEquivalentStressThreshold The maximum uniaxial stress of the linear behaviour
-     * @param rSlope The slope of the PlasticDiss-Threshold curve
-     * @param rValues Parameters of the constitutive law
-     */
-    static void CalculateEquivalentStressThresholdHardeningCurveInitialHardeningExponentialSoftening(
-        const double PlasticDissipation,
-        const double TensileIndicatorFactor,
-        const double CompressionIndicatorFactor,
-        double& rEquivalentStressThreshold,
-        double& rSlope,
-        ConstitutiveLaw::Parameters& rValues
-        )
-    {
-        GenericConstitutiveLawIntegratorPlasticity<YieldSurfaceType>::CalculateEquivalentStressThresholdHardeningCurveInitialHardeningExponentialSoftening(
-            PlasticDissipation,TensileIndicatorFactor,CompressionIndicatorFactor,rEquivalentStressThreshold,rSlope,rValues);
-    }
-
-    /**
-     * @brief This method computes the uniaxial threshold using a perfect plasticity law
-     * @param PlasticDissipation The internal variable of energy dissipation due to plasticity
-     * @param TensileIndicatorFactor The tensile indicator
-     * @param CompressionIndicatorFactor The compressive indicator
-     * @param rEquivalentStressThreshold The maximum uniaxial stress of the linear behaviour
-     * @param rSlope The slope of the PlasticDiss-Threshold curve
-     * @param rValues Parameters of the constitutive law
-     */
-    static void CalculateEquivalentStressThresholdHardeningCurvePerfectPlasticity(
-        const double PlasticDissipation,
-        const double TensileIndicatorFactor,
-        const double CompressionIndicatorFactor,
-        double& rEquivalentStressThreshold,
-        double& rSlope,
-        ConstitutiveLaw::Parameters& rValues
-        )
-    {
-        GenericConstitutiveLawIntegratorPlasticity<YieldSurfaceType>::CalculateEquivalentStressThresholdHardeningCurvePerfectPlasticity(
-            PlasticDissipation,TensileIndicatorFactor,CompressionIndicatorFactor,rEquivalentStressThreshold,rSlope,rValues);
-    }
-
-    /**
-     * @brief This method computes the uniaxial threshold using a perfect plasticity law
-     * @param PlasticDissipation The internal variable of energy dissipation due to plasticity
-     * @param TensileIndicatorFactor The tensile indicator
-     * @param CompressionIndicatorFactor The compressive indicator
-     * @param rEquivalentStressThreshold The maximum uniaxial stress of the linear behaviour
-     * @param rSlope The slope of the PlasticDiss-Threshold curve
-     * @param rValues Parameters of the constitutive law
-     * @param EquivalentPlasticStrain The Plastic Strain internal variable
-     * @param CharacteristicLength Characteristic length of the finite element
-     */
-    static void CalculateEquivalentStressThresholdCurveFittingHardening(
-        const double PlasticDissipation,
-        const double TensileIndicatorFactor,
-        const double CompressionIndicatorFactor,
-        double& rEquivalentStressThreshold,
-        double& rSlope,
-        ConstitutiveLaw::Parameters& rValues,
-        const double EquivalentPlasticStrain,
-        const double CharacteristicLength
-        )
-    {
-        GenericConstitutiveLawIntegratorPlasticity<YieldSurfaceType>::CalculateEquivalentStressThresholdCurveFittingHardening(
-            PlasticDissipation,TensileIndicatorFactor,CompressionIndicatorFactor,rEquivalentStressThreshold,rSlope,rValues,
-            EquivalentPlasticStrain,CharacteristicLength);
+            PlasticDissipation, TensileIndicatorFactor, CompressionIndicatorFactor, rEquivalentStressThreshold,
+            rSlope, rValues, EquivalentPlasticStrain, CharacteristicLength);
     }
 
     /**
@@ -603,7 +489,7 @@ class GenericConstitutiveLawIntegratorKinematicPlasticity
         )
     {
         GenericConstitutiveLawIntegratorPlasticity<YieldSurfaceType>::CalculateEquivalentPlasticStrain(
-            rStressVector,UniaxialStress,rPlasticStrain,r0,rValues,rEquivalentPlasticStrain);
+            rStressVector, UniaxialStress, rPlasticStrain, r0, rValues, rEquivalentPlasticStrain);
     }
 
     /**
@@ -631,7 +517,7 @@ class GenericConstitutiveLawIntegratorKinematicPlasticity
         )
     {
         GenericConstitutiveLawIntegratorPlasticity<YieldSurfaceType>::CalculateHardeningParameter(
-            rGFlux,SlopeThreshold,rHCapa,rHardeningParameter);
+            rGFlux, SlopeThreshold, rHCapa, rHardeningParameter);
     }
 
     /**

diff --git a/...ion/custom_constitutive/auxiliary_files/cl_integrators/generic_cl_integrator_plasticity.h b/...ion/custom_constitutive/auxiliary_files/cl_integrators/generic_cl_integrator_plasticity.h
@@ -247,7 +247,7 @@ class GenericConstitutiveLawIntegratorPlasticity
         array_1d<double, VoigtSize> h_capa = ZeroVector(VoigtSize);
         double J2, tensile_indicator_factor, compression_indicator_factor, slope, hardening_parameter, equivalent_plastic_strain, I1;
 
-        YieldSurfaceType::CalculateEquivalentStress( rPredictiveStressVector, rStrainVector, rUniaxialStress, rValues);
+        YieldSurfaceType::CalculateEquivalentStress(rPredictiveStressVector, rStrainVector, rUniaxialStress, rValues);
         AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateI1Invariant(rPredictiveStressVector, I1);
         AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateJ2Invariant(rPredictiveStressVector, I1, deviator, J2);
         CalculateFFluxVector(rPredictiveStressVector, deviator, J2, rFflux, rValues);
@@ -586,11 +586,22 @@ class GenericConstitutiveLawIntegratorPlasticity
 
         const double minimum_characteristic_fracture_energy_exponential_softening = (std::pow(yield_compression, 2)) / young_modulus;
 
+        const bool has_total_or_plastic_strain_space = r_material_properties.Has(TOTAL_OR_PLASTIC_STRAIN_SPACE);
+        const bool total_or_plastic_strain_space = has_total_or_plastic_strain_space ? r_material_properties[TOTAL_OR_PLASTIC_STRAIN_SPACE] : false; //Default value = plastic strain space
+
         double initial_threshold;
         GetInitialUniaxialThreshold(rValues, initial_threshold);
         KRATOS_ERROR_IF(characteristic_fracture_energy_compression < minimum_characteristic_fracture_energy_exponential_softening) << "The Fracture Energy is to low: " << characteristic_fracture_energy_compression << std::endl;
-        rEquivalentStressThreshold = initial_threshold * (1.0 - PlasticDissipation);
-        rSlope = - initial_threshold;
+        if (total_or_plastic_strain_space) { // Curve built in the total strain space
+            rEquivalentStressThreshold = (young_modulus / initial_threshold) * ((0.5 * std::pow(initial_threshold, 2.0) / young_modulus - characteristic_fracture_energy_compression)
+                                            + std::sqrt(std::pow((0.5 * std::pow(initial_threshold, 2.0) / young_modulus + characteristic_fracture_energy_compression), 2.0)
+                                            - 2.0 * std::pow(initial_threshold, 2.0) / young_modulus * characteristic_fracture_energy_compression * PlasticDissipation));
+            rSlope = - initial_threshold * characteristic_fracture_energy_compression * std::pow(std::pow((0.5 * std::pow(initial_threshold, 2.0) / young_modulus + characteristic_fracture_energy_compression), 2.0)
+                    - 2.0 * std::pow(initial_threshold, 2.0) / young_modulus * characteristic_fracture_energy_compression * PlasticDissipation, -0.5);
+        } else { // Curve built in the plastic strain space
+            rEquivalentStressThreshold = initial_threshold * (1.0 - PlasticDissipation);
+            rSlope = - initial_threshold;
+        }
     }
 
     /**
@@ -658,7 +669,7 @@ class GenericConstitutiveLawIntegratorPlasticity
     }
 
     /**
-     * @brief This method computes the uniaxial threshold using a perfect plasticity law
+     * @brief This method computes the uniaxial threshold using an isotropic hardening (polynomial) - softening (exponential) plasticity law
      * @param PlasticDissipation The internal variable of energy dissipation due to plasticity
      * @param TensileIndicatorFactor The tensile indicator
      * @param CompressionIndicatorFactor The compressive indicator
@@ -761,7 +772,7 @@ class GenericConstitutiveLawIntegratorPlasticity
     }
 
     /**
-     * @brief This method computes the uniaxial threshold using a linear-exponential softening, which changes from one to the other through the platic_dissipation_limit
+     * @brief This method computes the uniaxial threshold using an isotropic linear/exponential softening, which changes from one to the other through the platic_dissipation_limit
      * @param PlasticDissipation The internal variable of energy dissipation due to plasticity
      * @param TensileIndicatorFactor The tensile indicator
      * @param CompressionIndicatorFactor The compressive indicator