Skip to content

[ConstitutiveLawsApplication] Double Exponential Hardening Curve for Isotropic Damage Constitutive Law #14280

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
aslan1372-droid wants to merge 9 commits into
base: master
Choose a base branch
from
Open

[ConstitutiveLawsApplication] Double Exponential Hardening Curve for Isotropic Damage Constitutive Law #14280

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
e257419
Creating a new hardening curve for isotropic damage CL
aslan1372-droid Mar 2, 2026
c14d8d3
Merge remote-tracking branch 'origin/master'
aslan1372-droid Mar 2, 2026
33bb0b9
Merge branch 'master' into DoubleExponentialHardCurveIsoDamCL
aslan1372-droid Mar 2, 2026
15f60b4
Merge branch 'DoubleExponentialHardCurveIsoDamCL'
aslan1372-droid Mar 2, 2026
c5e2fc7
Merge branch 'master' into DoubleExponentialHardCurveIsoDamCL
aslan1372-droid Mar 2, 2026
ac7e149
Revert "Creating a new hardening curve for isotropic damage CL"
aslan1372-droid Mar 2, 2026
109da17
Merge branch 'DoubleExponentialHardCurveIsoDamCL'
aslan1372-droid Mar 2, 2026
c4e901c
Revert "Merge branch 'DoubleExponentialHardCurveIsoDamCL'"
aslan1372-droid Mar 2, 2026
caca1fd
Modifying description
Mar 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions applications/ConstitutiveLawsApplication/constitutive_laws_application.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
1 change: 1 addition & 0 deletions applications/ConstitutiveLawsApplication/constitutive_laws_application_variables.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
3 changes: 2 additions & 1 deletion applications/ConstitutiveLawsApplication/constitutive_laws_application_variables.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down Expand Up @@ -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)
Expand Down
43 changes: 43 additions & 0 deletions ...ication/custom_constitutive/auxiliary_files/cl_integrators/generic_cl_integrator_damage.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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];
Copy link
Copy Markdown
Member

@AlejandroCornejo AlejandroCornejo Mar 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you can check if this variable is present in the properties and assign a default if it is not, more convenient. Something like:

const bool satisfy_max_stress = r_mat_props.Has(SATISFY_MAXIMUM_STRESS) ? r_mat_props[SATISFY_MAXIMUM_STRESS] : true;

Copy link
Copy Markdown
Member

@loumalouomega loumalouomega Mar 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
const bool Max_Stress_Satisfied =r_mat_props[SATISFY_MAXIMUM_STRESS];
const bool max_stress_satisfied =r_mat_props[SATISFY_MAXIMUM_STRESS];

To follow the style


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));
Copy link
Copy Markdown
Member

@AlejandroCornejo AlejandroCornejo Mar 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we must check in advance that the denominator is not 0 to avoid nan

} 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));
Copy link
Copy Markdown
Member

@AlejandroCornejo AlejandroCornejo Mar 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please make sure of format the code so you add all spaces between operations, it improves readability. besides better use 1.0.


}

/**
* @brief This computes the damage variable according to linear softening
* @param UniaxialStress The equivalent uniaxial stress
Expand Down
2 changes: 2 additions & 0 deletions ...ations/ConstitutiveLawsApplication/custom_python/constitutive_laws_python_application.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading