diff --git a/applications/ConstitutiveLawsApplication/constitutive_laws_application.cpp b/applications/ConstitutiveLawsApplication/constitutive_laws_application.cpp
index 4d06bc1bf008..dd15c2fb2438 100644
--- a/applications/ConstitutiveLawsApplication/constitutive_laws_application.cpp
+++ b/applications/ConstitutiveLawsApplication/constitutive_laws_application.cpp
@@ -492,6 +492,7 @@ void KratosConstitutiveLawsApplication::Register()
     KRATOS_REGISTER_VARIABLE(DELAY_TIME)
     KRATOS_REGISTER_VARIABLE(MAXIMUM_STRESS)
     KRATOS_REGISTER_VARIABLE(MAXIMUM_STRESS_POSITION)
+    KRATOS_REGISTER_VARIABLE(SATISFY_MAXIMUM_STRESS)
     KRATOS_REGISTER_VARIABLE(PLASTIC_DISSIPATION_LIMIT_LINEAR_SOFTENING)
     KRATOS_REGISTER_VARIABLE(UNIAXIAL_STRESS)
     KRATOS_REGISTER_VARIABLE(FRICTION_ANGLE)
diff --git a/applications/ConstitutiveLawsApplication/constitutive_laws_application_variables.cpp b/applications/ConstitutiveLawsApplication/constitutive_laws_application_variables.cpp
index 0497f559cbb5..b1cc6749ec9f 100644
--- a/applications/ConstitutiveLawsApplication/constitutive_laws_application_variables.cpp
+++ b/applications/ConstitutiveLawsApplication/constitutive_laws_application_variables.cpp
@@ -75,6 +75,7 @@ KRATOS_CREATE_VARIABLE(double, VISCOUS_PARAMETER)
 KRATOS_CREATE_VARIABLE(double, DELAY_TIME)
 KRATOS_CREATE_VARIABLE(double, MAXIMUM_STRESS)
 KRATOS_CREATE_VARIABLE(double, MAXIMUM_STRESS_POSITION)
+KRATOS_CREATE_VARIABLE(bool, SATISFY_MAXIMUM_STRESS)
 KRATOS_CREATE_VARIABLE(double, PLASTIC_DISSIPATION_LIMIT_LINEAR_SOFTENING)
 KRATOS_CREATE_VARIABLE(double, UNIAXIAL_STRESS)
 KRATOS_CREATE_VARIABLE(double, FRICTION_ANGLE)
diff --git a/applications/ConstitutiveLawsApplication/constitutive_laws_application_variables.h b/applications/ConstitutiveLawsApplication/constitutive_laws_application_variables.h
index 5cb53e46f93e..1ed4c6961aaf 100644
--- a/applications/ConstitutiveLawsApplication/constitutive_laws_application_variables.h
+++ b/applications/ConstitutiveLawsApplication/constitutive_laws_application_variables.h
@@ -24,7 +24,7 @@
 
 namespace Kratos
 {
-    enum class SofteningType {Linear = 0, Exponential = 1, HardeningDamage = 2, CurveFittingDamage = 3};
+    enum class SofteningType {Linear = 0, Exponential = 1, HardeningDamage = 2, CurveFittingDamage = 3, DoubleExponentialHardeningDamage = 4};
 
     enum class TangentOperatorEstimation {Analytic = 0, FirstOrderPerturbation = 1, SecondOrderPerturbation = 2,
         Secant = 3, SecondOrderPerturbationV2 = 4, InitialStiffness = 5, OrthogonalSecant = 6};
@@ -93,6 +93,7 @@ namespace Kratos
     KRATOS_DEFINE_APPLICATION_VARIABLE(CONSTITUTIVE_LAWS_APPLICATION, double, DELAY_TIME)
     KRATOS_DEFINE_APPLICATION_VARIABLE(CONSTITUTIVE_LAWS_APPLICATION, double, MAXIMUM_STRESS)
     KRATOS_DEFINE_APPLICATION_VARIABLE(CONSTITUTIVE_LAWS_APPLICATION, double, MAXIMUM_STRESS_POSITION)
+    KRATOS_DEFINE_APPLICATION_VARIABLE(CONSTITUTIVE_LAWS_APPLICATION, bool, SATISFY_MAXIMUM_STRESS)
     KRATOS_DEFINE_APPLICATION_VARIABLE(CONSTITUTIVE_LAWS_APPLICATION, double, PLASTIC_DISSIPATION_LIMIT_LINEAR_SOFTENING)
     KRATOS_DEFINE_APPLICATION_VARIABLE(CONSTITUTIVE_LAWS_APPLICATION, double, UNIAXIAL_STRESS)
     KRATOS_DEFINE_APPLICATION_VARIABLE(CONSTITUTIVE_LAWS_APPLICATION, double, FRICTION_ANGLE)
diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/auxiliary_files/cl_integrators/generic_cl_integrator_damage.h b/applications/ConstitutiveLawsApplication/custom_constitutive/auxiliary_files/cl_integrators/generic_cl_integrator_damage.h
index f2e878aeb958..6cd2902300ba 100644
--- a/applications/ConstitutiveLawsApplication/custom_constitutive/auxiliary_files/cl_integrators/generic_cl_integrator_damage.h
+++ b/applications/ConstitutiveLawsApplication/custom_constitutive/auxiliary_files/cl_integrators/generic_cl_integrator_damage.h
@@ -151,6 +151,9 @@ class GenericConstitutiveLawIntegratorDamage
         case static_cast<int>(SofteningType::CurveFittingDamage):
             CalculateCurveFittingDamage(UniaxialStress, rThreshold, damage_parameter, CharacteristicLength, rValues, rDamage);
             break;
+        case static_cast<int>(SofteningType::DoubleExponentialHardeningDamage):
+            CalculateDoubleExponentialHardeningDamage(UniaxialStress, rThreshold, damage_parameter, CharacteristicLength, rValues, rDamage);
+            break;
         default:
             KRATOS_ERROR << "SOFTENING_TYPE not defined or wrong..." << softening_type << std::endl;
             break;
@@ -228,6 +231,46 @@ class GenericConstitutiveLawIntegratorDamage
         }
     }
 
+    /**
+     * @brief This computes the damage variable according to combined exponential hardening and softening
+     * @param UniaxialStress The equivalent uniaxial stress
+     * @param Threshold The maximum uniaxial stress achieved previously
+     * @param rDamage The internal variable of the damage model
+     * @param rValues Parameters of the constitutive law
+     * @param CharacteristicLength The equivalent length of the FE
+     * @param SATISFY_MAXIMUM_STRESS Select either maximum stress or tangency at yield point to be satisfied
+     */
+    static void CalculateDoubleExponentialHardeningDamage(
+        const double UniaxialStress,
+        const double Threshold,
+        const double DamageParameter,
+        const double CharacteristicLength,
+        ConstitutiveLaw::Parameters& rValues,
+        double& rDamage
+        )
+    {
+        const auto &r_mat_props = rValues.GetMaterialProperties();
+        const double max_stress = r_mat_props[MAXIMUM_STRESS];
+        const double Gf = r_mat_props[FRACTURE_ENERGY];
+        const double E = r_mat_props[YOUNG_MODULUS];
+        const bool Max_Stress_Satisfied =r_mat_props[SATISFY_MAXIMUM_STRESS];
+
+        double initial_threshold;
+        TYieldSurfaceType::GetInitialUniaxialThreshold(rValues, initial_threshold);
+
+        double X = 0.0;
+        if (Max_Stress_Satisfied) {
+            X = -std::sqrt(max_stress/(max_stress - initial_threshold));
+        } else {
+            X = (initial_threshold*initial_threshold/E + Gf + std::sqrt(initial_threshold*initial_threshold/E*(5.0/4.0*initial_threshold*initial_threshold/E+2.0*Gf))) / (0.5*initial_threshold*initial_threshold/E-Gf);
+        }
+
+        const double A = ((3.0*X+1.0)*initial_threshold*(UniaxialStress-initial_threshold))/((X+1.0)*(0.5*initial_threshold*initial_threshold-E*Gf));
+        
+        rDamage = 1-initial_threshold/(UniaxialStress*(X+1.0)) * (2.0*X*std::exp(A/2.0)+(1.0-X)*std::exp(A));
+
+    }
+
     /**
      * @brief This computes the damage variable according to linear softening
      * @param UniaxialStress The equivalent uniaxial stress
diff --git a/applications/ConstitutiveLawsApplication/custom_python/constitutive_laws_python_application.cpp b/applications/ConstitutiveLawsApplication/custom_python/constitutive_laws_python_application.cpp
index 0e1b74d88a4f..66ba5716604a 100644
--- a/applications/ConstitutiveLawsApplication/custom_python/constitutive_laws_python_application.cpp
+++ b/applications/ConstitutiveLawsApplication/custom_python/constitutive_laws_python_application.cpp
@@ -171,6 +171,8 @@ PYBIND11_MODULE(KratosConstitutiveLawsApplication,m)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, MODE_TWO_FRACTURE_ENERGY)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, TENSILE_INTERFACE_MODULUS)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, SHEAR_INTERFACE_MODULUS)
+    KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, MAXIMUM_STRESS)
+    KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, SATISFY_MAXIMUM_STRESS)
 
 
     // D+D- Damage Constitutive laws variables, additional Masonry 2D & 3D