diff --git a/applications/DEMApplication/CMakeLists.txt b/applications/DEMApplication/CMakeLists.txt index a967659baa3d..ac0ce65062ee 100755 --- a/applications/DEMApplication/CMakeLists.txt +++ b/applications/DEMApplication/CMakeLists.txt @@ -91,6 +91,8 @@ set(KRATOS_DEM_APPLICATION_CORE ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_parallel_bond_CL.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_smooth_joint_CL.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_parallel_bond_for_membrane_CL.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_parallel_bond_bilinear_damage_CL.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_parallel_bond_bilinear_damage_mixed_CL.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/inlet.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/fast_filling_creator.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/force_based_inlet.cpp diff --git a/applications/DEMApplication/DEM_application.cpp b/applications/DEMApplication/DEM_application.cpp index 6342e52721f6..72eaf8dc3d09 100644 --- a/applications/DEMApplication/DEM_application.cpp +++ b/applications/DEMApplication/DEM_application.cpp @@ -49,6 +49,8 @@ #include "custom_constitutive/DEM_parallel_bond_CL.h" #include "custom_constitutive/DEM_smooth_joint_CL.h" #include "custom_constitutive/DEM_parallel_bond_for_membrane_CL.h" +#include "custom_constitutive/DEM_parallel_bond_bilinear_damage_CL.h" +#include "custom_constitutive/DEM_parallel_bond_bilinear_damage_mixed_CL.h" #include "custom_constitutive/DEM_rolling_friction_model.h" #include "custom_constitutive/DEM_rolling_friction_model_constant_torque.h" #include "custom_constitutive/DEM_rolling_friction_model_viscous_torque.h" @@ -107,6 +109,7 @@ KRATOS_CREATE_VARIABLE(bool, GLOBAL_COORDINATION_NUMBER_OPTION) KRATOS_CREATE_VARIABLE(bool, AUTOMATIC_SKIN_COMPUTATION) KRATOS_CREATE_VARIABLE(double, SKIN_FACTOR_RADIUS) KRATOS_CREATE_VARIABLE(int, CLEAN_INDENT_OPTION) +KRATOS_CREATE_VARIABLE(int, CLEAN_INDENT_V2_OPTION) KRATOS_CREATE_VARIABLE(int, TRIHEDRON_OPTION) KRATOS_CREATE_VARIABLE(int, ROLLING_FRICTION_OPTION) KRATOS_CREATE_VARIABLE(int, POISSON_EFFECT_OPTION) @@ -134,6 +137,7 @@ KRATOS_CREATE_VARIABLE(double, DEM_ENGINE_PERFORMANCE) KRATOS_CREATE_VARIABLE(double, DEM_DRAG_CONSTANT_X) KRATOS_CREATE_VARIABLE(double, DEM_DRAG_CONSTANT_Y) KRATOS_CREATE_VARIABLE(double, DEM_DRAG_CONSTANT_Z) +KRATOS_CREATE_VARIABLE(bool, ENERGY_CALCULATION_OPTION) KRATOS_CREATE_VARIABLE(double, INITIAL_VELOCITY_X_VALUE) KRATOS_CREATE_VARIABLE(double, INITIAL_VELOCITY_Y_VALUE) @@ -253,6 +257,9 @@ KRATOS_CREATE_VARIABLE(bool, DEBUG_PRINTING_OPTION) KRATOS_CREATE_VARIABLE(int, DEBUG_PRINTING_ID_1) KRATOS_CREATE_VARIABLE(int, DEBUG_PRINTING_ID_2) KRATOS_CREATE_VARIABLE(double, FRACTURE_ENERGY) +KRATOS_CREATE_VARIABLE(double, FRACTURE_ENERGY_NORMAL) +KRATOS_CREATE_VARIABLE(double, FRACTURE_ENERGY_TANGENTIAL) +KRATOS_CREATE_VARIABLE(double, FRACTURE_ENERGY_EXPONENT) KRATOS_CREATE_VARIABLE(double, SIGMA_SLOPE_CHANGE_THRESHOLD) KRATOS_CREATE_VARIABLE(double, INTERNAL_FRICTION_AFTER_THRESHOLD) KRATOS_CREATE_VARIABLE(double, SEARCH_RADIUS_INCREMENT_FOR_BONDS_CREATION) @@ -325,6 +332,12 @@ KRATOS_CREATE_VARIABLE(Quaternion, AUX_ORIENTATION) KRATOS_CREATE_3D_VARIABLE_WITH_COMPONENTS(LOCAL_AUX_ANGULAR_VELOCITY) // ******************* Quaternion Integration END ******************* +// ****************Radius expansion method BEGIN******************* +KRATOS_CREATE_VARIABLE(bool, IS_RADIUS_EXPANSION) +KRATOS_CREATE_VARIABLE(double, RADIUS_EXPANSION_RATE) +KRATOS_CREATE_VARIABLE(double, RADIUS_MULTIPLIER_MAX) +// *****************Radius expansion method END******************** + // FORCE AND MOMENTUM KRATOS_CREATE_3D_VARIABLE_WITH_COMPONENTS(CONTACT_IMPULSE) KRATOS_CREATE_3D_VARIABLE_WITH_COMPONENTS(PARTICLE_MOMENT) @@ -347,6 +360,7 @@ KRATOS_CREATE_VARIABLE(double, PARTICLE_GRAVITATIONAL_ENERGY) KRATOS_CREATE_VARIABLE(double, PARTICLE_INELASTIC_VISCODAMPING_ENERGY) KRATOS_CREATE_VARIABLE(double, PARTICLE_INELASTIC_FRICTIONAL_ENERGY) KRATOS_CREATE_VARIABLE(double, PARTICLE_INELASTIC_ROLLING_RESISTANCE_ENERGY) +KRATOS_CREATE_VARIABLE(double, PARTICLE_MAX_NORMAL_BALL_TO_BALL_FORCE_TIMES_RADIUS) KRATOS_CREATE_VARIABLE(int, COMPUTE_ENERGY_OPTION) KRATOS_CREATE_VARIABLE(double, GLOBAL_DAMPING) KRATOS_CREATE_VARIABLE(double, GLOBAL_VISCOUS_DAMPING) @@ -539,6 +553,14 @@ template class DEM_compound_constitutive_law_for_PBM; template class DEM_compound_constitutive_law_for_PBM; +template class DEM_compound_constitutive_law_for_PBM; +template class DEM_compound_constitutive_law_for_PBM; +template class DEM_compound_constitutive_law_for_PBM; + +template class DEM_compound_constitutive_law_for_PBM; +template class DEM_compound_constitutive_law_for_PBM; +template class DEM_compound_constitutive_law_for_PBM; + void KratosDEMApplication::Register() { KRATOS_INFO("") << "\n" << " KRATOS | _ \\| ____| \\/ | _ \\ __ _ ___| | __ \n" @@ -593,6 +615,7 @@ void KratosDEMApplication::Register() { KRATOS_REGISTER_VARIABLE(AUTOMATIC_SKIN_COMPUTATION) KRATOS_REGISTER_VARIABLE(SKIN_FACTOR_RADIUS) KRATOS_REGISTER_VARIABLE(CLEAN_INDENT_OPTION) + KRATOS_REGISTER_VARIABLE(CLEAN_INDENT_V2_OPTION) KRATOS_REGISTER_VARIABLE(TRIHEDRON_OPTION) KRATOS_REGISTER_VARIABLE(NEIGH_INITIALIZED) KRATOS_REGISTER_VARIABLE(TRIAXIAL_TEST_OPTION) @@ -616,6 +639,7 @@ void KratosDEMApplication::Register() { KRATOS_REGISTER_VARIABLE(DEM_DRAG_CONSTANT_X) KRATOS_REGISTER_VARIABLE(DEM_DRAG_CONSTANT_Y) KRATOS_REGISTER_VARIABLE(DEM_DRAG_CONSTANT_Z) + KRATOS_REGISTER_VARIABLE(ENERGY_CALCULATION_OPTION) KRATOS_REGISTER_VARIABLE(INITIAL_VELOCITY_X_VALUE) KRATOS_REGISTER_VARIABLE(INITIAL_VELOCITY_Y_VALUE) @@ -729,6 +753,9 @@ void KratosDEMApplication::Register() { KRATOS_REGISTER_VARIABLE(DEBUG_PRINTING_ID_1) KRATOS_REGISTER_VARIABLE(DEBUG_PRINTING_ID_2) KRATOS_REGISTER_VARIABLE(FRACTURE_ENERGY) + KRATOS_REGISTER_VARIABLE(FRACTURE_ENERGY_NORMAL) + KRATOS_REGISTER_VARIABLE(FRACTURE_ENERGY_TANGENTIAL) + KRATOS_REGISTER_VARIABLE(FRACTURE_ENERGY_EXPONENT) KRATOS_REGISTER_VARIABLE(SIGMA_SLOPE_CHANGE_THRESHOLD) KRATOS_REGISTER_VARIABLE(INTERNAL_FRICTION_AFTER_THRESHOLD) KRATOS_REGISTER_VARIABLE(SEARCH_RADIUS_INCREMENT_FOR_BONDS_CREATION) @@ -808,6 +835,12 @@ void KratosDEMApplication::Register() { KRATOS_REGISTER_VARIABLE(LOCAL_AUX_ANGULAR_VELOCITY) // ******************* Quaternion Integration END ******************* + // ****************Radius expansion method BEGIN******************* + KRATOS_REGISTER_VARIABLE(IS_RADIUS_EXPANSION) + KRATOS_REGISTER_VARIABLE(RADIUS_EXPANSION_RATE) + KRATOS_REGISTER_VARIABLE(RADIUS_MULTIPLIER_MAX) + // *****************Radius expansion method END******************** + // FORCE AND MOMENTUM KRATOS_REGISTER_3D_VARIABLE_WITH_COMPONENTS(CONTACT_IMPULSE) KRATOS_REGISTER_3D_VARIABLE_WITH_COMPONENTS(PARTICLE_MOMENT) @@ -830,6 +863,7 @@ void KratosDEMApplication::Register() { KRATOS_REGISTER_VARIABLE(PARTICLE_INELASTIC_VISCODAMPING_ENERGY) KRATOS_REGISTER_VARIABLE(PARTICLE_INELASTIC_FRICTIONAL_ENERGY) KRATOS_REGISTER_VARIABLE(PARTICLE_INELASTIC_ROLLING_RESISTANCE_ENERGY) + KRATOS_REGISTER_VARIABLE(PARTICLE_MAX_NORMAL_BALL_TO_BALL_FORCE_TIMES_RADIUS) KRATOS_REGISTER_VARIABLE(COMPUTE_ENERGY_OPTION) KRATOS_REGISTER_VARIABLE(GLOBAL_DAMPING) KRATOS_REGISTER_VARIABLE(GLOBAL_VISCOUS_DAMPING) @@ -1026,6 +1060,8 @@ void KratosDEMApplication::Register() { Serializer::Register("DEM_parallel_bond", DEM_parallel_bond()); Serializer::Register("DEM_smooth_joint", DEM_smooth_joint()); Serializer::Register("DEM_parallel_bond_for_membrane", DEM_parallel_bond_for_membrane()); + Serializer::Register("DEM_parallel_bond_bilinear_damage", DEM_parallel_bond_bilinear_damage()); + Serializer::Register("DEM_parallel_bond_bilinear_damage_mixed", DEM_parallel_bond_bilinear_damage_mixed()); Serializer::Register("DEMRollingFrictionModelConstantTorque", DEMRollingFrictionModelConstantTorque()); Serializer::Register("DEMRollingFrictionModelViscousTorque", DEMRollingFrictionModelViscousTorque()); Serializer::Register("DEMRollingFrictionModelBounded", DEMRollingFrictionModelBounded()); diff --git a/applications/DEMApplication/DEM_application_variables.h b/applications/DEMApplication/DEM_application_variables.h index 33c934f34b6f..baf94b2340d5 100644 --- a/applications/DEMApplication/DEM_application_variables.h +++ b/applications/DEMApplication/DEM_application_variables.h @@ -76,6 +76,7 @@ namespace Kratos KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, bool, AUTOMATIC_SKIN_COMPUTATION) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, SKIN_FACTOR_RADIUS) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, CLEAN_INDENT_OPTION) + KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, CLEAN_INDENT_V2_OPTION) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, TRIHEDRON_OPTION) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, ROLLING_FRICTION_OPTION) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, POISSON_EFFECT_OPTION) @@ -102,6 +103,7 @@ namespace Kratos KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, DEM_DRAG_CONSTANT_X) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, DEM_DRAG_CONSTANT_Y) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, DEM_DRAG_CONSTANT_Z) + KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, bool, ENERGY_CALCULATION_OPTION) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, INITIAL_VELOCITY_X_VALUE) @@ -217,6 +219,9 @@ namespace Kratos KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, DEBUG_PRINTING_ID_1) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, DEBUG_PRINTING_ID_2) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, FRACTURE_ENERGY) + KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, FRACTURE_ENERGY_NORMAL) + KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, FRACTURE_ENERGY_TANGENTIAL) + KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, FRACTURE_ENERGY_EXPONENT) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, SIGMA_SLOPE_CHANGE_THRESHOLD) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, INTERNAL_FRICTION_AFTER_THRESHOLD) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, SEARCH_RADIUS_INCREMENT_FOR_BONDS_CREATION) @@ -301,6 +306,12 @@ namespace Kratos KRATOS_DEFINE_3D_APPLICATION_VARIABLE_WITH_COMPONENTS(DEM_APPLICATION, LOCAL_AUX_ANGULAR_VELOCITY) // ******************* Quaternion Integration END ******************* + // ****************Radius expansion method BEGIN******************* + KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, bool, IS_RADIUS_EXPANSION) + KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, RADIUS_EXPANSION_RATE) + KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, RADIUS_MULTIPLIER_MAX) + // *****************Radius expansion method END******************** + // FORCE AND MOMENTUM KRATOS_DEFINE_3D_APPLICATION_VARIABLE_WITH_COMPONENTS(DEM_APPLICATION, PARTICLE_MOMENT) @@ -324,6 +335,7 @@ namespace Kratos KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, PARTICLE_INELASTIC_VISCODAMPING_ENERGY) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, PARTICLE_INELASTIC_FRICTIONAL_ENERGY) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, PARTICLE_INELASTIC_ROLLING_RESISTANCE_ENERGY) + KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, PARTICLE_MAX_NORMAL_BALL_TO_BALL_FORCE_TIMES_RADIUS) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, COMPUTE_ENERGY_OPTION) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, GLOBAL_DAMPING) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, GLOBAL_VISCOUS_DAMPING) diff --git a/applications/DEMApplication/custom_constitutive/DEM_D_Hertz_viscous_Coulomb_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_D_Hertz_viscous_Coulomb_CL.cpp index 07a290d88426..1150c9a0ec29 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_D_Hertz_viscous_Coulomb_CL.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_D_Hertz_viscous_Coulomb_CL.cpp @@ -161,6 +161,7 @@ namespace Kratos { ViscoDampingLocalContactForce[0] = - equiv_visco_damp_coeff_tangential * LocalRelVel[0]; ViscoDampingLocalContactForce[1] = - equiv_visco_damp_coeff_tangential * LocalRelVel[1]; ViscoDampingLocalContactForce[2] = - equiv_visco_damp_coeff_normal * LocalRelVel[2]; + } ///////////////////////// @@ -392,4 +393,8 @@ namespace Kratos { inelastic_viscodamping_energy += viscodamping_energy; } + double DEM_D_Hertz_viscous_Coulomb::GetTangentialStiffness() { + return mKt; + } + } // namespace Kratos diff --git a/applications/DEMApplication/custom_constitutive/DEM_D_Hertz_viscous_Coulomb_CL.h b/applications/DEMApplication/custom_constitutive/DEM_D_Hertz_viscous_Coulomb_CL.h index 7ae4116895b9..71a9d10bbab3 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_D_Hertz_viscous_Coulomb_CL.h +++ b/applications/DEMApplication/custom_constitutive/DEM_D_Hertz_viscous_Coulomb_CL.h @@ -129,6 +129,7 @@ namespace Kratos { void CalculateInelasticViscodampingEnergyFEM(double& inelastic_viscodamping_energy, double ViscoDampingLocalContactForce[3], double LocalDeltDisp[3]); + double GetTangentialStiffness() override; private: diff --git a/applications/DEMApplication/custom_constitutive/DEM_D_Linear_classic_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_D_Linear_classic_CL.cpp index 9961888818f4..8629a018dc8e 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_D_Linear_classic_CL.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_D_Linear_classic_CL.cpp @@ -19,13 +19,15 @@ namespace Kratos{ return Kratos::make_unique(); } - void DEM_D_Linear_classic::InitializeContact(SphericParticle* const element1, SphericParticle* const element2, const double indentation) { + void DEM_D_Linear_classic::InitializeContact(SphericParticle* const element1, SphericParticle* const element2, const double indentation_particle) { //Get equivalent Radius const double my_radius = element1->GetRadius(); const double other_radius = element2->GetRadius(); const double radius_sum = my_radius + other_radius; const double min_radius = std::min(my_radius, other_radius); + //const double radius_sum_inv = 1.0 / radius_sum; + //const double equiv_radius = my_radius * other_radius * radius_sum_inv; //Get equivalent Young's Modulus const double my_young = element1->GetYoung(); @@ -40,9 +42,10 @@ namespace Kratos{ //const double other_shear_modulus = 0.5 * other_young / (1.0 + other_poisson); //const double equiv_shear = 1.0 / ((2.0 - my_poisson)/my_shear_modulus + (2.0 - other_poisson)/other_shear_modulus); - //Literature [Cundall, 2004, "A bonded particle model for rock"] [PFC 7.0 manual] + //Literature [Cundall, 2004, "A bonded particle model for rock"] [PFC 7.0 manual][mKt Mindlin model] mKn = equiv_young * Globals::Pi * min_radius * min_radius / radius_sum; mKt = mKn / (2.0 * (equiv_poisson + 1.0)); + //mKt = 8.0 * equiv_shear * sqrt(equiv_radius * indentation_particle); } void DEM_D_Linear_classic::InitializeContactWithFEM(SphericParticle* const element, Condition* const wall, const double indentation, const double ini_delta) { @@ -68,4 +71,8 @@ namespace Kratos{ mKn = equiv_young * Globals::Pi * effective_radius; mKt = mKn / (2.0 * (equiv_poisson + 1.0)); } + + double DEM_D_Linear_classic::GetTangentialStiffness() { + return mKt; + } } // namespace Kratos \ No newline at end of file diff --git a/applications/DEMApplication/custom_constitutive/DEM_D_Linear_classic_CL.h b/applications/DEMApplication/custom_constitutive/DEM_D_Linear_classic_CL.h index 0ea21af91bf6..65d1c0e6e77a 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_D_Linear_classic_CL.h +++ b/applications/DEMApplication/custom_constitutive/DEM_D_Linear_classic_CL.h @@ -33,6 +33,8 @@ namespace Kratos { void InitializeContactWithFEM(SphericParticle* const element, Condition* const wall, const double indentation, const double ini_delta = 0.0) override; + double GetTangentialStiffness() override; + protected: }; diff --git a/applications/DEMApplication/custom_constitutive/DEM_D_Linear_viscous_Coulomb_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_D_Linear_viscous_Coulomb_CL.cpp index 014254cbc5da..98dce56cb679 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_D_Linear_viscous_Coulomb_CL.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_D_Linear_viscous_Coulomb_CL.cpp @@ -170,6 +170,7 @@ namespace Kratos { ViscoDampingLocalContactForce[0] = - equiv_visco_damp_coeff_tangential * LocalRelVel[0]; ViscoDampingLocalContactForce[1] = - equiv_visco_damp_coeff_tangential * LocalRelVel[1]; ViscoDampingLocalContactForce[2] = - equiv_visco_damp_coeff_normal * LocalRelVel[2]; + } ///////////////////////// diff --git a/applications/DEMApplication/custom_constitutive/DEM_D_Quadratic_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_D_Quadratic_CL.cpp index 73fd90086155..5f0df4f827e2 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_D_Quadratic_CL.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_D_Quadratic_CL.cpp @@ -98,4 +98,8 @@ namespace Kratos { mKt = mKn / (2.0 * (equiv_poisson + 1.0)); } + double DEM_D_Quadratic::GetTangentialStiffness() { + return mKt; + } + } //namespace Kratos diff --git a/applications/DEMApplication/custom_constitutive/DEM_D_Quadratic_CL.h b/applications/DEMApplication/custom_constitutive/DEM_D_Quadratic_CL.h index 69902b8c5e43..cdd0497169cb 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_D_Quadratic_CL.h +++ b/applications/DEMApplication/custom_constitutive/DEM_D_Quadratic_CL.h @@ -45,6 +45,8 @@ namespace Kratos{ void InitializeContactWithFEM(SphericParticle* const element, Condition* const wall, const double indentation, const double ini_delta = 0.0) override; + double GetTangentialStiffness() override; + private: }; //CLASS DEM_D_QUADRATIC diff --git a/applications/DEMApplication/custom_constitutive/DEM_Dempack_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_Dempack_CL.cpp index 94df780ed168..0f83ec39d884 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_Dempack_CL.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_Dempack_CL.cpp @@ -146,6 +146,7 @@ namespace Kratos { double equiv_young, double equiv_shear, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -164,6 +165,7 @@ namespace Kratos { kn_el, equiv_young, indentation, + indentation_particle, calculation_area, acumulated_damage, element1, @@ -214,6 +216,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -505,7 +508,8 @@ namespace Kratos { double ElasticLocalRotationalMoment[3], double ViscoLocalRotationalMoment[3], double equiv_poisson, - double indentation, + double indentation, + double indentation_particle, double LocalElasticContactForce[3], double normalLocalContactForce, double GlobalElasticContactForces[3], diff --git a/applications/DEMApplication/custom_constitutive/DEM_Dempack_CL.h b/applications/DEMApplication/custom_constitutive/DEM_Dempack_CL.h index a5e19cddbd2c..fe7436b93169 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_Dempack_CL.h +++ b/applications/DEMApplication/custom_constitutive/DEM_Dempack_CL.h @@ -80,6 +80,7 @@ namespace Kratos { double equiv_young, double equiv_shear, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -96,6 +97,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -134,7 +136,8 @@ namespace Kratos { double ElasticLocalRotationalMoment[3], double ViscoLocalRotationalMoment[3], double equiv_poisson, - double indentation, + double indentation, + double indentation_particle, double LocalElasticContactForce[3], double normalLocalContactForce, double GlobalElasticContactForces[3], diff --git a/applications/DEMApplication/custom_constitutive/DEM_Dempack_dev_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_Dempack_dev_CL.cpp index 7f917a3d0c0f..e1151e6dba4b 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_Dempack_dev_CL.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_Dempack_dev_CL.cpp @@ -77,6 +77,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, diff --git a/applications/DEMApplication/custom_constitutive/DEM_Dempack_dev_CL.h b/applications/DEMApplication/custom_constitutive/DEM_Dempack_dev_CL.h index 0858ab52ed89..89f1ff68a197 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_Dempack_dev_CL.h +++ b/applications/DEMApplication/custom_constitutive/DEM_Dempack_dev_CL.h @@ -69,6 +69,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, diff --git a/applications/DEMApplication/custom_constitutive/DEM_ExponentialHC_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_ExponentialHC_CL.cpp index 0dbe265f4865..9113707b6577 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_ExponentialHC_CL.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_ExponentialHC_CL.cpp @@ -35,6 +35,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, diff --git a/applications/DEMApplication/custom_constitutive/DEM_ExponentialHC_CL.h b/applications/DEMApplication/custom_constitutive/DEM_ExponentialHC_CL.h index 684a4de3e46c..d66a0667a1f0 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_ExponentialHC_CL.h +++ b/applications/DEMApplication/custom_constitutive/DEM_ExponentialHC_CL.h @@ -37,6 +37,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, diff --git a/applications/DEMApplication/custom_constitutive/DEM_KDEM_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_KDEM_CL.cpp index 5fb4a0e5693a..c3387ceb62d3 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_KDEM_CL.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_KDEM_CL.cpp @@ -237,6 +237,7 @@ namespace Kratos { double equiv_young, double equiv_shear, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -254,6 +255,7 @@ namespace Kratos { kn_el, equiv_young, indentation, + indentation_particle, calculation_area, acumulated_damage, element1, @@ -304,6 +306,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -494,7 +497,8 @@ namespace Kratos { double ElasticLocalRotationalMoment[3], double ViscoLocalRotationalMoment[3], double equiv_poisson, - double indentation, + double indentation, + double indentation_particle, double LocalElasticContactForce[3], double normalLocalContactForce, double GlobalElasticContactForces[3], diff --git a/applications/DEMApplication/custom_constitutive/DEM_KDEM_CL.h b/applications/DEMApplication/custom_constitutive/DEM_KDEM_CL.h index 39472f5db819..ba772c25bdb7 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_KDEM_CL.h +++ b/applications/DEMApplication/custom_constitutive/DEM_KDEM_CL.h @@ -56,6 +56,7 @@ namespace Kratos { double equiv_young, double equiv_shear, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -72,6 +73,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -130,7 +132,8 @@ namespace Kratos { double ElasticLocalRotationalMoment[3], double ViscoLocalRotationalMoment[3], double equiv_poisson, - double indentation, + double indentation, + double indentation_particle, double LocalElasticContactForce[3], double normalLocalContactForce, double GlobalElasticContactForces[3], diff --git a/applications/DEMApplication/custom_constitutive/DEM_KDEM_Rankine_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_KDEM_Rankine_CL.cpp index 80a6e7ccf2b1..3b07ab61dbcb 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_KDEM_Rankine_CL.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_KDEM_Rankine_CL.cpp @@ -63,6 +63,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, diff --git a/applications/DEMApplication/custom_constitutive/DEM_KDEM_Rankine_CL.h b/applications/DEMApplication/custom_constitutive/DEM_KDEM_Rankine_CL.h index a88e2599c3ff..2480ebc0822c 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_KDEM_Rankine_CL.h +++ b/applications/DEMApplication/custom_constitutive/DEM_KDEM_Rankine_CL.h @@ -39,6 +39,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, diff --git a/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_CL.cpp index f246ebacfb06..320af3610344 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_CL.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_CL.cpp @@ -38,6 +38,7 @@ namespace Kratos { double equiv_young, double equiv_shear, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -56,6 +57,7 @@ namespace Kratos { kn_el, equiv_young, indentation, + indentation_particle, calculation_area, acumulated_damage, element1, @@ -108,6 +110,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, diff --git a/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_CL.h b/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_CL.h index 56e303de22e6..15802c40d4c3 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_CL.h +++ b/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_CL.h @@ -43,6 +43,7 @@ namespace Kratos { double equiv_young, double equiv_shear, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -59,6 +60,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, diff --git a/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_parallel_bond_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_parallel_bond_CL.cpp index 75ff53c4e270..090330b99961 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_parallel_bond_CL.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_parallel_bond_CL.cpp @@ -115,6 +115,7 @@ namespace Kratos { double equiv_young, double equiv_shear, double indentation, + double indentation_particle, double calculation_area, double& accumulated_damage, SphericContinuumParticle* element1, @@ -133,6 +134,7 @@ namespace Kratos { kn_el, equiv_young, indentation, + indentation_particle, calculation_area, accumulated_damage, element1, @@ -285,6 +287,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& accumulated_damage, SphericContinuumParticle* element1, diff --git a/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_parallel_bond_CL.h b/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_parallel_bond_CL.h index d9d4a7523c7d..de8ca73dfdc7 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_parallel_bond_CL.h +++ b/applications/DEMApplication/custom_constitutive/DEM_KDEM_with_damage_parallel_bond_CL.h @@ -47,6 +47,7 @@ namespace Kratos { double equiv_young, double equiv_shear, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -78,6 +79,7 @@ namespace Kratos { const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, diff --git a/applications/DEMApplication/custom_constitutive/DEM_compound_constitutive_law_for_PBM.h b/applications/DEMApplication/custom_constitutive/DEM_compound_constitutive_law_for_PBM.h index a575af42edca..4c7113dee002 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_compound_constitutive_law_for_PBM.h +++ b/applications/DEMApplication/custom_constitutive/DEM_compound_constitutive_law_for_PBM.h @@ -41,6 +41,10 @@ class KRATOS_API(DEM_APPLICATION) DEM_compound_constitutive_law_for_PBM : public return mCCL.CalculateViscoDampingForce(LocalRelVel, UnbondedViscoDampingLocalContactForce, element1, element2); } + double GetTangentialStiffness() override { + return mCCL.GetTangentialStiffness(); + } + private: UnbondCL mCCL; diff --git a/applications/DEMApplication/custom_constitutive/DEM_continuum_constitutive_law.cpp b/applications/DEMApplication/custom_constitutive/DEM_continuum_constitutive_law.cpp index e6115faac95b..780368e2c5bf 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_continuum_constitutive_law.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_continuum_constitutive_law.cpp @@ -93,4 +93,9 @@ namespace Kratos { return false; } + double DEMContinuumConstitutiveLaw::GetTangentialStiffness() { + KRATOS_ERROR << "This function (DEMContinuumConstitutiveLaw::LocalMaxSearchDistance) shouldn't be accessed, use derived class instead"<Has(DEBUG_PRINTING_OPTION)) { + mDebugPrintingOption = false; + } else { + mDebugPrintingOption = bool((*mpProperties)[DEBUG_PRINTING_OPTION]); + } + if (mDebugPrintingOption) { + if (!mpProperties->Has(DEBUG_PRINTING_ID_1) || !mpProperties->Has(DEBUG_PRINTING_ID_2)) { + KRATOS_WARNING("DEM") << "\nWARNING: We are currently in DEBUG PRINTING mode, so the ids of the two particles involved must be given.\n\n"; + } + } + + #pragma omp critical + { + srand(element1->GetId()); + mBondSigmaMax = rand_normal((*mpProperties)[BOND_SIGMA_MAX], (*mpProperties)[BOND_SIGMA_MAX_DEVIATION]); + srand(element1->GetId()); + mBondTauZero = rand_normal((*mpProperties)[BOND_TAU_ZERO], (*mpProperties)[BOND_TAU_ZERO_DEVIATION]); + } + +} + +double DEM_parallel_bond::rand_normal(const double mean, const double stddev) { + + KRATOS_TRY + + if (!stddev) return mean; + double x, y, r; + + do { + x = 2.0 * rand() / RAND_MAX - 1; + y = 2.0 * rand() / RAND_MAX - 1; + r = x*x + y*y; + } while (r == 0.0 || r > 1.0); + + double d = sqrt(- 2.0 * log(r) / r); + return x * d * stddev + mean; + + KRATOS_CATCH("") +} + void DEM_parallel_bond::TransferParametersToProperties(const Parameters& parameters, Properties::Pointer pProp){ BaseClassType::TransferParametersToProperties(parameters, pProp); } @@ -203,13 +248,17 @@ void DEM_parallel_bond::GetContactArea(const double radius, const double other_r } // Here we calculate the mKn and mKt -void DEM_parallel_bond::InitializeContact(SphericParticle* const element1, SphericParticle* const element2, const double indentation) { +void DEM_parallel_bond::InitializeContact(SphericParticle* const element1, SphericParticle* const element2, const double indentation_particle) { + KRATOS_ERROR << "This function shouldn't be accessed here, use basic contact model instead."<mIniNeighbourFailureId[i_neighbour_count]; const double bonded_indentation = indentation - mInitialIndentationForBondedPart; - double BondedLocalElasticContactForce2 = 0.0; - if (!failure_type){ //if the bond is not broken - BondedLocalElasticContactForce2 = kn_el * bonded_indentation; + mBondedLocalElasticContactForce2 = kn_el * bonded_indentation; } else { //else the bond is broken - //if the bond is broken, we still calculate the normal compressive force but not the normal tensile force + //if the bond is broken, we still calculate the normal compressive force but not the normal tensile force //TODO:If the bond disappears, this will not work if (bonded_indentation > 0.0){ - BondedLocalElasticContactForce2 = kn_el * bonded_indentation; + mBondedLocalElasticContactForce2 = kn_el * bonded_indentation; } else { - BondedLocalElasticContactForce2 = 0.0; + mBondedLocalElasticContactForce2 = 0.0; } + //mBondedLocalElasticContactForce2 = 0.0; } - if (indentation > 0.0) { - mUnbondedLocalElasticContactForce2 = ComputeNormalUnbondedForce(indentation); + if (indentation_particle > 0.0) { + mUnbondedLocalElasticContactForce2 = ComputeNormalUnbondedForce(indentation_particle); } else { mUnbondedLocalElasticContactForce2 = 0.0; } - if(calculation_area){ - contact_sigma = BondedLocalElasticContactForce2 / calculation_area; - } + //if(calculation_area){ + // contact_sigma = mBondedLocalElasticContactForce2 / calculation_area; + //} - LocalElasticContactForce[2] = BondedLocalElasticContactForce2 + mUnbondedLocalElasticContactForce2; + LocalElasticContactForce[2] = mBondedLocalElasticContactForce2 + mUnbondedLocalElasticContactForce2; if (LocalElasticContactForce[2]) { - mBondedScalingFactor[2] = BondedLocalElasticContactForce2 / LocalElasticContactForce[2]; + mBondedScalingFactor[2] = mBondedLocalElasticContactForce2 / LocalElasticContactForce[2]; } else { mBondedScalingFactor[2] = 0.0; } + //for debug + if (mDebugPrintingOption) { + + const long unsigned int& sphere_id = (*mpProperties)[DEBUG_PRINTING_ID_1]; + const long unsigned int& neigh_sphere_id = (*mpProperties)[DEBUG_PRINTING_ID_2]; + + if ((element1->Id() == sphere_id) && (element2->Id() == neigh_sphere_id)) { + std::ofstream normal_forces_file("delta_stress_normal.txt", std::ios_base::out | std::ios_base::app); + normal_forces_file << r_process_info[TIME] << " " << bonded_indentation/*2*/ << " " << contact_sigma/*3*/ << " " << indentation <<'\n'; + normal_forces_file.flush(); + normal_forces_file.close(); + } + } + KRATOS_CATCH("") } // CalculateNormalForces @@ -446,7 +514,7 @@ void DEM_parallel_bond::CalculateViscoDampingCoeff(double &equiv_visco_damp_coef equiv_visco_damp_coeff_normal = 2.0 * equiv_gamma * sqrt(equiv_mass * kn_el); equiv_visco_damp_coeff_tangential = 2.0 * equiv_gamma * sqrt(equiv_mass * kt_el); - KRATOS_CATCH("") + KRATOS_CATCH("") } void DEM_parallel_bond::CalculateUnbondedViscoDampingForce(double LocalRelVel[3], @@ -460,7 +528,7 @@ void DEM_parallel_bond::CalculateUnbondedViscoDampingForce(double LocalRelVel[3] void DEM_parallel_bond::CalculateViscoDamping(double LocalRelVel[3], double ViscoDampingLocalContactForce[3], - double indentation, + double indentation_particle, double equiv_visco_damp_coeff_normal, double equiv_visco_damp_coeff_tangential, bool& sliding, @@ -477,7 +545,7 @@ void DEM_parallel_bond::CalculateViscoDamping(double LocalRelVel[3], mBondedViscoDampingLocalContactForce[1] = 0.0; mBondedViscoDampingLocalContactForce[2] = 0.0; - if (indentation > 0) { + if (indentation_particle > 0) { CalculateUnbondedViscoDampingForce(LocalRelVel, mUnbondedViscoDampingLocalContactForce, element1, element2); } @@ -498,6 +566,8 @@ void DEM_parallel_bond::CalculateViscoDamping(double LocalRelVel[3], ViscoDampingLocalContactForce[2] = mUnbondedViscoDampingLocalContactForce[2] + mBondedViscoDampingLocalContactForce[2]; } + mBondedLocalContactForce[2] = mBondedLocalElasticContactForce2 + mBondedViscoDampingLocalContactForce[2]; + KRATOS_CATCH("") } @@ -510,8 +580,9 @@ void DEM_parallel_bond::CalculateTangentialForces(double OldLocalElasticContactF double LocalRelVel[3], const double kt_el, const double equiv_shear, + double& contact_sigma, double& contact_tau, - double indentation, + double indentation_particle, double calculation_area, double& failure_criterion_state, SphericContinuumParticle* element1, @@ -529,27 +600,34 @@ void DEM_parallel_bond::CalculateTangentialForces(double OldLocalElasticContactF double current_tangential_force_module = 0.0; // The [mBondedScalingFactor] is divided into two direction [0] and [1]. June, 2022 + double OldBondedLocalElasticContactForce[2] = {0.0}; + OldBondedLocalElasticContactForce[0] = mBondedScalingFactor[0] * OldLocalElasticContactForce[0]; + OldBondedLocalElasticContactForce[1] = mBondedScalingFactor[1] * OldLocalElasticContactForce[1]; + double BondedLocalElasticContactForce[2] = {0.0}; // bond force if (!failure_type) { - mAccumulatedBondedTangentialLocalDisplacement[0] += LocalDeltDisp[0]; - mAccumulatedBondedTangentialLocalDisplacement[1] += LocalDeltDisp[1]; - BondedLocalElasticContactForce[0] -= kt_el * mAccumulatedBondedTangentialLocalDisplacement[0]; // 0: first tangential - BondedLocalElasticContactForce[1] -= kt_el * mAccumulatedBondedTangentialLocalDisplacement[1]; // 1: second tangential + //mAccumulatedBondedTangentialLocalDisplacement[0] += LocalDeltDisp[0]; + //mAccumulatedBondedTangentialLocalDisplacement[1] += LocalDeltDisp[1]; + //BondedLocalElasticContactForce[0] -= kt_el * mAccumulatedBondedTangentialLocalDisplacement[0]; // 0: first tangential + //BondedLocalElasticContactForce[1] -= kt_el * mAccumulatedBondedTangentialLocalDisplacement[1]; // 1: second tangential + BondedLocalElasticContactForce[0] = OldBondedLocalElasticContactForce[0] - kt_el * LocalDeltDisp[0]; // 0: first tangential + BondedLocalElasticContactForce[1] = OldBondedLocalElasticContactForce[1] - kt_el * LocalDeltDisp[1]; // 1: second tangential } else { //TODO: maybe a friction force due to the broekn bond should be added here BondedLocalElasticContactForce[0] = 0.0; // 0: first tangential BondedLocalElasticContactForce[1] = 0.0; // 1: second tangential } - current_tangential_force_module = sqrt(BondedLocalElasticContactForce[0] * BondedLocalElasticContactForce[0] - + BondedLocalElasticContactForce[1] * BondedLocalElasticContactForce[1]); + //current_tangential_force_module = sqrt(BondedLocalElasticContactForce[0] * BondedLocalElasticContactForce[0] + // + BondedLocalElasticContactForce[1] * BondedLocalElasticContactForce[1]); - if (calculation_area){ - contact_tau = current_tangential_force_module / calculation_area; - } + //if (calculation_area){ + // contact_tau = current_tangential_force_module / calculation_area; + // contact_sigma = mBondedLocalContactForce / calculation_area; + //} //unbonded force double OldUnbondedLocalElasticContactForce[2] = {0.0}; @@ -558,9 +636,11 @@ void DEM_parallel_bond::CalculateTangentialForces(double OldLocalElasticContactF double max_admissible_shear_force = 0.0; double fraction = 0.0; - if (indentation <= 0.0) { + if (indentation_particle <= 0.0) { UnbondedLocalElasticContactForce[0] = 0.0; UnbondedLocalElasticContactForce[1] = 0.0; + //ViscoDampingLocalContactForce[0] = 0.0; + //ViscoDampingLocalContactForce[1] = 0.0; } else { // Here, unBondedScalingFactor[] = 1 - mBondedScalingFactor[] OldUnbondedLocalElasticContactForce[0] = (1 - mBondedScalingFactor[0]) * OldLocalElasticContactForce[0]; @@ -636,6 +716,7 @@ void DEM_parallel_bond::CalculateTangentialForces(double OldLocalElasticContactF mUnbondedViscoDampingLocalContactForce[1] = 0.0; } } + ViscoDampingLocalContactForce[0] = mBondedViscoDampingLocalContactForce[0] + mUnbondedViscoDampingLocalContactForce[0]; ViscoDampingLocalContactForce[1] = mBondedViscoDampingLocalContactForce[1] + mUnbondedViscoDampingLocalContactForce[1]; sliding = true; @@ -646,11 +727,41 @@ void DEM_parallel_bond::CalculateTangentialForces(double OldLocalElasticContactF LocalElasticContactForce[1] = BondedLocalElasticContactForce[1] + UnbondedLocalElasticContactForce[1]; // Here, we only calculate the BondedScalingFactor and [unBondedScalingFactor = 1 - BondedScalingFactor]. - if (LocalElasticContactForce[0] && LocalElasticContactForce[1]) { + if (LocalElasticContactForce[0]) { mBondedScalingFactor[0] = BondedLocalElasticContactForce[0] / LocalElasticContactForce[0]; + } else { + mBondedScalingFactor[0] = 0.0; + } + + if (LocalElasticContactForce[1]) { mBondedScalingFactor[1] = BondedLocalElasticContactForce[1] / LocalElasticContactForce[1]; } else { - mBondedScalingFactor[0] = mBondedScalingFactor[1] = 0.0; + mBondedScalingFactor[1] = 0.0; + } + + current_tangential_force_module = sqrt((BondedLocalElasticContactForce[0] + mBondedViscoDampingLocalContactForce[0]) + * (BondedLocalElasticContactForce[0] + mBondedViscoDampingLocalContactForce[0]) + + (BondedLocalElasticContactForce[1] + mBondedViscoDampingLocalContactForce[1]) + * (BondedLocalElasticContactForce[1] + mBondedViscoDampingLocalContactForce[1])); + + if (calculation_area){ + contact_tau = current_tangential_force_module / calculation_area; + contact_sigma = mBondedLocalContactForce[2] / calculation_area; + } + + //for debug + if (mDebugPrintingOption) { + + const long unsigned int& sphere_id = (*mpProperties)[DEBUG_PRINTING_ID_1]; + const long unsigned int& neigh_sphere_id = (*mpProperties)[DEBUG_PRINTING_ID_2]; + const double AccumulatedBondedTangentialLocalDisplacementModulus = sqrt(mAccumulatedBondedTangentialLocalDisplacement[0]*mAccumulatedBondedTangentialLocalDisplacement[0] + mAccumulatedBondedTangentialLocalDisplacement[1]*mAccumulatedBondedTangentialLocalDisplacement[1]); + + if ((element1->Id() == sphere_id) && (element2->Id() == neigh_sphere_id)) { + std::ofstream normal_forces_file("delta_stress_tangential.txt", std::ios_base::out | std::ios_base::app); + normal_forces_file << r_process_info[TIME] << " " << AccumulatedBondedTangentialLocalDisplacementModulus/*4*/ << " " << contact_tau/*5*/ << '\n'; + normal_forces_file.flush(); + normal_forces_file.close(); + } } KRATOS_CATCH("") @@ -669,10 +780,11 @@ void DEM_parallel_bond::CalculateMoments(SphericContinuumParticle* element, double ElasticLocalRotationalMoment[3], double ViscoLocalRotationalMoment[3], double equiv_poisson, - double indentation, - double LocalElasticContactForce[3], + double indentation, + double indentation_particle, + double LocalContactForce[3], double normalLocalContactForce, - double GlobalElasticContactForces[3], + double GlobalContactForce[3], double LocalCoordSystem_2[3], const int i_neighbor_count) { @@ -681,32 +793,26 @@ void DEM_parallel_bond::CalculateMoments(SphericContinuumParticle* element, int failure_type = element->mIniNeighbourFailureId[i_neighbor_count]; if (failure_type == 0) { - ComputeParticleRotationalMoments(element, - neighbor, - equiv_young, - distance, - calculation_area, - LocalCoordSystem, - ElasticLocalRotationalMoment, - ViscoLocalRotationalMoment, - equiv_poisson, - indentation, - LocalElasticContactForce); - } - - double LocalUnbondElasticContactForce[3] = {0.0}; - double GlobalUnbondElasticContactForce[3] = {0.0}; - LocalUnbondElasticContactForce[0] = (1 - mBondedScalingFactor[0]) * LocalElasticContactForce[0]; - LocalUnbondElasticContactForce[1] = (1 - mBondedScalingFactor[1]) * LocalElasticContactForce[1]; - LocalUnbondElasticContactForce[2] = (1 - mBondedScalingFactor[2]) * LocalElasticContactForce[2]; - GeometryFunctions::VectorLocal2Global(LocalCoordSystem, LocalUnbondElasticContactForce, GlobalUnbondElasticContactForce); + ComputeParticleRotationalMoments(element, + neighbor, + equiv_young, + distance, + calculation_area, + LocalCoordSystem, + ElasticLocalRotationalMoment, + ViscoLocalRotationalMoment, + equiv_poisson, + indentation, + LocalContactForce); + CalculateBondRotationalDamping(element, neighbor, LocalCoordSystem, ViscoLocalRotationalMoment); + } DemContact::ComputeParticleContactMoments(normalLocalContactForce, - GlobalUnbondElasticContactForce, + GlobalContactForce, LocalCoordSystem_2, element, neighbor, - indentation, + indentation_particle, i_neighbor_count); KRATOS_CATCH("") @@ -717,6 +823,14 @@ double DEM_parallel_bond::GetYoungModulusForComputingRotationalMoments(const dou return bond_equiv_young; } +double DEM_parallel_bond::GetBondKn(double bond_equiv_young, double calculation_area, double distance){ + KRATOS_TRY + + return (bond_equiv_young * calculation_area / distance); + + KRATOS_CATCH("") +} + void DEM_parallel_bond::ComputeParticleRotationalMoments(SphericContinuumParticle* element, SphericContinuumParticle* neighbor, double equiv_young, @@ -727,10 +841,11 @@ void DEM_parallel_bond::ComputeParticleRotationalMoments(SphericContinuumParticl double ViscoLocalRotationalMoment[3], double equiv_poisson, double indentation, - double LocalElasticContactForce[3]) { + double LocalContactForce[3]) { KRATOS_TRY - + //const double& bond_rotational_moment_coefficient_normal =(*mpProperties)[BOND_ROTATIONAL_MOMENT_COEFFICIENT_NORMAL]; + //const double& bond_rotational_moment_coefficient_tangential =(*mpProperties)[BOND_ROTATIONAL_MOMENT_COEFFICIENT_TANGENTIAL]; double LocalDeltaRotatedAngle[3] = {0.0}; double LocalDeltaAngularVelocity[3] = {0.0}; @@ -745,29 +860,28 @@ void DEM_parallel_bond::ComputeParticleRotationalMoments(SphericContinuumParticl GeometryFunctions::VectorGlobal2Local(LocalCoordSystem, GlobalDeltaAngularVelocity, LocalDeltaAngularVelocity); const double equivalent_radius = std::sqrt(calculation_area / Globals::Pi); - const double element_mass = element->GetMass(); - const double neighbor_mass = neighbor->GetMass(); - const double equiv_mass = element_mass * neighbor_mass / (element_mass + neighbor_mass); + //const double element_mass = element->GetMass(); + //const double neighbor_mass = neighbor->GetMass(); + //const double equiv_mass = element_mass * neighbor_mass / (element_mass + neighbor_mass); const double bond_equiv_young = GetYoungModulusForComputingRotationalMoments(equiv_young); - double kn_el = bond_equiv_young * calculation_area / distance; + //double kn_el = bond_equiv_young * calculation_area / distance; + double kn_el = GetBondKn(bond_equiv_young, calculation_area, distance); double kt_el = kn_el / (*mpProperties)[BOND_KNKS_RATIO]; const double Inertia_I = 0.25 * Globals::Pi * equivalent_radius * equivalent_radius * equivalent_radius * equivalent_radius; const double Inertia_J = 2.0 * Inertia_I; // This is the polar inertia - const double& damping_gamma = (*mpProperties)[DAMPING_GAMMA]; + //const double& damping_gamma = (*mpProperties)[DAMPING_GAMMA]; //Viscous parameter taken from Olmedo et al., 'Discrete element model of the dynamic response of fresh wood stems to impact' - array_1d visc_param; + /*array_1d visc_param; visc_param[0] = 2.0 * damping_gamma * std::sqrt(equiv_mass * bond_equiv_young * Inertia_I / distance); // OLMEDO visc_param[1] = 2.0 * damping_gamma * std::sqrt(equiv_mass * bond_equiv_young * Inertia_I / distance); // OLMEDO - visc_param[2] = 2.0 * damping_gamma * std::sqrt(equiv_mass * bond_equiv_young * Inertia_J / distance); // OLMEDO + visc_param[2] = 2.0 * damping_gamma * std::sqrt(equiv_mass * bond_equiv_young * Inertia_J / distance); // OLMEDO*/ - double aux = 0.0; - aux = (element->GetRadius() + neighbor->GetRadius()) / distance; // This is necessary because if spheres are not tangent the DeltaAngularVelocity has to be interpolated - + double aux = (element->GetRadius() + neighbor->GetRadius()) / distance; // This is necessary because if spheres are not tangent the DeltaAngularVelocity has to be interpolated array_1d LocalEffDeltaRotatedAngle; LocalEffDeltaRotatedAngle[0] = LocalDeltaRotatedAngle[0] * aux; @@ -783,6 +897,22 @@ void DEM_parallel_bond::ComputeParticleRotationalMoments(SphericContinuumParticl ElasticLocalRotationalMoment[1] = -kn_el / calculation_area * Inertia_I * LocalEffDeltaRotatedAngle[1]; ElasticLocalRotationalMoment[2] = -kt_el / calculation_area * Inertia_J * LocalEffDeltaRotatedAngle[2]; + /*ViscoLocalRotationalMoment[0] = -visc_param[0] * LocalEffDeltaAngularVelocity[0]; + ViscoLocalRotationalMoment[1] = -visc_param[1] * LocalEffDeltaAngularVelocity[1]; + ViscoLocalRotationalMoment[2] = -visc_param[2] * LocalEffDeltaAngularVelocity[2];*/ + + //DEM_MULTIPLY_BY_SCALAR_3(ElasticLocalRotationalMoment, rotational_moment_coeff); + //DEM_MULTIPLY_BY_SCALAR_3(ViscoLocalRotationalMoment, rotational_moment_coeff); + + KRATOS_CATCH("") +}//ComputeParticleRotationalMoments + +void DEM_parallel_bond::CalculateBondRotationalDamping(SphericContinuumParticle* element, + SphericContinuumParticle* neighbor, + double LocalCoordSystem[3][3], + double ViscoLocalRotationalMoment[3]){ + KRATOS_TRY + // Bond rotational 'friction' based on particle rolling fricton //Not damping but simple implementation to help energy dissipation @@ -805,30 +935,24 @@ void DEM_parallel_bond::ComputeParticleRotationalMoments(SphericContinuumParticl element1AngularVelocity_normalise[2] = LocalElement1AngularVelocity[2] / element1AngularVelocity_modulus; Properties& properties_of_this_contact = element->GetProperties().GetSubProperties(neighbor->GetProperties().Id()); - - double BondedLocalElasticContactForce[3] = {0.0}; - BondedLocalElasticContactForce[0] = mBondedScalingFactor[0] * LocalElasticContactForce[0]; - BondedLocalElasticContactForce[1] = mBondedScalingFactor[1] * LocalElasticContactForce[1]; - BondedLocalElasticContactForce[2] = mBondedScalingFactor[2] * LocalElasticContactForce[2]; - ViscoLocalRotationalMoment[0] = - element1AngularVelocity_normalise[0] * std::abs(BondedLocalElasticContactForce[2]) * bond_center_point_to_element1_mass_center_distance + ViscoLocalRotationalMoment[0] = - element1AngularVelocity_normalise[0] * std::abs(mBondedLocalContactForce[2]) * bond_center_point_to_element1_mass_center_distance * properties_of_this_contact[ROLLING_FRICTION]; - ViscoLocalRotationalMoment[1] = - element1AngularVelocity_normalise[1] * std::abs(BondedLocalElasticContactForce[2]) * bond_center_point_to_element1_mass_center_distance + ViscoLocalRotationalMoment[1] = - element1AngularVelocity_normalise[1] * std::abs(mBondedLocalContactForce[2]) * bond_center_point_to_element1_mass_center_distance * properties_of_this_contact[ROLLING_FRICTION]; - ViscoLocalRotationalMoment[2] = - element1AngularVelocity_normalise[2] * std::abs(BondedLocalElasticContactForce[2]) * bond_center_point_to_element1_mass_center_distance + ViscoLocalRotationalMoment[2] = - element1AngularVelocity_normalise[2] * std::abs(mBondedLocalContactForce[2]) * bond_center_point_to_element1_mass_center_distance * properties_of_this_contact[ROLLING_FRICTION]; } else { ViscoLocalRotationalMoment[0] = 0.0; ViscoLocalRotationalMoment[1] = 0.0; ViscoLocalRotationalMoment[2] = 0.0; - } - - KRATOS_CATCH("") -}//ComputeParticleRotationalMoments + } + KRATOS_CATCH("") +} //************************************* // Bond failure checking @@ -851,18 +975,18 @@ void DEM_parallel_bond::CheckFailure(const int i_neighbour_count, if (failure_type == 0) { //parameters - const double& bond_sigma_max = (*mpProperties)[BOND_SIGMA_MAX]; + //const double& bond_sigma_max = (*mpProperties)[BOND_SIGMA_MAX]; //const double& bond_sigma_max_deviation = (*mpProperties)[BOND_SIGMA_MAX_DEVIATION]; - const double& bond_tau_zero = (*mpProperties)[BOND_TAU_ZERO]; + //const double& bond_tau_zero = (*mpProperties)[BOND_TAU_ZERO]; //const double& bond_tau_zero_deviation = (*mpProperties)[BOND_TAU_ZERO_DEVIATION]; const double& bond_interanl_friction = (*mpProperties)[BOND_INTERNAL_FRICC]; const double& bond_rotational_moment_coefficient_normal =(*mpProperties)[BOND_ROTATIONAL_MOMENT_COEFFICIENT_NORMAL]; const double& bond_rotational_moment_coefficient_tangential =(*mpProperties)[BOND_ROTATIONAL_MOMENT_COEFFICIENT_TANGENTIAL]; double bond_rotational_moment[3] = {0.0}; - bond_rotational_moment[0] = ElasticLocalRotationalMoment[0]; - bond_rotational_moment[1] = ElasticLocalRotationalMoment[1]; - bond_rotational_moment[2] = ElasticLocalRotationalMoment[2]; + bond_rotational_moment[0] = ElasticLocalRotationalMoment[0] + ViscoLocalRotationalMoment[0]; + bond_rotational_moment[1] = ElasticLocalRotationalMoment[1] + ViscoLocalRotationalMoment[1]; + bond_rotational_moment[2] = ElasticLocalRotationalMoment[2] + ViscoLocalRotationalMoment[2]; double bond_rotational_moment_normal_modulus = 0.0; double bond_rotational_moment_tangential_modulus = 0.0; @@ -878,14 +1002,15 @@ void DEM_parallel_bond::CheckFailure(const int i_neighbour_count, const double I = 0.25 * Globals::Pi * bond_radius * bond_radius * bond_radius * bond_radius; const double J = 2.0 * I; // This is the polar inertia - double bond_current_tau_max = bond_tau_zero; + double bond_current_tau_max = mBondTauZero; if (contact_sigma >= 0) { bond_current_tau_max += tan(bond_interanl_friction * Globals::Pi / 180.0) * contact_sigma; } + //bond_current_tau_max += tan(bond_interanl_friction * Globals::Pi / 180.0) * contact_sigma; - if (contact_sigma < 0.0 /*break only in tension*/ - && (-1 * contact_sigma + bond_rotational_moment_coefficient_normal * bond_rotational_moment_tangential_modulus * bond_radius / I > bond_sigma_max) + if (contact_sigma < 0.0 /*break in tension*/ + && ((-1 * contact_sigma + bond_rotational_moment_coefficient_normal * bond_rotational_moment_tangential_modulus * bond_radius / I) > mBondSigmaMax) && !(*mpProperties)[IS_UNBREAKABLE]) { //for normal failure_type = 4; // failure in tension @@ -900,12 +1025,14 @@ void DEM_parallel_bond::CheckFailure(const int i_neighbour_count, ElasticLocalRotationalMoment[0] = 0.0; ElasticLocalRotationalMoment[1] = 0.0; ElasticLocalRotationalMoment[2] = 0.0; - ViscoLocalRotationalMoment[0] = 0.0; - ViscoLocalRotationalMoment[1] = 0.0; - ViscoLocalRotationalMoment[2] = 0.0; - } - else if(( std::abs(contact_tau) + bond_rotational_moment_coefficient_tangential * bond_rotational_moment_normal_modulus * bond_radius / J > bond_current_tau_max) - && !(*mpProperties)[IS_UNBREAKABLE]) + //ViscoLocalRotationalMoment[0] = 0.0; + //ViscoLocalRotationalMoment[1] = 0.0; + //ViscoLocalRotationalMoment[2] = 0.0; + mBondedScalingFactor[0] = 0.0; + mBondedScalingFactor[1] = 0.0; + mBondedScalingFactor[2] = 0.0; + } else if(((std::abs(contact_tau) + bond_rotational_moment_coefficient_tangential * bond_rotational_moment_normal_modulus * bond_radius / J) > bond_current_tau_max) + && !(*mpProperties)[IS_UNBREAKABLE]) { //for tangential failure_type = 2; // failure in shear contact_sigma = 0.0; @@ -913,16 +1040,19 @@ void DEM_parallel_bond::CheckFailure(const int i_neighbour_count, //If bond break in shear, the normal compressive force should still be there like before LocalElasticContactForce[0] *= (1 - mBondedScalingFactor[0]); LocalElasticContactForce[1] *= (1 - mBondedScalingFactor[1]); - //LocalElasticContactForce[2] = mUnbondedLocalElasticContactForce2; + LocalElasticContactForce[2] = mUnbondedLocalElasticContactForce2; ViscoDampingLocalContactForce[0] = mUnbondedViscoDampingLocalContactForce[0]; ViscoDampingLocalContactForce[1] = mUnbondedViscoDampingLocalContactForce[1]; - //ViscoDampingLocalContactForce[2] = mUnbondedViscoDampingLocalContactForce[2]; + ViscoDampingLocalContactForce[2] = mUnbondedViscoDampingLocalContactForce[2]; ElasticLocalRotationalMoment[0] = 0.0; ElasticLocalRotationalMoment[1] = 0.0; ElasticLocalRotationalMoment[2] = 0.0; - ViscoLocalRotationalMoment[0] = 0.0; - ViscoLocalRotationalMoment[1] = 0.0; - ViscoLocalRotationalMoment[2] = 0.0; + //ViscoLocalRotationalMoment[0] = 0.0; + //ViscoLocalRotationalMoment[1] = 0.0; + //ViscoLocalRotationalMoment[2] = 0.0; + mBondedScalingFactor[0] = 0.0; + mBondedScalingFactor[1] = 0.0; + mBondedScalingFactor[2] = 0.0; } } diff --git a/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_CL.h b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_CL.h index c1808e855317..69c95e3917a2 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_CL.h +++ b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_CL.h @@ -27,6 +27,7 @@ namespace Kratos{ void TransferParametersToProperties(const Parameters& parameters, Properties::Pointer pProp) override; std::string GetTypeOfLaw() override; void Check(Properties::Pointer pProp) const override; + void Initialize(SphericContinuumParticle* element1, SphericContinuumParticle* element2, Properties::Pointer pProps) override; ~DEM_parallel_bond() {} @@ -37,7 +38,7 @@ namespace Kratos{ void GetContactArea(const double radius, const double other_radius, const Vector& vector_of_initial_areas, const int neighbour_position, double& calculation_area) override; void CalculateElasticConstants(double& kn_el, double& kt_el, double initial_dist, double equiv_young, double equiv_poisson, double calculation_area, SphericContinuumParticle* element1, SphericContinuumParticle* element2, double indentation) override; - virtual void InitializeContact(SphericParticle* const element1, SphericParticle* const element2, const double indentation); + virtual void InitializeContact(SphericParticle* const element1, SphericParticle* const element2, const double indentation_particle); virtual void CalculateUnbondedViscoDampingForce(double LocalRelVel[3], double UnbondedViscoDampingLocalContactForce[3], @@ -52,8 +53,9 @@ namespace Kratos{ //TODO:CHECK virtual double GetYoungModulusForComputingRotationalMoments(const double& equiv_young); + virtual double GetBondKn(double bond_equiv_young, double calculation_area, double distance); - virtual void CheckFailure(const int i_neighbour_count, + void CheckFailure(const int i_neighbour_count, SphericContinuumParticle* element1, SphericContinuumParticle* element2, double& contact_sigma, @@ -77,6 +79,7 @@ namespace Kratos{ double equiv_young, double equiv_shear, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -91,10 +94,11 @@ namespace Kratos{ virtual double ComputeNormalUnbondedForce(double unbonded_indentation); - void CalculateNormalForces(double LocalElasticContactForce[3], + virtual void CalculateNormalForces(double LocalElasticContactForce[3], const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -104,7 +108,7 @@ namespace Kratos{ const ProcessInfo& r_process_info, double& contact_sigma); - virtual void CalculateViscoDampingCoeff(double &equiv_visco_damp_coeff_normal, + void CalculateViscoDampingCoeff(double &equiv_visco_damp_coeff_normal, double &equiv_visco_damp_coeff_tangential, SphericContinuumParticle* element1, SphericContinuumParticle* element2, @@ -113,7 +117,7 @@ namespace Kratos{ void CalculateViscoDamping(double LocalRelVel[3], double ViscoDampingLocalContactForce[3], - double indentation, + double indentation_particle, double equiv_visco_damp_coeff_normal, double equiv_visco_damp_coeff_tangential, bool& sliding, @@ -131,15 +135,16 @@ namespace Kratos{ double LocalRelVel[3], const double kt_el, const double equiv_shear, + double& contact_sigma, double& contact_tau, - double indentation, + double indentation_particle, double calculation_area, double& failure_criterion_state, SphericContinuumParticle* element1, SphericContinuumParticle* element2, int i_neighbour_count, bool& sliding, - const ProcessInfo& r_process_info); + const ProcessInfo& r_process_info) override; void CalculateMoments(SphericContinuumParticle* element, SphericContinuumParticle* neighbor, @@ -151,6 +156,7 @@ namespace Kratos{ double ViscoLocalRotationalMoment[3], double equiv_poisson, double indentation, + double indentation_particle, double LocalElasticContactForce[3], double normalLocalContactForce, double GlobalElasticContactForces[3], @@ -167,7 +173,12 @@ namespace Kratos{ double ViscoLocalRotationalMoment[3], double equiv_poisson, double indentation, - double LocalElasticContactForce[3]) override; + double LocalElasticContactForce[3]) override; + + void CalculateBondRotationalDamping(SphericContinuumParticle* element, + SphericContinuumParticle* neighbor, + double LocalCoordSystem[3][3], + double ViscoLocalRotationalMoment[3]); void AddContributionOfShearStrainParallelToBond(double OldLocalElasticContactForce[3], double LocalElasticExtraContactForce[3], @@ -178,21 +189,30 @@ namespace Kratos{ SphericContinuumParticle* element1, SphericContinuumParticle* element2); + double GetTangentialStiffness() override; + + double rand_normal(const double mean, const double stddev); double mUnbondedLocalElasticContactForce2 = 0.0; double mUnbondedNormalElasticConstant = 0.0; double mUnbondedTangentialElasticConstant = 0.0; double mUnbondedViscoDampingLocalContactForce[3] = {0.0}; double mBondedViscoDampingLocalContactForce[3] = {0.0}; + double mBondedLocalContactForce[3] = {0.0}; + double mBondedLocalElasticContactForce2 = 0.0; double mBondedScalingFactor[3] = {0.0}; double mUnbondedEquivViscoDampCoeffTangential = 0.0; double mUnbondedEquivViscoDampCoeffNormal = 0.0; double mInitialIndentationForBondedPart = 0.0; double mAccumulatedBondedTangentialLocalDisplacement[2] = {0.0}; - double mBondedLocalContactNormalTorque[3] = {0.0}; - double mBondedLocalContactTangentTorque[3] = {0.0}; + //double mBondedLocalContactNormalTorque[3] = {0.0}; + //double mBondedLocalContactTangentTorque[3] = {0.0}; double mKn; double mKt; + bool mDebugPrintingOption = false; + double mBondSigmaMax = 0.0; + double mBondTauZero = 0.0; + protected: diff --git a/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_bilinear_damage_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_bilinear_damage_CL.cpp new file mode 100644 index 000000000000..f8b8a56b9dbf --- /dev/null +++ b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_bilinear_damage_CL.cpp @@ -0,0 +1,533 @@ +///////////////////////////////////////////////// +// Main author: Chengshun Shang (CIMNE) +// Email: cshang@cimne.upc.edu +// Date: Oct 2023 +///////////////////////////////////////////////// + +// System includes +#include +#include +#include +#include +#include + +// Project includes +#include "DEM_parallel_bond_bilinear_damage_CL.h" +#include "custom_elements/spheric_continuum_particle.h" + +namespace Kratos{ + +DEMContinuumConstitutiveLaw::Pointer DEM_parallel_bond_bilinear_damage::Clone() const{ + DEMContinuumConstitutiveLaw::Pointer p_clone(new DEM_parallel_bond_bilinear_damage(*this)); + return p_clone; +} + +//************************************* +// Parameters preparation +//************************************* + +void DEM_parallel_bond_bilinear_damage::Check(Properties::Pointer pProp) const { + + BaseClassType::Check(pProp); + + if (!pProp->Has(FRACTURE_ENERGY_NORMAL)) { + KRATOS_WARNING("DEM")<<"\nWARNING: Variable FRACTURE_ENERGY_NORMAL was not found in cthe Properties when using DEM_parallel_bond_bilinear_damage. A default value of 0.0 was assigned.\n\n"; + pProp->GetValue(FRACTURE_ENERGY_NORMAL) = 0.0; + } + + if (!pProp->Has(FRACTURE_ENERGY_TANGENTIAL)) { + KRATOS_WARNING("DEM")<<"\nWARNING: Variable FRACTURE_ENERGY_TANGENTIAL was not found in cthe Properties when using DEM_parallel_bond_bilinear_damage. A default value of 0.0 was assigned.\n\n"; + pProp->GetValue(FRACTURE_ENERGY_TANGENTIAL) = 0.0; + } +} + +//************************************* +// Force calculation +//************************************* + +double DEM_parallel_bond_bilinear_damage::ComputeNormalUnbondedForce(double indentation){ + + KRATOS_TRY + + KRATOS_ERROR << "This function shouldn't be accessed here, use basic contact model instead."<mIniNeighbourFailureId[i_neighbour_count]; + const double bonded_indentation = indentation - mInitialIndentationForBondedPart; + + mBondedLocalElasticContactForce2 = 0.0; //normal forces + //double BondedLocalElasticContactForce[2] = {0.0}; + //double current_tangential_force_module = 0.0; + + //const double bond_sigma_max = (*mpProperties)[BOND_SIGMA_MAX]; //tension limit + const double bond_sigma_max = mBondSigmaMax; + const double fracture_energy_normal = (*mpProperties)[FRACTURE_ENERGY_NORMAL]; + //const double delta_at_undamaged_peak_normal = bond_sigma_max * calculation_area / kn_el; + //const double delta_at_failure_point_normal = (2.0 * fracture_energy_normal) / bond_sigma_max; + + const double initial_limit_force = bond_sigma_max * calculation_area; + double DamageEnergyCoeffNormal = 0.0; + + if (bond_sigma_max) { + DamageEnergyCoeffNormal = 2.0 * fracture_energy_normal * kn_el / (calculation_area * bond_sigma_max * bond_sigma_max) - 1.0; + } else { + DamageEnergyCoeffNormal = 0.0; + } + + if (DamageEnergyCoeffNormal > 30.0){ + double suggested_fracture_energy = 31.0 * (calculation_area * bond_sigma_max * bond_sigma_max) / (2.0 * kn_el); + KRATOS_INFO("DEM") << "The normal damage energy is too big! It shouble be smaller than " << suggested_fracture_energy << std::endl; + KRATOS_ERROR<< "The [normal damage energy] is too small!" << std::endl; + } + //KRATOS_ERROR_IF(DamageEnergyCoeffNormal > 30.0) << "Normal damage energy is too big!" << std::endl; + + if (DamageEnergyCoeffNormal < 0.0) { + DamageEnergyCoeffNormal = 0.0; + } + + double k_softening = 0.0; + double limit_force = 0.0; + + if (DamageEnergyCoeffNormal) { + k_softening = kn_el / DamageEnergyCoeffNormal; + } + + double kn_updated = (1.0 - mDamageNormal) * kn_el; + + double delta_accumulated = 0.0; + + double returned_by_mapping_force = 0.0; + + double current_normal_force_module = 0.0; + + if (bonded_indentation >= 0.0) { //COMPRESSION + if (!failure_type) { + mBondedLocalElasticContactForce2 = kn_updated * bonded_indentation; + delta_accumulated = bonded_indentation; + } else { + mBondedLocalElasticContactForce2 = 0.0; + } + } else { //TENSION + + if (!failure_type) { + if (DamageEnergyCoeffNormal) { + limit_force = initial_limit_force * (1.0 + k_softening / kn_el) * kn_updated / (kn_updated + k_softening); + } else { + limit_force = initial_limit_force; + } + + current_normal_force_module = fabs(kn_updated * bonded_indentation); + + mBondedLocalElasticContactForce2 = kn_updated * bonded_indentation; + + delta_accumulated = current_normal_force_module / kn_updated; + + returned_by_mapping_force = current_normal_force_module; + + if ((current_normal_force_module > limit_force) && !(*mpProperties)[IS_UNBREAKABLE]) { + + if (DamageEnergyCoeffNormal) { // the material can sustain further damage, not failure yet + + const double delta_at_undamaged_peak = initial_limit_force / kn_el; + + returned_by_mapping_force = initial_limit_force - k_softening * (delta_accumulated - delta_at_undamaged_peak); + + if (returned_by_mapping_force < 0.0) { + returned_by_mapping_force = 0.0; + } + + mBondedLocalElasticContactForce2 = -returned_by_mapping_force; + + mDamageNormal = 1.0 - (returned_by_mapping_force / delta_accumulated) / kn_el; + + if (mDamageNormal > mDamageThresholdTolerance) { + failure_type = 4; // failure by tension + mBondedLocalElasticContactForce2 = 0.0; + mDamageNormal = 1.0; + } + } else { // Fully fragile behaviour + + failure_type = 4; // failure by tension + mBondedLocalElasticContactForce2 = 0.0; + mDamageNormal = 1.0; + } + } + } else { + mBondedLocalElasticContactForce2 = 0.0; + } + } + + //unbonded force + if (indentation_particle > 0.0) { + mUnbondedLocalElasticContactForce2 = ComputeNormalUnbondedForce(indentation_particle); + } else { + mUnbondedLocalElasticContactForce2 = 0.0; + } + + LocalElasticContactForce[2] = mBondedLocalElasticContactForce2 + mUnbondedLocalElasticContactForce2; + + if (LocalElasticContactForce[2]) { + mBondedScalingFactor[2] = mBondedLocalElasticContactForce2 / LocalElasticContactForce[2]; + } else { + mBondedScalingFactor[2] = 0.0; + } + + KRATOS_CATCH("") + +} // CalculateNormalForces + +void DEM_parallel_bond_bilinear_damage::CalculateViscoDampingCoeff(double &equiv_visco_damp_coeff_normal, + double &equiv_visco_damp_coeff_tangential, + SphericContinuumParticle* element1, + SphericContinuumParticle* element2, + const double kn_el, + const double kt_el) { + + KRATOS_TRY + + const double my_mass = element1->GetMass(); + const double other_mass = element2->GetMass(); + + const double equiv_mass = 1.0 / (1.0/my_mass + 1.0/other_mass); + + const double equiv_gamma = (*mpProperties)[DAMPING_GAMMA]; + + double kn_el_update = (1 - mDamageReal) * kn_el; + double kt_el_update = (1 - mDamageReal) * kt_el; + + equiv_visco_damp_coeff_normal = 2.0 * equiv_gamma * sqrt(equiv_mass * kn_el_update); + equiv_visco_damp_coeff_tangential = 2.0 * equiv_gamma * sqrt(equiv_mass * kt_el_update); + + KRATOS_CATCH("") +} + +void DEM_parallel_bond_bilinear_damage::CalculateTangentialForces(double OldLocalElasticContactForce[3], + double LocalElasticContactForce[3], + double LocalElasticExtraContactForce[3], + double ViscoDampingLocalContactForce[3], + double LocalCoordSystem[3][3], + double LocalDeltDisp[3], + double LocalRelVel[3], + const double kt_el, + const double equiv_shear, + double& contact_sigma, + double& contact_tau, + double indentation_particle, + double calculation_area, + double& failure_criterion_state, + SphericContinuumParticle* element1, + SphericContinuumParticle* element2, + int i_neighbour_count, + bool& sliding, + const ProcessInfo& r_process_info) { + + KRATOS_TRY + + int& failure_type = element1->mIniNeighbourFailureId[i_neighbour_count]; + //const double bond_tau_zero = (*mpProperties)[BOND_TAU_ZERO]; + const double bond_tau_zero = mBondTauZero; + const double internal_friction = (*mpProperties)[BOND_INTERNAL_FRICC]; + const double fracture_energy_tangential = (*mpProperties)[FRACTURE_ENERGY_TANGENTIAL]; + + // The [mBondedScalingFactor] is divided into two direction [0] and [1]. June, 2022 + double OldBondedLocalElasticContactForce[2] = {0.0}; + OldBondedLocalElasticContactForce[0] = mBondedScalingFactor[0] * OldLocalElasticContactForce[0]; + OldBondedLocalElasticContactForce[1] = mBondedScalingFactor[1] * OldLocalElasticContactForce[1]; + + double BondedLocalElasticContactForce[2] = {0.0}; + double current_tangential_force_module = 0.0; + double k_softening = 0.0; + double tau_strength = bond_tau_zero; + double DamageEnergyCoeffTangential = 0.0; + + if (tau_strength) { + DamageEnergyCoeffTangential = 2.0 * fracture_energy_tangential * kt_el / (calculation_area * bond_tau_zero * bond_tau_zero) - 1.0; + } else { + DamageEnergyCoeffTangential = 0.0; + } + + if (DamageEnergyCoeffTangential > 30.0){ + double suggested_fracture_energy = 31.0 * (calculation_area * bond_tau_zero * bond_tau_zero) / (2.0 * kt_el); + KRATOS_INFO("DEM") << "The tangential damage energy is too big! It shouble be smaller than " << suggested_fracture_energy << std::endl; + KRATOS_ERROR<< "The [tangential damage energy] is too small!" << std::endl; + } + //KRATOS_ERROR_IF(DamageEnergyCoeffTangential > 30.0) << "Tangential damage energy is too big!" << std::endl; + + if (DamageEnergyCoeffTangential < 0.0) { + DamageEnergyCoeffTangential = 0.0; + } + + if (DamageEnergyCoeffTangential) { + k_softening = kt_el / DamageEnergyCoeffTangential; + } + + double kt_updated = (1.0 - mDamageTangential) * kt_el; + + if (!failure_type) { + BondedLocalElasticContactForce[0] = OldBondedLocalElasticContactForce[0] - kt_updated * LocalDeltDisp[0]; // 0: first tangential + BondedLocalElasticContactForce[1] = OldBondedLocalElasticContactForce[1] - kt_updated * LocalDeltDisp[1]; // 1: second tangential + } else { + BondedLocalElasticContactForce[0] = 0.0; // 0: first tangential + BondedLocalElasticContactForce[1] = 0.0; // 1: second tangential + } + + double delta_accumulated = 0.0; + double returned_by_mapping_force = 0.0; + + if (!failure_type) { // This means it has not broken yet + + current_tangential_force_module = sqrt(BondedLocalElasticContactForce[0] * BondedLocalElasticContactForce[0] + + BondedLocalElasticContactForce[1] * BondedLocalElasticContactForce[1]); + + delta_accumulated = current_tangential_force_module / kt_updated; + + returned_by_mapping_force = current_tangential_force_module; + + contact_sigma = LocalElasticContactForce[2] / calculation_area; + contact_tau = current_tangential_force_module / calculation_area; + + double updated_max_tau_strength = bond_tau_zero; + + if (contact_sigma >= 0) { + updated_max_tau_strength += internal_friction * contact_sigma; + } + tau_strength = updated_max_tau_strength * (1.0 + k_softening / kt_el) * kt_updated / (kt_updated + k_softening); + + if ((contact_tau > tau_strength) && !(*mpProperties)[IS_UNBREAKABLE]) { // damage + + if (DamageEnergyCoeffTangential) { // the material can sustain further damage, not failure yet + + const double delta_at_undamaged_peak = updated_max_tau_strength * calculation_area / kt_el; + + delta_accumulated = current_tangential_force_module / kt_updated; + + returned_by_mapping_force = updated_max_tau_strength * calculation_area - k_softening * (delta_accumulated - delta_at_undamaged_peak); + + if (returned_by_mapping_force < 0.0) { + returned_by_mapping_force = 0.0; + } + + if (current_tangential_force_module) { + BondedLocalElasticContactForce[0] = (returned_by_mapping_force / current_tangential_force_module) * BondedLocalElasticContactForce[0]; + BondedLocalElasticContactForce[1] = (returned_by_mapping_force / current_tangential_force_module) * BondedLocalElasticContactForce[1]; + } + + mDamageTangential = 1.0 - (returned_by_mapping_force / delta_accumulated) / kt_el; // This is 1 - quotient of K_damaged/K_intact + + if (mDamageTangential > mDamageThresholdTolerance) { + failure_type = 2; // failure by shear + BondedLocalElasticContactForce[0] = 0.0; + BondedLocalElasticContactForce[1] = 0.0; + mDamageTangential = 1.0; + } + } else { // Fully fragile behaviour + + failure_type = 2; // failure by shear + BondedLocalElasticContactForce[0] = 0.0; + BondedLocalElasticContactForce[1] = 0.0; + mDamageTangential = 1.0; + } + } + } + + current_tangential_force_module = sqrt(BondedLocalElasticContactForce[0] * BondedLocalElasticContactForce[0] + + BondedLocalElasticContactForce[1] * BondedLocalElasticContactForce[1]); + + if(calculation_area){ + contact_sigma = mBondedLocalElasticContactForce2 / calculation_area; + contact_tau = current_tangential_force_module / calculation_area; + } + + + //unbonded tangential force + double OldUnbondedLocalElasticContactForce[2] = {0.0}; + double UnbondedLocalElasticContactForce[2] = {0.0}; + double ActualTotalShearForce = 0.0; + double max_admissible_shear_force = 0.0; + double fraction = 0.0; + + if (indentation_particle <= 0.0) { + UnbondedLocalElasticContactForce[0] = 0.0; + UnbondedLocalElasticContactForce[1] = 0.0; + } else { + // Here, unBondedScalingFactor[] = 1 - mBondedScalingFactor[] + OldUnbondedLocalElasticContactForce[0] = (1 - mBondedScalingFactor[0]) * OldLocalElasticContactForce[0]; + OldUnbondedLocalElasticContactForce[1] = (1 - mBondedScalingFactor[1]) * OldLocalElasticContactForce[1]; + + UnbondedLocalElasticContactForce[0] = OldUnbondedLocalElasticContactForce[0] - mKt * LocalDeltDisp[0]; + UnbondedLocalElasticContactForce[1] = OldUnbondedLocalElasticContactForce[1] - mKt * LocalDeltDisp[1]; + + const double& equiv_tg_of_static_fri_ang = (*mpProperties)[STATIC_FRICTION]; + const double& equiv_tg_of_dynamic_fri_ang = (*mpProperties)[DYNAMIC_FRICTION]; + const double& equiv_friction_decay_coefficient = (*mpProperties)[FRICTION_DECAY]; + + const double ShearRelVel = sqrt(LocalRelVel[0] * LocalRelVel[0] + LocalRelVel[1] * LocalRelVel[1]); + // TODO: where this equation from? [equiv_friction] + double equiv_friction = equiv_tg_of_dynamic_fri_ang + (equiv_tg_of_static_fri_ang - equiv_tg_of_dynamic_fri_ang) * exp(-equiv_friction_decay_coefficient * ShearRelVel); + + max_admissible_shear_force = (mUnbondedLocalElasticContactForce2 + mUnbondedViscoDampingLocalContactForce[2]) * equiv_friction; + + if (equiv_tg_of_static_fri_ang < 0.0 || equiv_tg_of_dynamic_fri_ang < 0.0) { + KRATOS_ERROR << "The averaged friction is negative for one contact of element with Id: "<< element1->Id()< max_admissible_shear_force) { + const double ActualElasticShearForce = sqrt(UnbondedLocalElasticContactForce[0] * UnbondedLocalElasticContactForce[0] + + UnbondedLocalElasticContactForce[1] * UnbondedLocalElasticContactForce[1]); + + const double dot_product = UnbondedLocalElasticContactForce[0] * mUnbondedViscoDampingLocalContactForce[0] + + UnbondedLocalElasticContactForce[1] * mUnbondedViscoDampingLocalContactForce[1]; + const double ViscoDampingLocalContactForceModule = sqrt(mUnbondedViscoDampingLocalContactForce[0] * mUnbondedViscoDampingLocalContactForce[0] + + mUnbondedViscoDampingLocalContactForce[1] * mUnbondedViscoDampingLocalContactForce[1]); + + if (dot_product >= 0.0) { + if (ActualElasticShearForce > max_admissible_shear_force) { + // if the ActualElasticShearForce is too big, it should be reduced + if(ActualElasticShearForce){ + fraction = max_admissible_shear_force / ActualElasticShearForce; + } + UnbondedLocalElasticContactForce[0] = UnbondedLocalElasticContactForce[0] * fraction; + UnbondedLocalElasticContactForce[1] = UnbondedLocalElasticContactForce[1] * fraction; + mUnbondedViscoDampingLocalContactForce[0] = 0.0; + mUnbondedViscoDampingLocalContactForce[1] = 0.0; + } + else { + const double ActualViscousShearForce = max_admissible_shear_force - ActualElasticShearForce; + if(ViscoDampingLocalContactForceModule){ + fraction = ActualViscousShearForce / ViscoDampingLocalContactForceModule; + } + mUnbondedViscoDampingLocalContactForce[0] *= fraction; + mUnbondedViscoDampingLocalContactForce[1] *= fraction; + } + } + else { + if (ViscoDampingLocalContactForceModule >= ActualElasticShearForce) { + if(ViscoDampingLocalContactForceModule){ + fraction = (max_admissible_shear_force + ActualElasticShearForce) / ViscoDampingLocalContactForceModule; + } + mUnbondedViscoDampingLocalContactForce[0] *= fraction; + mUnbondedViscoDampingLocalContactForce[1] *= fraction; + } + else { + if(ActualElasticShearForce){ + fraction = max_admissible_shear_force / ActualElasticShearForce; + } + UnbondedLocalElasticContactForce[0] = UnbondedLocalElasticContactForce[0] * fraction; + UnbondedLocalElasticContactForce[1] = UnbondedLocalElasticContactForce[1] * fraction; + mUnbondedViscoDampingLocalContactForce[0] = 0.0; + mUnbondedViscoDampingLocalContactForce[1] = 0.0; + } + } + ViscoDampingLocalContactForce[0] = mBondedViscoDampingLocalContactForce[0] + mUnbondedViscoDampingLocalContactForce[0]; + ViscoDampingLocalContactForce[1] = mBondedViscoDampingLocalContactForce[1] + mUnbondedViscoDampingLocalContactForce[1]; + sliding = true; + } + } + + LocalElasticContactForce[0] = BondedLocalElasticContactForce[0] + UnbondedLocalElasticContactForce[0]; + LocalElasticContactForce[1] = BondedLocalElasticContactForce[1] + UnbondedLocalElasticContactForce[1]; + + // Here, we only calculate the BondedScalingFactor and [unBondedScalingFactor = 1 - BondedScalingFactor]. + if (LocalElasticContactForce[0]) { + mBondedScalingFactor[0] = BondedLocalElasticContactForce[0] / LocalElasticContactForce[0]; + } else { + mBondedScalingFactor[0] = 0.0; + } + + if (LocalElasticContactForce[1]) { + mBondedScalingFactor[1] = BondedLocalElasticContactForce[1] / LocalElasticContactForce[1]; + } else { + mBondedScalingFactor[1] = 0.0; + } + + //for debug + if (mDebugPrintingOption) { + + const long unsigned int& sphere_id = (*mpProperties)[DEBUG_PRINTING_ID_1]; + const long unsigned int& neigh_sphere_id = (*mpProperties)[DEBUG_PRINTING_ID_2]; + + if ((element1->Id() == sphere_id) && (element2->Id() == neigh_sphere_id)) { + std::ofstream normal_forces_file("delta_stress.txt", std::ios_base::out | std::ios_base::app); + normal_forces_file << r_process_info[TIME] << " " << contact_sigma/*1*/ << " " << contact_tau/*2*/ << " " << mDamageNormal/*3*/ << " " << mDamageTangential/*4*/ << " " << '\n'; + normal_forces_file.flush(); + normal_forces_file.close(); + } + } + + //real damage + if (mDamageNormal > mDamageTangential){ + mDamageReal = mDamageNormal; + mDamageTangential = mDamageReal; + } else { + mDamageReal = mDamageTangential; + mDamageNormal = mDamageReal; + } + + KRATOS_CATCH("") +} + +//************************************* +// Moment calculation +//************************************* + +double DEM_parallel_bond_bilinear_damage::GetBondKn(double bond_equiv_young, double calculation_area, double distance){ + KRATOS_TRY + + return ((1 - mDamageReal) * bond_equiv_young * calculation_area / distance); + + KRATOS_CATCH("") +} + + +//************************************* +// Bond failure checking +//************************************* + +void DEM_parallel_bond_bilinear_damage::CheckFailure(const int i_neighbour_count, + SphericContinuumParticle* element1, + SphericContinuumParticle* element2, + double& contact_sigma, + double& contact_tau, + double LocalElasticContactForce[3], + double ViscoDampingLocalContactForce[3], + double ElasticLocalRotationalMoment[3], + double ViscoLocalRotationalMoment[3]){ + + KRATOS_TRY + + //The moments do not contribute to failure in current version + + KRATOS_CATCH("") +}//CheckFailure + + +} //namespace Kratos \ No newline at end of file diff --git a/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_bilinear_damage_CL.h b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_bilinear_damage_CL.h new file mode 100644 index 000000000000..4e061bc90ffa --- /dev/null +++ b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_bilinear_damage_CL.h @@ -0,0 +1,100 @@ +///////////////////////////////////////////////// +// Main author: Chengshun Shang (CIMNE) +// Email: cshang@cimne.upc.edu +// Date: Oct 2023 +///////////////////////////////////////////////// + +#if !defined(DEM_PARALLEL_BOND_BILINEAR_DAMAGE_CL_H_INCLUDE) +#define DEM_PARALLEL_BOND_BILINEAR_DAMAGE_CL_H_INCLUDE + +// Project includes +#include "DEM_parallel_bond_CL.h" + +namespace Kratos{ + + class KRATOS_API(DEM_APPLICATION) DEM_parallel_bond_bilinear_damage : public DEM_parallel_bond { + + typedef DEM_parallel_bond BaseClassType; + + public: + + KRATOS_CLASS_POINTER_DEFINITION(DEM_parallel_bond_bilinear_damage); + + DEM_parallel_bond_bilinear_damage() {} + + void Check(Properties::Pointer pProp) const override; + + ~DEM_parallel_bond_bilinear_damage() {} + + DEMContinuumConstitutiveLaw::Pointer Clone() const override; + double ComputeNormalUnbondedForce(double unbonded_indentation) override; + double GetBondKn(double bond_equiv_young, double calculation_area, double distance) override; + + void CalculateNormalForces(double LocalElasticContactForce[3], + const double kn_el, + double equiv_young, + double indentation, + double indentation_particle, + double calculation_area, + double& acumulated_damage, + SphericContinuumParticle* element1, + SphericContinuumParticle* element2, + int i_neighbour_count, + int time_steps, + const ProcessInfo& r_process_info, + double& contact_sigma) override; + + void CalculateViscoDampingCoeff(double &equiv_visco_damp_coeff_normal, + double &equiv_visco_damp_coeff_tangential, + SphericContinuumParticle* element1, + SphericContinuumParticle* element2, + const double kn_el, + const double kt_el) override; + + void CalculateTangentialForces(double OldLocalElasticContactForce[3], + double LocalElasticContactForce[3], + double LocalElasticExtraContactForce[3], + double ViscoDampingLocalContactForce[3], + double LocalCoordSystem[3][3], + double LocalDeltDisp[3], + double LocalRelVel[3], + const double kt_el, + const double equiv_shear, + double& contact_sigma, + double& contact_tau, + double indentation_particle, + double calculation_area, + double& failure_criterion_state, + SphericContinuumParticle* element1, + SphericContinuumParticle* element2, + int i_neighbour_count, + bool& sliding, + const ProcessInfo& r_process_info) override; + + void CheckFailure(const int i_neighbour_count, + SphericContinuumParticle* element1, + SphericContinuumParticle* element2, + double& contact_sigma, + double& contact_tau, + double LocalElasticContactForce[3], + double ViscoDampingLocalContactForce[3], + double ElasticLocalRotationalMoment[3], + double ViscoLocalRotationalMoment[3]) override; + + double mDamageNormal = 0.0; + double mDamageTangential = 0.0; + double mDamageMoment = 0.0; + const double mDamageThresholdTolerance = 0.9999; + double mDamageReal = 0.0; + bool mDebugPrintingOption = false; + double mInitialIndentationForBondedPart = 0.0; + + protected: + + private: + + }; + +} // namespace Kratos + +#endif /*DEM_PARALLEL_BOND_BILINEAR_DAMAGE_CL_H_INCLUDE defined*/ \ No newline at end of file diff --git a/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_bilinear_damage_mixed_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_bilinear_damage_mixed_CL.cpp new file mode 100644 index 000000000000..09e1e726d6c0 --- /dev/null +++ b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_bilinear_damage_mixed_CL.cpp @@ -0,0 +1,477 @@ +///////////////////////////////////////////////// +// Main author: Chengshun Shang (CIMNE) +// Email: cshang@cimne.upc.edu +// Date: Oct 2023 +///////////////////////////////////////////////// + +//TODO: +//!!!!!!!!!!!!!!!!!!!!!!!!! This model hasn't been well validated! Be careful when you use it + +// System includes +#include +#include +#include +#include +#include + +// Project includes +#include "DEM_parallel_bond_bilinear_damage_mixed_CL.h" +#include "custom_elements/spheric_continuum_particle.h" + +namespace Kratos{ + +DEMContinuumConstitutiveLaw::Pointer DEM_parallel_bond_bilinear_damage_mixed::Clone() const{ + DEMContinuumConstitutiveLaw::Pointer p_clone(new DEM_parallel_bond_bilinear_damage_mixed(*this)); + return p_clone; +} + +void DEM_parallel_bond_bilinear_damage_mixed::Check(Properties::Pointer pProp) const { + + BaseClassType::Check(pProp); + + if (!pProp->Has(FRACTURE_ENERGY_EXPONENT)) { + KRATOS_ERROR<<"\nWARNING: Variable FRACTURE_ENERGY_EXPONENT was not found in cthe Properties when using DEM_parallel_bond_bilinear_damage. A default value of 0.0 was assigned.\n\n"; + pProp->GetValue(FRACTURE_ENERGY_EXPONENT) = 1.0; + } +} + +//************************************* +// Force calculation +//************************************* + +double DEM_parallel_bond_bilinear_damage_mixed::ComputeNormalUnbondedForce(double indentation){ + + KRATOS_TRY + + KRATOS_ERROR << "This function shouldn't be accessed here, use basic contact model instead."<mIniNeighbourFailureId[i_neighbour_count]; + const double bonded_indentation = indentation - mInitialIndentationForBondedPart; + + double BondedLocalElasticContactForce2 = 0.0; //normal forces + // The [mBondedScalingFactor] is divided into two direction [0] and [1]. June, 2022 + double BondedLocalElasticContactForce[2] = {0.0}; + double current_tangential_force_module = 0.0; + + const double bond_sigma_max = (*mpProperties)[BOND_SIGMA_MAX]; //tension limit + const double bond_tau_zero = (*mpProperties)[BOND_TAU_ZERO]; + const double internal_friction = (*mpProperties)[BOND_INTERNAL_FRICC]; + const double fracture_energy_normal = (*mpProperties)[FRACTURE_ENERGY_NORMAL]; + const double fracture_energy_tangential = (*mpProperties)[FRACTURE_ENERGY_TANGENTIAL]; + const double fracture_energy_exponent = (*mpProperties)[FRACTURE_ENERGY_EXPONENT]; + const double delta_at_undamaged_peak_normal = bond_sigma_max * calculation_area / kn_el; + const double delta_at_failure_point_normal = (2.0 * fracture_energy_normal) / bond_sigma_max; + + if (std::abs(delta_at_failure_point_normal) < std::abs(delta_at_undamaged_peak_normal) && fracture_energy_normal){ + double suggested_fracture_energy = delta_at_undamaged_peak_normal * bond_sigma_max / 2.0; + KRATOS_INFO("DEM") << "The [Fracture_energy_normal] is too small! It shouble be bigger than " << suggested_fracture_energy << std::endl; + KRATOS_ERROR<< "The [Fracture_energy_normal] is too small!" << std::endl; + } + + //********************************************** + + if (!failure_type){ //if the bond is not broken + + BondedLocalElasticContactForce2 = (1 - mDamageNormal) * kn_el * bonded_indentation; + //const double OldAccumulatedBondedTangentialLocalDisplacementModulus = sqrt(mAccumulatedBondedTangentialLocalDisplacement[0]*mAccumulatedBondedTangentialLocalDisplacement[0] + mAccumulatedBondedTangentialLocalDisplacement[1]*mAccumulatedBondedTangentialLocalDisplacement[1]); + mAccumulatedBondedTangentialLocalDisplacement[0] += LocalDeltDisp[0]; + mAccumulatedBondedTangentialLocalDisplacement[1] += LocalDeltDisp[1]; + const double AccumulatedBondedTangentialLocalDisplacementModulus = sqrt(mAccumulatedBondedTangentialLocalDisplacement[0]*mAccumulatedBondedTangentialLocalDisplacement[0] + mAccumulatedBondedTangentialLocalDisplacement[1]*mAccumulatedBondedTangentialLocalDisplacement[1]); + + if (bonded_indentation <= 0.0){ // in tension-shear mode + + // See literature [Gao, 2021, https://link.springer.com/article/10.1007/s00603-021-02671-0] + + // Here, we use a mixed mode for considering both mode I (tension) and mode II (shear) + + double delta_equivalent = sqrt(AccumulatedBondedTangentialLocalDisplacementModulus * AccumulatedBondedTangentialLocalDisplacementModulus + bonded_indentation * bonded_indentation); + + double beta = 0.0; + if (bonded_indentation){ + beta = AccumulatedBondedTangentialLocalDisplacementModulus / std::abs(bonded_indentation); + } + double beta_zero = (kn_el * (bond_tau_zero + bond_sigma_max * tan(internal_friction * Globals::Pi / 180.0)))/(bond_sigma_max * kt_el); + + double delta_equivalent_0 = 0.0; + if (beta > beta_zero){ + delta_equivalent_0 = (bond_tau_zero * sqrt(1 + beta * beta)) + / (beta * kt_el / calculation_area + kn_el / calculation_area * tan(internal_friction * Globals::Pi / 180.0)); + } else { + delta_equivalent_0 = delta_at_undamaged_peak_normal * sqrt(1 + beta * beta); + } + + //double deta_equivalent_0_normal = delta_equivalent_0 / sqrt(1 + beta * beta); + //double deta_equivalent_0_tangential = beta * deta_equivalent_0_normal; + double k_bar_equivalent = (kn_el / calculation_area + beta * beta * kt_el / calculation_area) / (1 + beta * beta); + + //yield check + if (fracture_energy_normal && fracture_energy_tangential && !(*mpProperties)[IS_UNBREAKABLE]){ // the material can sustain further damage, not failure yet + + double delta_equivalent_failure = 2 / (k_bar_equivalent * delta_equivalent_0) + * (fracture_energy_normal + (fracture_energy_tangential - fracture_energy_normal) + * std::pow((kt_el / calculation_area * beta * beta)/(kn_el / calculation_area + kt_el / calculation_area * beta * beta), fracture_energy_exponent)); + + if (std::abs(delta_equivalent_failure) < std::abs(delta_equivalent_0)){ + KRATOS_INFO("DEM") << "The [Fracture_energy] is too small! " << std::endl; + KRATOS_ERROR<< "The [Fracture_energy] is too small!" << std::endl; + } + + double mDamageReal_trial = (delta_equivalent_failure * (delta_equivalent - delta_equivalent_0)) / (delta_equivalent * (delta_equivalent_failure - delta_equivalent_0)); + + //real damage calculation + if (mDamageReal_trial > mDamageReal){ + mDamageReal = mDamageReal_trial; + } + if (mDamageReal > 1.0){ + mDamageReal = 1.0; + } + mDamageNormal = mDamageReal; + mDamageTangential = mDamageReal; + + // real force in current step + BondedLocalElasticContactForce2 = (1 - mDamageNormal) * kn_el * bonded_indentation; + BondedLocalElasticContactForce[0] = -1 * (1 - mDamageTangential) * kt_el * mAccumulatedBondedTangentialLocalDisplacement[0]; // 0: first tangential + BondedLocalElasticContactForce[1] = -1 * (1 - mDamageTangential) * kt_el * mAccumulatedBondedTangentialLocalDisplacement[1]; // 1: second tangential + + if (mDamageReal > mDamageThresholdTolerance) { + failure_type = 3; // failure in normal and tangential + BondedLocalElasticContactForce2 = 0.0; + BondedLocalElasticContactForce[0] = 0.0; + BondedLocalElasticContactForce[1] = 0.0; + mDamageReal = 1.0; + } + + } + + if (!(fracture_energy_normal && fracture_energy_tangential) && !(*mpProperties)[IS_UNBREAKABLE]){ // Fully fragile behaviour + + if (delta_equivalent > delta_equivalent_0) { + failure_type = 3; // failure by mixed mode + BondedLocalElasticContactForce2 = 0.0; + BondedLocalElasticContactForce[0] = 0.0; + BondedLocalElasticContactForce[1] = 0.0; + mDamageReal = 1.0; + } + + } + } else if (bonded_indentation > 0.0) { // in compression-shear mode + + // See literature [Gao, 2021, https://link.springer.com/article/10.1007/s00603-021-02671-0] + + double current_sigma = BondedLocalElasticContactForce2 / calculation_area; + + double delta_residual_tangential = current_sigma * tan(internal_friction * Globals::Pi * calculation_area / 180.0) / kt_el; + double delta_at_failure_point_tangential = (2.0 * fracture_energy_tangential) / bond_tau_zero + delta_residual_tangential; + + double max_tau_peak = bond_tau_zero + current_sigma * tan(internal_friction * Globals::Pi / 180.0); + double delta_at_undamaged_peak_tangential = max_tau_peak * calculation_area / kt_el; + + if (std::abs(delta_at_failure_point_tangential) < std::abs(delta_at_undamaged_peak_tangential) && fracture_energy_tangential){ + double suggested_fracture_energy = delta_at_undamaged_peak_tangential * bond_tau_zero / 2.0; + KRATOS_INFO("DEM") << "The [Fracture_energy_tangential] is too small! It shouble be bigger than " << suggested_fracture_energy << std::endl; + KRATOS_ERROR<< "The [Fracture_energy_tangential] is too small!" << std::endl; + } + + double DamageTangentialOnlyForCalculation = 0.0; + if (AccumulatedBondedTangentialLocalDisplacementModulus){ + DamageTangentialOnlyForCalculation = (1 - delta_residual_tangential / AccumulatedBondedTangentialLocalDisplacementModulus) * mDamageTangential; + } else { + DamageTangentialOnlyForCalculation = mDamageTangential; + } + BondedLocalElasticContactForce[0] = -1 * (1 - DamageTangentialOnlyForCalculation) * kt_el * mAccumulatedBondedTangentialLocalDisplacement[0]; // 0: first tangential + BondedLocalElasticContactForce[1] = -1 * (1 - DamageTangentialOnlyForCalculation) * kt_el * mAccumulatedBondedTangentialLocalDisplacement[1]; // 1: second tangential + + current_tangential_force_module = sqrt(BondedLocalElasticContactForce[0] * BondedLocalElasticContactForce[0] + + BondedLocalElasticContactForce[1] * BondedLocalElasticContactForce[1]); + + double current_tau = current_tangential_force_module / calculation_area; + double max_tau = (1 - DamageTangentialOnlyForCalculation) * (bond_tau_zero + current_sigma * tan(internal_friction * Globals::Pi / 180.0)); + + //yield check + if (fracture_energy_normal && fracture_energy_tangential && !(*mpProperties)[IS_UNBREAKABLE]){ // the material can sustain further damage, not failure yet + + if (current_tau > max_tau) { + if ((AccumulatedBondedTangentialLocalDisplacementModulus - delta_residual_tangential) && (delta_at_failure_point_tangential - delta_at_undamaged_peak_tangential)){ + mDamageTangential = ((delta_at_failure_point_tangential - delta_residual_tangential) / (AccumulatedBondedTangentialLocalDisplacementModulus - delta_residual_tangential)) * (AccumulatedBondedTangentialLocalDisplacementModulus - delta_at_undamaged_peak_tangential) / (delta_at_failure_point_tangential - delta_at_undamaged_peak_tangential); + } + } + + //real damage calculation + if (mDamageTangential > mDamageReal){ + mDamageReal = mDamageTangential; + } + if (mDamageReal > 1.0){ + mDamageReal = 1.0; + } + mDamageNormal = mDamageReal; + mDamageTangential = mDamageReal; + + // real force in current step + BondedLocalElasticContactForce2 = (1 - mDamageNormal) * kn_el * bonded_indentation; + if (AccumulatedBondedTangentialLocalDisplacementModulus){ + DamageTangentialOnlyForCalculation = (1 - delta_residual_tangential / AccumulatedBondedTangentialLocalDisplacementModulus) * mDamageTangential; + } else { + DamageTangentialOnlyForCalculation = mDamageTangential; + } + + BondedLocalElasticContactForce[0] = -1 * (1 - DamageTangentialOnlyForCalculation) * kt_el * mAccumulatedBondedTangentialLocalDisplacement[0]; // 0: first tangential + BondedLocalElasticContactForce[1] = -1 * (1 - DamageTangentialOnlyForCalculation) * kt_el * mAccumulatedBondedTangentialLocalDisplacement[1]; // 1: second tangential + + if (mDamageReal > mDamageThresholdTolerance) { + failure_type = 3; // failure in normal and tangential + BondedLocalElasticContactForce2 = 0.0; + BondedLocalElasticContactForce[0] = 0.0; + BondedLocalElasticContactForce[1] = 0.0; + mDamageReal = 1.0; + } + + } + + if (!(fracture_energy_normal && fracture_energy_tangential) && !(*mpProperties)[IS_UNBREAKABLE]){ // Fully fragile behaviour + + if (current_tau > max_tau) { + failure_type = 2; // failure by shear + BondedLocalElasticContactForce[0] = 0.0; + BondedLocalElasticContactForce[1] = 0.0; + mDamageTangential = 1.0; + } + + } + } + + } else { //else the bond is broken + BondedLocalElasticContactForce2 = 0.0; + BondedLocalElasticContactForce[0] = 0.0; + BondedLocalElasticContactForce[1] = 0.0; + } + + if (indentation_particle > 0.0) { + mUnbondedLocalElasticContactForce2 = ComputeNormalUnbondedForce(indentation_particle); + } else { + mUnbondedLocalElasticContactForce2 = 0.0; + } + + current_tangential_force_module = sqrt(BondedLocalElasticContactForce[0] * BondedLocalElasticContactForce[0] + + BondedLocalElasticContactForce[1] * BondedLocalElasticContactForce[1]); + + if(calculation_area){ + contact_sigma = BondedLocalElasticContactForce2 / calculation_area; + contact_tau = current_tangential_force_module / calculation_area; + } + + LocalElasticContactForce[2] = BondedLocalElasticContactForce2 + mUnbondedLocalElasticContactForce2; + + if (LocalElasticContactForce[2]) { + mBondedScalingFactor[2] = BondedLocalElasticContactForce2 / LocalElasticContactForce[2]; + } else { + mBondedScalingFactor[2] = 0.0; + } + + //********************Calculate ViscoDampingCoeff*************************** + if (mDamageNormal < 1.0 || mDamageTangential < 1.0){ + const double my_mass = element1->GetMass(); + const double other_mass = element2->GetMass(); + + const double equiv_mass = 1.0 / (1.0/my_mass + 1.0/other_mass); + + // TODO: check the value of gamma + const double equiv_gamma = (*mpProperties)[DAMPING_GAMMA]; + + double kn_el_update = (1 - mDamageNormal) * kn_el; + double kt_el_update = (1 - mDamageTangential) * kt_el; + + equiv_visco_damp_coeff_normal = 2.0 * equiv_gamma * sqrt(equiv_mass * kn_el_update); + equiv_visco_damp_coeff_tangential = 2.0 * equiv_gamma * sqrt(equiv_mass * kt_el_update); + } else { + equiv_visco_damp_coeff_normal = 0.0; + equiv_visco_damp_coeff_tangential = 0.0; + } + + + //************************Calculate ViscoDamping***************************** + + DEM_parallel_bond::CalculateViscoDamping(LocalRelVel, + ViscoDampingLocalContactForce, + indentation_particle, + equiv_visco_damp_coeff_normal, + equiv_visco_damp_coeff_tangential, + sliding, + element1->mIniNeighbourFailureId[i_neighbour_count], + i_neighbour_count, + element1, + element2); + + // Tangential forces are calculated after the viscodamping because the frictional limit bounds the sum of elastic plus viscous forces + + //************************Calculate tangential forces***************************** + + //unbonded force + double OldUnbondedLocalElasticContactForce[2] = {0.0}; + double UnbondedLocalElasticContactForce[2] = {0.0}; + double ActualTotalShearForce = 0.0; + double max_admissible_shear_force = 0.0; + double fraction = 0.0; + + if (indentation_particle <= 0.0) { + UnbondedLocalElasticContactForce[0] = 0.0; + UnbondedLocalElasticContactForce[1] = 0.0; + } else { + // Here, unBondedScalingFactor[] = 1 - mBondedScalingFactor[] + OldUnbondedLocalElasticContactForce[0] = (1 - mBondedScalingFactor[0]) * OldLocalElasticContactForce[0]; + OldUnbondedLocalElasticContactForce[1] = (1 - mBondedScalingFactor[1]) * OldLocalElasticContactForce[1]; + + UnbondedLocalElasticContactForce[0] = OldUnbondedLocalElasticContactForce[0] - mKt * LocalDeltDisp[0]; + UnbondedLocalElasticContactForce[1] = OldUnbondedLocalElasticContactForce[1] - mKt * LocalDeltDisp[1]; + + const double& equiv_tg_of_static_fri_ang = (*mpProperties)[STATIC_FRICTION]; + const double& equiv_tg_of_dynamic_fri_ang = (*mpProperties)[DYNAMIC_FRICTION]; + const double& equiv_friction_decay_coefficient = (*mpProperties)[FRICTION_DECAY]; + + const double ShearRelVel = sqrt(LocalRelVel[0] * LocalRelVel[0] + LocalRelVel[1] * LocalRelVel[1]); + // TODO: where this equation from? [equiv_friction] + double equiv_friction = equiv_tg_of_dynamic_fri_ang + (equiv_tg_of_static_fri_ang - equiv_tg_of_dynamic_fri_ang) * exp(-equiv_friction_decay_coefficient * ShearRelVel); + + max_admissible_shear_force = (mUnbondedLocalElasticContactForce2 + mUnbondedViscoDampingLocalContactForce[2]) * equiv_friction; + + if (equiv_tg_of_static_fri_ang < 0.0 || equiv_tg_of_dynamic_fri_ang < 0.0) { + KRATOS_ERROR << "The averaged friction is negative for one contact of element with Id: "<< element1->Id()< max_admissible_shear_force) { + const double ActualElasticShearForce = sqrt(UnbondedLocalElasticContactForce[0] * UnbondedLocalElasticContactForce[0] + + UnbondedLocalElasticContactForce[1] * UnbondedLocalElasticContactForce[1]); + + const double dot_product = UnbondedLocalElasticContactForce[0] * mUnbondedViscoDampingLocalContactForce[0] + + UnbondedLocalElasticContactForce[1] * mUnbondedViscoDampingLocalContactForce[1]; + const double ViscoDampingLocalContactForceModule = sqrt(mUnbondedViscoDampingLocalContactForce[0] * mUnbondedViscoDampingLocalContactForce[0] + + mUnbondedViscoDampingLocalContactForce[1] * mUnbondedViscoDampingLocalContactForce[1]); + + if (dot_product >= 0.0) { + if (ActualElasticShearForce > max_admissible_shear_force) { + // if the ActualElasticShearForce is too big, it should be reduced + if(ActualElasticShearForce){ + fraction = max_admissible_shear_force / ActualElasticShearForce; + } + UnbondedLocalElasticContactForce[0] = UnbondedLocalElasticContactForce[0] * fraction; + UnbondedLocalElasticContactForce[1] = UnbondedLocalElasticContactForce[1] * fraction; + mUnbondedViscoDampingLocalContactForce[0] = 0.0; + mUnbondedViscoDampingLocalContactForce[1] = 0.0; + } + else { + const double ActualViscousShearForce = max_admissible_shear_force - ActualElasticShearForce; + if(ViscoDampingLocalContactForceModule){ + fraction = ActualViscousShearForce / ViscoDampingLocalContactForceModule; + } + mUnbondedViscoDampingLocalContactForce[0] *= fraction; + mUnbondedViscoDampingLocalContactForce[1] *= fraction; + } + } + else { + if (ViscoDampingLocalContactForceModule >= ActualElasticShearForce) { + if(ViscoDampingLocalContactForceModule){ + fraction = (max_admissible_shear_force + ActualElasticShearForce) / ViscoDampingLocalContactForceModule; + } + mUnbondedViscoDampingLocalContactForce[0] *= fraction; + mUnbondedViscoDampingLocalContactForce[1] *= fraction; + } + else { + if(ActualElasticShearForce){ + fraction = max_admissible_shear_force / ActualElasticShearForce; + } + UnbondedLocalElasticContactForce[0] = UnbondedLocalElasticContactForce[0] * fraction; + UnbondedLocalElasticContactForce[1] = UnbondedLocalElasticContactForce[1] * fraction; + mUnbondedViscoDampingLocalContactForce[0] = 0.0; + mUnbondedViscoDampingLocalContactForce[1] = 0.0; + } + } + ViscoDampingLocalContactForce[0] = mBondedViscoDampingLocalContactForce[0] + mUnbondedViscoDampingLocalContactForce[0]; + ViscoDampingLocalContactForce[1] = mBondedViscoDampingLocalContactForce[1] + mUnbondedViscoDampingLocalContactForce[1]; + sliding = true; + } + } + + LocalElasticContactForce[0] = BondedLocalElasticContactForce[0] + UnbondedLocalElasticContactForce[0]; + LocalElasticContactForce[1] = BondedLocalElasticContactForce[1] + UnbondedLocalElasticContactForce[1]; + + // Here, we only calculate the BondedScalingFactor and [unBondedScalingFactor = 1 - BondedScalingFactor]. + if (LocalElasticContactForce[0]) { + mBondedScalingFactor[0] = BondedLocalElasticContactForce[0] / LocalElasticContactForce[0]; + } else { + mBondedScalingFactor[0] = 0.0; + } + + if (LocalElasticContactForce[1]) { + mBondedScalingFactor[1] = BondedLocalElasticContactForce[1] / LocalElasticContactForce[1]; + } else { + mBondedScalingFactor[1] = 0.0; + } + + //for debug + if (mDebugPrintingOption) { + + const long unsigned int& sphere_id = (*mpProperties)[DEBUG_PRINTING_ID_1]; + const long unsigned int& neigh_sphere_id = (*mpProperties)[DEBUG_PRINTING_ID_2]; + const double AccumulatedBondedTangentialLocalDisplacementModulus = sqrt(mAccumulatedBondedTangentialLocalDisplacement[0]*mAccumulatedBondedTangentialLocalDisplacement[0] + mAccumulatedBondedTangentialLocalDisplacement[1]*mAccumulatedBondedTangentialLocalDisplacement[1]); + + if ((element1->Id() == sphere_id) && (element2->Id() == neigh_sphere_id)) { + std::ofstream normal_forces_file("delta_stress.txt", std::ios_base::out | std::ios_base::app); + normal_forces_file << r_process_info[TIME] << " " << bonded_indentation/*2*/ << " " << contact_sigma/*3*/ << " " + << AccumulatedBondedTangentialLocalDisplacementModulus/*4*/ << " " << contact_tau/*5*/ << " " + << mDamageReal/*6*/ << " " << mDamageNormal << " " << mDamageTangential << " " << '\n'; + normal_forces_file.flush(); + normal_forces_file.close(); + } + } + + KRATOS_CATCH("") +} + +} //namespace Kratos \ No newline at end of file diff --git a/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_bilinear_damage_mixed_CL.h b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_bilinear_damage_mixed_CL.h new file mode 100644 index 000000000000..66ae89b59b6f --- /dev/null +++ b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_bilinear_damage_mixed_CL.h @@ -0,0 +1,77 @@ +///////////////////////////////////////////////// +// Main author: Chengshun Shang (CIMNE) +// Email: cshang@cimne.upc.edu +// Date: Oct 2023 +///////////////////////////////////////////////// + +#if !defined(DEM_PARALLEL_BOND_BILINEAR_DAMAGE_MIXED_CL_H_INCLUDE) +#define DEM_PARALLEL_BOND_BILINEAR_DAMAGE_MIXED_CL_H_INCLUDE + +// Project includes +#include "DEM_parallel_bond_bilinear_damage_CL.h" + +namespace Kratos{ + + class KRATOS_API(DEM_APPLICATION) DEM_parallel_bond_bilinear_damage_mixed : public DEM_parallel_bond_bilinear_damage { + + typedef DEM_parallel_bond_bilinear_damage BaseClassType; + + public: + + KRATOS_CLASS_POINTER_DEFINITION(DEM_parallel_bond_bilinear_damage_mixed); + + DEM_parallel_bond_bilinear_damage_mixed() {} + + void Check(Properties::Pointer pProp) const override; + + ~DEM_parallel_bond_bilinear_damage_mixed() {} + + DEMContinuumConstitutiveLaw::Pointer Clone() const override; + + double ComputeNormalUnbondedForce(double unbonded_indentation) override; + + void CalculateForces(const ProcessInfo& r_process_info, + double OldLocalElasticContactForce[3], + double LocalElasticContactForce[3], + double LocalElasticExtraContactForce[3], + double LocalCoordSystem[3][3], + double LocalDeltDisp[3], + const double kn_el, + const double kt_el, + double& contact_sigma, + double& contact_tau, + double& failure_criterion_state, + double equiv_young, + double equiv_shear, + double indentation, + double indentation_particle, + double calculation_area, + double& acumulated_damage, + SphericContinuumParticle* element1, + SphericContinuumParticle* element2, + int i_neighbour_count, + int time_steps, + bool& sliding, + double &equiv_visco_damp_coeff_normal, + double &equiv_visco_damp_coeff_tangential, + double LocalRelVel[3], + double ViscoDampingLocalContactForce[3]) override; + + /* + double mDamageNormal = 0.0; + double mDamageTangential = 0.0; + double mDamageMoment = 0.0; + const double mDamageThresholdTolerance = 0.999; + double mDamageReal = 0.0; + bool mDebugPrintingOption = false; + double mInitialIndentationForBondedPart = 0.0;*/ + + protected: + + private: + + }; + +} // namespace Kratos + +#endif /*DEM_PARALLEL_BOND_BILINEAR_DAMAGE_MIXED_CL_H_INCLUDE defined*/ \ No newline at end of file diff --git a/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_for_membrane_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_for_membrane_CL.cpp index 7440d0bdf19d..5f2b9762295c 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_for_membrane_CL.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_for_membrane_CL.cpp @@ -88,6 +88,28 @@ double DEM_parallel_bond_for_membrane::ComputeNormalUnbondedForce(double indenta KRATOS_CATCH("") } +void DEM_parallel_bond_for_membrane::CalculateUnbondedViscoDampingForce(double LocalRelVel[3], + double UnbondedViscoDampingLocalContactForce[3], + SphericParticle* const element1, + SphericParticle* const element2){ + KRATOS_TRY + const double my_mass = element1->GetMass(); + const double other_mass = element2->GetMass(); + + const double equiv_mass = 1.0 / (1.0/my_mass + 1.0/other_mass); + + Properties& properties_of_this_contact = element1->GetProperties().GetSubProperties(element2->GetProperties().Id()); + const double damping_gamma = properties_of_this_contact[DAMPING_GAMMA]; + + const double equiv_visco_damp_coeff_normal = 2.0 * damping_gamma * sqrt(equiv_mass * mKn); + const double equiv_visco_damp_coeff_tangential = 2.0 * damping_gamma * sqrt(equiv_mass * mKt); + + UnbondedViscoDampingLocalContactForce[0] = - equiv_visco_damp_coeff_tangential * LocalRelVel[0]; + UnbondedViscoDampingLocalContactForce[1] = - equiv_visco_damp_coeff_tangential * LocalRelVel[1]; + UnbondedViscoDampingLocalContactForce[2] = - equiv_visco_damp_coeff_normal * LocalRelVel[2]; + KRATOS_CATCH("") +} + //************************************* // Moment calculation diff --git a/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_for_membrane_CL.h b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_for_membrane_CL.h index 929efb077107..836b86bb5123 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_for_membrane_CL.h +++ b/applications/DEMApplication/custom_constitutive/DEM_parallel_bond_for_membrane_CL.h @@ -32,6 +32,11 @@ namespace Kratos{ double ComputeNormalUnbondedForce(double unbonded_indentation) override; + void CalculateUnbondedViscoDampingForce(double LocalRelVel[3], + double UnbondedViscoDampingLocalContactForce[3], + SphericParticle* const element1, + SphericParticle* const element2) override; + void ComputeParticleRotationalMoments(SphericContinuumParticle* element, SphericContinuumParticle* neighbor, double equiv_young, diff --git a/applications/DEMApplication/custom_constitutive/DEM_smooth_joint_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_smooth_joint_CL.cpp index d7d031b50f61..b8e6383d84bf 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_smooth_joint_CL.cpp +++ b/applications/DEMApplication/custom_constitutive/DEM_smooth_joint_CL.cpp @@ -224,21 +224,25 @@ void DEM_smooth_joint::CalculateElasticConstants(double& kn_el, double& kt_el, d kn_el = (*mpProperties)[JOINT_NORMAL_STIFFNESS]; kt_el = (*mpProperties)[JOINT_TANGENTIAL_STIFFNESS]; - double GlobalJointNormal[3] = {0.0}; - GlobalJointNormal[0] = (*mpProperties)[JOINT_NORMAL_DIRECTION_X]; - GlobalJointNormal[1] = (*mpProperties)[JOINT_NORMAL_DIRECTION_Y]; - GlobalJointNormal[2] = (*mpProperties)[JOINT_NORMAL_DIRECTION_Z]; + mGlobalJointNormal[0] = (*mpProperties)[JOINT_NORMAL_DIRECTION_X]; // X = -1 * sin(alpha) + mGlobalJointNormal[1] = (*mpProperties)[JOINT_NORMAL_DIRECTION_Y]; // Y = cos(alpha) + mGlobalJointNormal[2] = (*mpProperties)[JOINT_NORMAL_DIRECTION_Z]; // Z = 0 array_1d OtherToMeVector; noalias(OtherToMeVector) = element1->GetGeometry()[0].Coordinates() - element2->GetGeometry()[0].Coordinates(); double LocalCoordSystem[3][3]; double Distance = DEM_MODULUS_3(OtherToMeVector); GeometryFunctions::ComputeContactLocalCoordSystem(OtherToMeVector, Distance, LocalCoordSystem); - GeometryFunctions::VectorGlobal2Local(LocalCoordSystem, GlobalJointNormal, mLocalJointNormal); - + GeometryFunctions::VectorGlobal2Local(LocalCoordSystem, mGlobalJointNormal, mLocalJointNormal); KRATOS_CATCH("") } +void DEM_smooth_joint::GetGlobalJointNormal(double GlobalJointNormal[3]){ + GlobalJointNormal[0] = mGlobalJointNormal[0]; + GlobalJointNormal[1] = mGlobalJointNormal[1]; + GlobalJointNormal[2] = mGlobalJointNormal[2]; +} + double DEM_smooth_joint::LocalMaxSearchDistance(const int i, SphericContinuumParticle* element1, SphericContinuumParticle* element2) { @@ -299,6 +303,7 @@ void DEM_smooth_joint::CalculateForces(const ProcessInfo& r_process_info, double equiv_young, double equiv_shear, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -317,6 +322,7 @@ void DEM_smooth_joint::CalculateForces(const ProcessInfo& r_process_info, kn_el, equiv_young, indentation, + indentation_particle, calculation_area, acumulated_damage, element1, @@ -346,7 +352,52 @@ void DEM_smooth_joint::CalculateForces(const ProcessInfo& r_process_info, sliding, r_process_info, time_steps); + + /* + bool is_normal_force_smaller_than_zero = false; + if (LocalElasticContactForce[2] < 0.0){ + is_normal_force_smaller_than_zero = true; + } + + double temp_GlobalElasticContactForce[3] = {0.0}; + GeometryFunctions::VectorLocal2Global(LocalCoordSystem, LocalElasticContactForce, temp_GlobalElasticContactForce); + + double temp_GlobalElasticContactForce_Normal = 0.0; + temp_GlobalElasticContactForce_Normal = GeometryFunctions::DotProduct(temp_GlobalElasticContactForce, mGlobalJointNormal); + double temp_GlobalElasticContactForce_Tangential[2] = {0.0}; + + //temp_GlobalElasticContactForce_Tangential[0] = temp_GlobalElasticContactForce[0] - temp_GlobalElasticContactForce_Normal * mGlobalJointNormal[0]; + //temp_GlobalElasticContactForce_Tangential[1] = temp_GlobalElasticContactForce[2] - temp_GlobalElasticContactForce_Normal * mGlobalJointNormal[2]; + temp_GlobalElasticContactForce_Tangential[0] = temp_GlobalElasticContactForce[0]; + temp_GlobalElasticContactForce_Tangential[1] = temp_GlobalElasticContactForce[2]; + + //LocalElasticContactForce is already a global vector + LocalElasticContactForce[0] = temp_GlobalElasticContactForce_Tangential[0]; + LocalElasticContactForce[1] = temp_GlobalElasticContactForce_Normal; + LocalElasticContactForce[2] = temp_GlobalElasticContactForce_Tangential[1]; + + if(calculation_area){ + contact_sigma = temp_GlobalElasticContactForce_Normal / calculation_area; + + bool is_global_normal_force_smaller_than_zero = false; + if (temp_GlobalElasticContactForce_Normal < 0.0){ + is_global_normal_force_smaller_than_zero = true; + } + if ((is_normal_force_smaller_than_zero == false && is_global_normal_force_smaller_than_zero == true) || (is_normal_force_smaller_than_zero == true && is_global_normal_force_smaller_than_zero == false)) { + contact_sigma *= -1; + } + } + + double temp_current_tangential_force_module = 0.0; + temp_current_tangential_force_module = sqrt(temp_GlobalElasticContactForce_Tangential[0] * temp_GlobalElasticContactForce_Tangential[0] + + temp_GlobalElasticContactForce_Tangential[1] * temp_GlobalElasticContactForce_Tangential[1]); + + if (calculation_area){ + contact_tau = temp_current_tangential_force_module / calculation_area; + } + */ + KRATOS_CATCH("") } @@ -354,6 +405,7 @@ void DEM_smooth_joint::CalculateNormalForces(double LocalElasticContactForce[3], const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -448,7 +500,7 @@ void DEM_smooth_joint::CalculateTangentialForces(double OldLocalElasticContactFo JointLocalElasticContactForce[1] -= kt_el * mAccumulatedJointTangentialLocalDisplacement[1]; // 1: second tangential current_tangential_force_module = sqrt(JointLocalElasticContactForce[0] * JointLocalElasticContactForce[0] + JointLocalElasticContactForce[1] * JointLocalElasticContactForce[1]); - double friction_force = 0.5 * LocalElasticContactForce[2]; + double friction_force = (*mpProperties)[JOINT_FRICTION_COEFF] * LocalElasticContactForce[2]; if (current_tangential_force_module > friction_force){ if (current_tangential_force_module > 0.0){ const double fraction = friction_force / current_tangential_force_module; @@ -458,9 +510,13 @@ void DEM_smooth_joint::CalculateTangentialForces(double OldLocalElasticContactFo } } - current_tangential_force_module = sqrt(JointLocalElasticContactForce[0] * JointLocalElasticContactForce[0] - + JointLocalElasticContactForce[1] * JointLocalElasticContactForce[1]); + if (mGlobalJointNormal[0] == 0.0 && mGlobalJointNormal[1] == 1.0 && mGlobalJointNormal[2] == 0.0) { + JointLocalElasticContactForce[0] = 0.0; + JointLocalElasticContactForce[1] = 0.0; + } + current_tangential_force_module = sqrt(JointLocalElasticContactForce[0] * JointLocalElasticContactForce[0] + + JointLocalElasticContactForce[1] * JointLocalElasticContactForce[1]); if (calculation_area){ contact_tau = current_tangential_force_module / calculation_area; @@ -485,7 +541,8 @@ void DEM_smooth_joint::CalculateMoments(SphericContinuumParticle* element, double ElasticLocalRotationalMoment[3], double ViscoLocalRotationalMoment[3], double equiv_poisson, - double indentation, + double indentation, + double indentation_particle, double LocalElasticContactForce[3], double normalLocalContactForce, double GlobalElasticContactForces[3], @@ -549,7 +606,10 @@ void DEM_smooth_joint::CheckFailure(const int i_neighbour_count, failure_type = 2; // failure in shear contact_sigma = 0.0; contact_tau = 0.0; + LocalElasticContactForce[0] = 0.0; + LocalElasticContactForce[1] = 0.0; LocalElasticContactForce[2] = 0.0; + /* double current_tangential_force_module = sqrt(LocalElasticContactForce[0] * LocalElasticContactForce[0] + LocalElasticContactForce[1] * LocalElasticContactForce[1]); double friction_force = (*mpProperties)[JOINT_FRICTION_COEFF] * LocalElasticContactForce[2]; @@ -559,10 +619,14 @@ void DEM_smooth_joint::CheckFailure(const int i_neighbour_count, LocalElasticContactForce[0] *= fraction; LocalElasticContactForce[1] *= fraction; } + } else { + LocalElasticContactForce[0] = 0.0; + LocalElasticContactForce[1] = 0.0; } - } - - } + LocalElasticContactForce[2] = 0.0; + */ + } + } KRATOS_CATCH("") }//CheckFailure diff --git a/applications/DEMApplication/custom_constitutive/DEM_smooth_joint_CL.h b/applications/DEMApplication/custom_constitutive/DEM_smooth_joint_CL.h index 0b129238fb6c..ab04d1f4871f 100644 --- a/applications/DEMApplication/custom_constitutive/DEM_smooth_joint_CL.h +++ b/applications/DEMApplication/custom_constitutive/DEM_smooth_joint_CL.h @@ -39,6 +39,7 @@ namespace Kratos{ void CalculateElasticConstants(double& kn_el, double& kt_el, double initial_dist, double equiv_young, double equiv_poisson, double calculation_area, SphericContinuumParticle* element1, SphericContinuumParticle* element2, double indentation) override; //virtual void InitializeContact(SphericParticle* const element1, SphericParticle* const element2, const double indentation); + void GetGlobalJointNormal(double GlobalJointNormal[3]) override; double LocalMaxSearchDistance(const int i, SphericContinuumParticle* element1, @@ -70,6 +71,7 @@ namespace Kratos{ double equiv_young, double equiv_shear, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -86,6 +88,7 @@ namespace Kratos{ const double kn_el, double equiv_young, double indentation, + double indentation_particle, double calculation_area, double& acumulated_damage, SphericContinuumParticle* element1, @@ -124,7 +127,8 @@ namespace Kratos{ double ElasticLocalRotationalMoment[3], double ViscoLocalRotationalMoment[3], double equiv_poisson, - double indentation, + double indentation, + double indentation_particle, double LocalElasticContactForce[3], double normalLocalContactForce, double GlobalElasticContactForces[3], @@ -145,6 +149,7 @@ namespace Kratos{ double mKn; double mKt; double mLocalJointNormal[3] = {0.0}; + double mGlobalJointNormal[3] = {0.0}; double mInitialDistanceJoint = 0.0; protected: diff --git a/applications/DEMApplication/custom_elements/spheric_continuum_particle.cpp b/applications/DEMApplication/custom_elements/spheric_continuum_particle.cpp index 9ac0eb9d8edc..29da1b4e2311 100644 --- a/applications/DEMApplication/custom_elements/spheric_continuum_particle.cpp +++ b/applications/DEMApplication/custom_elements/spheric_continuum_particle.cpp @@ -241,6 +241,22 @@ namespace Kratos { double initial_dist = radius_sum - initial_delta; double indentation = initial_dist - data_buffer.mDistance; + double indentation_particle = radius_sum - data_buffer.mDistance; + + if (this->mIndentationInitialOption) { + + if (data_buffer.mTime < (0.0 + 2.0 * data_buffer.mDt)){ + if (indentation_particle > 0.0) { + this->mIndentationInitial[data_buffer.mpOtherParticle->Id()] = indentation_particle; + } + } + + if (this->mIndentationInitial.find(data_buffer.mpOtherParticle->Id()) != this->mIndentationInitial.end()){ + indentation_particle -= this->mIndentationInitial[data_buffer.mpOtherParticle->Id()]; + } + + } + double myYoung = GetYoung(); double myPoisson = GetPoisson(); @@ -270,7 +286,11 @@ namespace Kratos { if (i < (int)mContinuumInitialNeighborsSize) { mContinuumConstitutiveLawArray[i]->GetContactArea(GetRadius(), other_radius, cont_ini_neigh_area, i, calculation_area); //some Constitutive Laws get a value, some others calculate the value. - mContinuumConstitutiveLawArray[i]->CalculateElasticConstants(kn_el, kt_el, initial_dist, equiv_young, equiv_poisson, calculation_area, this, neighbour_iterator, indentation); + if (mContinuumConstitutiveLawArray[i]->GetTypeOfLaw() == "parallel_bond_CL") { + mContinuumConstitutiveLawArray[i]->CalculateElasticConstants(kn_el, kt_el, initial_dist, equiv_young, equiv_poisson, calculation_area, this, neighbour_iterator, indentation_particle); + } else { + mContinuumConstitutiveLawArray[i]->CalculateElasticConstants(kn_el, kt_el, initial_dist, equiv_young, equiv_poisson, calculation_area, this, neighbour_iterator, indentation); + } } EvaluateDeltaDisplacement(data_buffer, DeltDisp, RelVel, data_buffer.mLocalCoordSystem, data_buffer.mOldLocalCoordSystem, vel, delta_displ); @@ -315,7 +335,7 @@ namespace Kratos { double GlobalContactForce[3] = {0.0}; GeometryFunctions::VectorGlobal2Local(data_buffer.mLocalCoordSystem, RelVel, LocalRelVel); - + //*******************Forces calculation**************** if (i < (int)mContinuumInitialNeighborsSize) { @@ -333,6 +353,7 @@ namespace Kratos { equiv_young, equiv_shear, indentation, + indentation_particle, calculation_area, acumulated_damage, this, @@ -345,11 +366,11 @@ namespace Kratos { LocalRelVel, ViscoDampingLocalContactForce); - } else if (indentation > 0.0) { - const double previous_indentation = indentation + LocalDeltDisp[2]; + } else if (indentation_particle > 0.0) { + const double previous_indentation = indentation_particle + LocalDeltDisp[2]; mDiscontinuumConstitutiveLaw = pCloneDiscontinuumConstitutiveLawWithNeighbour(data_buffer.mpOtherParticle); mDiscontinuumConstitutiveLaw->CalculateForces(r_process_info, OldLocalElasticContactForce, LocalElasticContactForce, - LocalDeltDisp, LocalRelVel, indentation, previous_indentation, + LocalDeltDisp, LocalRelVel, indentation_particle, previous_indentation, ViscoDampingLocalContactForce, cohesive_force, this, data_buffer.mpOtherParticle, sliding, data_buffer.mLocalCoordSystem); } else { //Not bonded and no idata_buffer.mpOtherParticlendentation LocalElasticContactForce[0] = 0.0; LocalElasticContactForce[1] = 0.0; LocalElasticContactForce[2] = 0.0; @@ -361,13 +382,26 @@ namespace Kratos { ComputeOtherBallToBallForces(other_ball_to_ball_forces); if (i < (int)mContinuumInitialNeighborsSize && mContinuumConstitutiveLawArray[i]->GetTypeOfLaw() == "smooth_joint_CL") { - double JointLocalCoordSystem[3][3]; - GeometryFunctions::RotateCoordToDirection(data_buffer.mLocalCoordSystem, mLocalJointNormal, JointLocalCoordSystem); - GeometryFunctions::VectorLocal2Global(JointLocalCoordSystem, LocalElasticContactForce, GlobalElasticContactForce); + //double JointLocalCoordSystem[3][3]; + double GlobalJointNormal[3]; + mContinuumConstitutiveLawArray[i]->GetGlobalJointNormal(GlobalJointNormal); + //GeometryFunctions::RotateCoordToDirection(data_buffer.mLocalCoordSystem, LocalJointNormal, JointLocalCoordSystem); + //GeometryFunctions::VectorLocal2Global(JointLocalCoordSystem, LocalElasticContactForce, GlobalElasticContactForce); + GeometryFunctions::VectorLocal2Global(data_buffer.mLocalCoordSystem, LocalElasticContactForce, GlobalElasticContactForce); + //GlobalElasticContactForce[0] = LocalElasticContactForce[0]; //LocalElasticContactForce is acturlly a global force for "smooth_joint_CL" + //GlobalElasticContactForce[1] = LocalElasticContactForce[1]; + //GlobalElasticContactForce[2] = LocalElasticContactForce[2]; + GeometryFunctions::RotateVectorToVector(GlobalElasticContactForce, GlobalJointNormal, GlobalElasticContactForce); + } else { GeometryFunctions::VectorLocal2Global(data_buffer.mLocalCoordSystem, LocalElasticContactForce, GlobalElasticContactForce); } + //*******************Add up forces and moments************** + AddUpForcesAndProject(data_buffer.mOldLocalCoordSystem, data_buffer.mLocalCoordSystem, LocalContactForce, LocalElasticContactForce, LocalElasticExtraContactForce, GlobalContactForce, + GlobalElasticContactForce, GlobalElasticExtraContactForce, TotalGlobalElasticContactForce, ViscoDampingLocalContactForce, 0.0, other_ball_to_ball_forces, rElasticForce, rContactForce, i, r_process_info); //TODO: replace the 0.0 with an actual cohesive force for discontinuum neighbours + + //******************Moments calculation start**************** if (this->Is(DEMFlags::HAS_ROTATION)) { if (i < (int)mContinuumInitialNeighborsSize) { @@ -381,16 +415,17 @@ namespace Kratos { ViscoLocalRotationalMoment, equiv_poisson, indentation, - LocalElasticContactForce, + indentation_particle, + LocalContactForce, LocalContactForce[2], - GlobalElasticContactForce, + GlobalContactForce, data_buffer.mLocalCoordSystem[2], i); if (this->Is(DEMFlags::HAS_ROLLING_FRICTION) && !data_buffer.mMultiStageRHS) { if (mRollingFrictionModel->CheckIfThisModelRequiresRecloningForEachNeighbour()){ mRollingFrictionModel = pCloneRollingFrictionModelWithNeighbour(data_buffer.mpOtherParticle); - mRollingFrictionModel->ComputeRollingFriction(this, data_buffer.mpOtherParticle, r_process_info, LocalContactForce, indentation, mContactMoment); + mRollingFrictionModel->ComputeRollingFriction(this, data_buffer.mpOtherParticle, r_process_info, LocalContactForce, indentation_particle, mContactMoment); } else { if ((i >= (int)mContinuumInitialNeighborsSize) || mIniNeighbourFailureId[i]) { @@ -401,15 +436,12 @@ namespace Kratos { } else { //for unbonded particles - double GlobalElasticContactForce[3] = {0.0}; - GeometryFunctions::VectorLocal2Global(data_buffer.mLocalCoordSystem, LocalElasticContactForce, GlobalElasticContactForce); - ComputeMoments(LocalContactForce[2], GlobalElasticContactForce, data_buffer.mLocalCoordSystem[2], data_buffer.mpOtherParticle, indentation, i); - + ComputeMoments(LocalContactForce[2], GlobalContactForce, data_buffer.mLocalCoordSystem[2], data_buffer.mpOtherParticle, indentation_particle, i); if (this->Is(DEMFlags::HAS_ROLLING_FRICTION) && !data_buffer.mMultiStageRHS) { if (mRollingFrictionModel->CheckIfThisModelRequiresRecloningForEachNeighbour()){ mRollingFrictionModel = pCloneRollingFrictionModelWithNeighbour(data_buffer.mpOtherParticle); - mRollingFrictionModel->ComputeRollingFriction(this, data_buffer.mpOtherParticle, r_process_info, LocalContactForce, indentation, mContactMoment); + mRollingFrictionModel->ComputeRollingFriction(this, data_buffer.mpOtherParticle, r_process_info, LocalContactForce, indentation_particle, mContactMoment); } else { if ((i >= (int)mContinuumInitialNeighborsSize) || mIniNeighbourFailureId[i]) { @@ -432,14 +464,11 @@ namespace Kratos { ViscoDampingLocalContactForce, ElasticLocalRotationalMoment, ViscoLocalRotationalMoment); } - //*******************Add up forces and moments************** - AddUpForcesAndProject(data_buffer.mOldLocalCoordSystem, data_buffer.mLocalCoordSystem, LocalContactForce, LocalElasticContactForce, LocalElasticExtraContactForce, GlobalContactForce, - GlobalElasticContactForce, GlobalElasticExtraContactForce, TotalGlobalElasticContactForce, ViscoDampingLocalContactForce, 0.0, other_ball_to_ball_forces, rElasticForce, rContactForce, i, r_process_info); //TODO: replace the 0.0 with an actual cohesive force for discontinuum neighbours - if (this->Is(DEMFlags::HAS_ROTATION)) { AddUpMomentsAndProject(data_buffer.mLocalCoordSystem, ElasticLocalRotationalMoment, ViscoLocalRotationalMoment); } + if (r_process_info[CONTACT_MESH_OPTION] == 1 && (i < (int)mContinuumInitialNeighborsSize) && this->Id() < neighbour_iterator_id) { double total_local_elastic_contact_force[3] = {0.0}; total_local_elastic_contact_force[0] = LocalElasticContactForce[0] + LocalElasticExtraContactForce[0]; diff --git a/applications/DEMApplication/custom_elements/spheric_continuum_particle.h b/applications/DEMApplication/custom_elements/spheric_continuum_particle.h index c4039b8a0882..bf29b6e25c52 100644 --- a/applications/DEMApplication/custom_elements/spheric_continuum_particle.h +++ b/applications/DEMApplication/custom_elements/spheric_continuum_particle.h @@ -134,7 +134,7 @@ namespace Kratos unsigned int mInitialNeighborsSize; std::vector mContinuumConstitutiveLawArray; double mLocalRadiusAmplificationFactor = 1.0; - double mLocalJointNormal[3]; + //double mLocalJointNormal[3]; protected: diff --git a/applications/DEMApplication/custom_elements/spheric_particle.cpp b/applications/DEMApplication/custom_elements/spheric_particle.cpp index 89e06d445067..d558d560c9bb 100644 --- a/applications/DEMApplication/custom_elements/spheric_particle.cpp +++ b/applications/DEMApplication/custom_elements/spheric_particle.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include @@ -117,6 +118,7 @@ SphericParticle::~SphericParticle(){ SphericParticle& SphericParticle::operator=(const SphericParticle& rOther) { DiscreteElement::operator=(rOther); mElasticEnergy = rOther.mElasticEnergy; + mMaxNormalBallToBallForceTimesRadius = rOther.mMaxNormalBallToBallForceTimesRadius; mInelasticFrictionalEnergy = rOther.mInelasticFrictionalEnergy; mInelasticViscodampingEnergy = rOther.mInelasticViscodampingEnergy; mInelasticRollingResistanceEnergy = rOther.mInelasticRollingResistanceEnergy; @@ -189,6 +191,7 @@ void SphericParticle::Initialize(const ProcessInfo& r_process_info) KRATOS_TRY this->mInitializationTime = r_process_info[TIME]; + this->mIndentationInitialOption = r_process_info[CLEAN_INDENT_V2_OPTION]; SetValue(NEIGHBOUR_IDS, DenseVector()); @@ -250,6 +253,8 @@ void SphericParticle::Initialize(const ProcessInfo& r_process_info) inelastic_viscodamping_energy = 0.0; double& inelastic_rollingresistance_energy = this->GetInelasticRollingResistanceEnergy(); inelastic_rollingresistance_energy = 0.0; + double& max_normal_ball_to_ball_force_times_radius = this->GetMaxNormalBallToBallForceTimesRadius(); + max_normal_ball_to_ball_force_times_radius = 0.0; DEMIntegrationScheme::Pointer& translational_integration_scheme = GetProperties()[DEM_TRANSLATIONAL_INTEGRATION_SCHEME_POINTER]; DEMIntegrationScheme::Pointer& rotational_integration_scheme = GetProperties()[DEM_ROTATIONAL_INTEGRATION_SCHEME_POINTER]; @@ -280,6 +285,7 @@ void SphericParticle::CalculateRightHandSide(const ProcessInfo& r_process_info, NodeType& this_node = GetGeometry()[0]; data_buffer.mDt = dt; + data_buffer.mTime = r_process_info[TIME]; data_buffer.mMultiStageRHS = false; array_1d additional_forces = ZeroVector(3); @@ -544,6 +550,7 @@ void SphericParticle::EvaluateDeltaDisplacement(ParticleDataBuffer & data_buffer RelDeltDisp[0] = (delta_displ[0] - other_delta_displ[0]); RelDeltDisp[1] = (delta_displ[1] - other_delta_displ[1]); RelDeltDisp[2] = (delta_displ[2] - other_delta_displ[2]); + } void SphericParticle::RelativeDisplacementAndVelocityOfContactPointDueToRotation(const double indentation, @@ -597,6 +604,7 @@ void SphericParticle::RelativeDisplacementAndVelocityOfContactPointDueToRotation RelDeltDisp[0] += my_delta_disp_at_contact_point_due_to_rotation[0] - other_delta_disp_at_contact_point_due_to_rotation[0]; RelDeltDisp[1] += my_delta_disp_at_contact_point_due_to_rotation[1] - other_delta_disp_at_contact_point_due_to_rotation[1]; RelDeltDisp[2] += my_delta_disp_at_contact_point_due_to_rotation[2] - other_delta_disp_at_contact_point_due_to_rotation[2]; + } @@ -748,6 +756,7 @@ void SphericParticle::RelativeDisplacementAndVelocityOfContactPointDueToRotation RelDeltDisp[0] += (my_new_arm_vector[0] - other_new_arm_vector[0]) + (other_arm_vector[0] - my_arm_vector[0]); RelDeltDisp[1] += (my_new_arm_vector[1] - other_new_arm_vector[1]) + (other_arm_vector[1] - my_arm_vector[1]); RelDeltDisp[2] += (my_new_arm_vector[2] - other_new_arm_vector[2]) + (other_arm_vector[2] - my_arm_vector[2]); + } void SphericParticle::ComputeMoments(double NormalLocalContactForce, @@ -768,6 +777,7 @@ void SphericParticle::ComputeMoments(double NormalLocalContactForce, array_1d moment_of_this_neighbour; GeometryFunctions::CrossProduct(arm_vector, Force, moment_of_this_neighbour); noalias(mContactMoment) += moment_of_this_neighbour; + } void SphericParticle::ComputeMomentsWithWalls(double NormalLocalContactForce, @@ -882,6 +892,13 @@ void SphericParticle::ComputeBallToBallContactForceAndMoment(SphericParticle::Pa } } + if (r_process_info[ENERGY_CALCULATION_OPTION]){ + double NormalBallToBallForceTimesRadius = this->GetRadius() * LocalContactForce[2]; + if ( NormalBallToBallForceTimesRadius > mMaxNormalBallToBallForceTimesRadius) { + mMaxNormalBallToBallForceTimesRadius = NormalBallToBallForceTimesRadius; + } + } + // Store contact information needed for later processes StoreBallToBallContactInfo(r_process_info, data_buffer, GlobalContactForce, LocalContactForce, ViscoDampingLocalContactForce, sliding); @@ -926,7 +943,7 @@ void SphericParticle::EvaluateBallToBallForcesForPositiveIndentiations(SphericPa data_buffer.mLocalRelVel[1] = 0.0; data_buffer.mLocalRelVel[2] = 0.0; GeometryFunctions::VectorGlobal2Local(LocalCoordSystem, RelVel, data_buffer.mLocalRelVel); - + mDiscontinuumConstitutiveLaw = pCloneDiscontinuumConstitutiveLawWithNeighbour(p_neighbour_element); mDiscontinuumConstitutiveLaw->CalculateForces(r_process_info, OldLocalElasticContactForce, LocalElasticContactForce, LocalDeltDisp, data_buffer.mLocalRelVel, indentation, previous_indentation, @@ -987,6 +1004,20 @@ void SphericParticle::ComputeBallToRigidFaceContactForceAndMoment(SphericParticl if (ContactType == 1 || ContactType == 2 || ContactType == 3) { double indentation = -(DistPToB - GetInteractionRadius()) - ini_delta; + + if (this->mIndentationInitialOption) { + + if (data_buffer.mTime < (0.0 + 2.0 * data_buffer.mDt)){ + if (indentation > 0.0) { + this->mIndentationInitialWall[static_cast(rNeighbours[i]->Id())] = indentation; + } + } + + if (this->mIndentationInitialWall.find(static_cast(rNeighbours[i]->Id())) != this->mIndentationInitialWall.end()){ + indentation -= this->mIndentationInitialWall[static_cast(rNeighbours[i]->Id())]; + } + } + double DeltDisp[3] = {0.0}; double DeltVel [3] = {0.0}; @@ -1459,6 +1490,8 @@ void SphericParticle::InitializeSolutionStep(const ProcessInfo& r_process_info) central_node.FastGetSolutionStepValue(REPRESENTATIVE_VOLUME) = CalculateVolume(); double& elastic_energy = this->GetElasticEnergy(); elastic_energy = 0.0; + double& max_normal_ball_to_ball_force_times_radius = this->GetMaxNormalBallToBallForceTimesRadius(); + max_normal_ball_to_ball_force_times_radius = 0.0; if (this->Is(DEMFlags::HAS_STRESS_TENSOR)) { for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { @@ -1756,6 +1789,7 @@ void SphericParticle::ComputeAdditionalForces(array_1d& externally_ap const double vel_magnitude = DEM_MODULUS_3(vel); if (vel_magnitude != 0.0) { + mGlobalViscousDamping = r_process_info[GLOBAL_VISCOUS_DAMPING]; //update the coefficient of GLOBAL_VISCOUS_DAMPING const array_1d viscous_damping_force = -2.0 * mGlobalViscousDamping * sqrt(GetMass() * GetRadius() * GetYoung()) * vel; noalias(externally_applied_force) += viscous_damping_force; } @@ -1831,8 +1865,9 @@ void SphericParticle::AddUpMomentsAndProject(double LocalCoordSystem[3][3], } GeometryFunctions::VectorLocal2Global(LocalCoordSystem, LocalContactRotationalMoment, GlobalContactRotationalMoment); - + DEM_ADD_SECOND_TO_FIRST(mContactMoment, GlobalContactRotationalMoment) + } void SphericParticle::AddUpFEMForcesAndProject(double LocalCoordSystem[3][3], @@ -1976,6 +2011,12 @@ void SphericParticle::Calculate(const Variable& rVariable, double& Outpu } + if (rVariable == PARTICLE_MAX_NORMAL_BALL_TO_BALL_FORCE_TIMES_RADIUS) { + + Output = GetMaxNormalBallToBallForceTimesRadius(); + + } + if (rVariable == PARTICLE_INELASTIC_FRICTIONAL_ENERGY) { Output = GetInelasticFrictionalEnergy(); @@ -2060,7 +2101,7 @@ bool SphericParticle::CalculateRelativePositionsOrSkipContact(ParticleDataBuffer data_buffer.mOtherToMeVector[1] = data_buffer.mMyCoors[1] - data_buffer.mOtherCoors[1]; data_buffer.mOtherToMeVector[2] = data_buffer.mMyCoors[2] - data_buffer.mOtherCoors[2]; - data_buffer.mDistance = DEM_MODULUS_3(data_buffer.mOtherToMeVector); + data_buffer.mDistance = DEM_MODULUS_3(data_buffer.mOtherToMeVector); must_skip_contact_calculation = data_buffer.mDistance < std::numeric_limits::epsilon(); @@ -2072,6 +2113,19 @@ bool SphericParticle::CalculateRelativePositionsOrSkipContact(ParticleDataBuffer data_buffer.mRadiusSum = this->GetInteractionRadius() + data_buffer.mOtherRadius; data_buffer.mIndentation = data_buffer.mRadiusSum - data_buffer.mDistance; + if (this->mIndentationInitialOption) { + + if (data_buffer.mTime < (0.0 + 2.0 * data_buffer.mDt)){ + if (data_buffer.mIndentation > 0.0) { + this->mIndentationInitial[data_buffer.mpOtherParticle->Id()] = data_buffer.mIndentation; + } + } + + if (this->mIndentationInitial.find(data_buffer.mpOtherParticle->Id()) != this->mIndentationInitial.end()){ + data_buffer.mIndentation -= this->mIndentationInitial[data_buffer.mpOtherParticle->Id()]; + } + } + return data_buffer.mIndentation > 0.0; } @@ -2196,6 +2250,19 @@ double SphericParticle::GetRadius() double SphericParticle::CalculateVolume() { return 4.0 * Globals::Pi / 3.0 * mRadius * mRadius * mRadius; } void SphericParticle::SetRadius(double radius) { mRadius = radius; } void SphericParticle::SetRadius() { mRadius = GetGeometry()[0].FastGetSolutionStepValue(RADIUS); } +void SphericParticle::SetRadius(bool is_radius_expansion, double radius_expansion_rate, double radius_multiplier_max, double radius_multiplier, double radius_multiplier_old) +{ + if (is_radius_expansion){ + if (radius_multiplier_old >= 1.0) { + mRadius = (GetGeometry()[0].FastGetSolutionStepValue(RADIUS) / radius_multiplier_old) * radius_multiplier; + GetGeometry()[0].FastGetSolutionStepValue(RADIUS) = mRadius; + } else { + mRadius = GetGeometry()[0].FastGetSolutionStepValue(RADIUS); + } + } else { + mRadius = GetGeometry()[0].FastGetSolutionStepValue(RADIUS); + } +} double SphericParticle::GetInteractionRadius(const int radius_index) { return mRadius; } void SphericParticle::SetInteractionRadius(const double radius, const int radius_index) { mRadius = radius; GetGeometry()[0].FastGetSolutionStepValue(RADIUS) = radius;} double SphericParticle::GetSearchRadius() { return mSearchRadius; } @@ -2219,6 +2286,7 @@ double& SphericParticle::GetElasticEnergy() double& SphericParticle::GetInelasticFrictionalEnergy() { return mInelasticFrictionalEnergy; } double& SphericParticle::GetInelasticViscodampingEnergy() { return mInelasticViscodampingEnergy; } double& SphericParticle::GetInelasticRollingResistanceEnergy() { return mInelasticRollingResistanceEnergy; } +double& SphericParticle::GetMaxNormalBallToBallForceTimesRadius() { return mMaxNormalBallToBallForceTimesRadius; } void SphericParticle::SetYoungFromProperties(double* young) { GetFastProperties()->SetYoungFromProperties( young); } void SphericParticle::SetPoissonFromProperties(double* poisson) { GetFastProperties()->SetPoissonFromProperties( poisson); } diff --git a/applications/DEMApplication/custom_elements/spheric_particle.h b/applications/DEMApplication/custom_elements/spheric_particle.h index 44da78ce145d..1053d440d885 100644 --- a/applications/DEMApplication/custom_elements/spheric_particle.h +++ b/applications/DEMApplication/custom_elements/spheric_particle.h @@ -8,6 +8,7 @@ // System includes #include #include +#include // Project includes #include "includes/define.h" @@ -95,6 +96,7 @@ bool mDomainIsPeriodic; double mDistance; double mRadiusSum; double mDt; +double mTime; double mOtherRadius; double mIndentation; double mMyCoors[3]; @@ -189,6 +191,7 @@ void SetProgrammedDestructionTime(const double time); virtual double GetRadius(); virtual void SetRadius(double radius); virtual void SetRadius(); +virtual void SetRadius(bool is_radius_expansion, double radius_expansion_rate, double radius_multiplier_max, double radius_multiplier, double radius_multiplier_old); virtual double CalculateVolume(); virtual double GetInteractionRadius(const int radius_index = 0); virtual void SetInteractionRadius(const double radius, const int radius_index = 0); @@ -214,6 +217,7 @@ virtual double& GetElasticEnergy(); virtual double& GetInelasticFrictionalEnergy(); virtual double& GetInelasticViscodampingEnergy(); virtual double& GetInelasticRollingResistanceEnergy(); +virtual double& GetMaxNormalBallToBallForceTimesRadius(); PropertiesProxy* GetFastProperties(); void SetFastProperties(PropertiesProxy* pProps); @@ -243,6 +247,7 @@ double mInelasticFrictionalEnergy; double mInelasticViscodampingEnergy; double mInelasticRollingResistanceEnergy; double mPartialRepresentativeVolume; +double mMaxNormalBallToBallForceTimesRadius; std::vector mBondElements; std::vector mNeighbourElements; @@ -438,6 +443,9 @@ std::unique_ptr mDiscontinuumConstitutiveLaw; std::unique_ptr mRollingFrictionModel; double mInitializationTime; +double mIndentationInitialOption; +std::map mIndentationInitial; +std::map mIndentationInitialWall; double mProgrammedDestructionTime=-1.0; // set to a negative value, so that when marked TO_ERASE, elimination is by default. double mRadius; double mSearchRadius; @@ -465,6 +473,7 @@ virtual void save(Serializer& rSerializer) const override rSerializer.save("mInelasticViscodampingEnergy", mInelasticViscodampingEnergy); rSerializer.save("mInelasticRollingResistanceEnergy", mInelasticRollingResistanceEnergy); rSerializer.save("mPartialRepresentativeVolume", mPartialRepresentativeVolume); + rSerializer.save("mMaxNormalBallToBallForceTimesRadius", mMaxNormalBallToBallForceTimesRadius); rSerializer.save("mBondElements", mBondElements); rSerializer.save("mNeighbourElements", mNeighbourElements); rSerializer.save("mContactingNeighbourIds", mContactingNeighbourIds); @@ -509,6 +518,7 @@ virtual void load(Serializer& rSerializer) override rSerializer.load("mInelasticViscodampingEnergy", mInelasticViscodampingEnergy); rSerializer.load("mInelasticRollingResistanceEnergy", mInelasticRollingResistanceEnergy); rSerializer.load("mPartialRepresentativeVolume", mPartialRepresentativeVolume); + rSerializer.load("mMaxNormalBallToBallForceTimesRadius", mMaxNormalBallToBallForceTimesRadius); rSerializer.load("mBondElements", mBondElements); rSerializer.load("mNeighbourElements", mNeighbourElements); rSerializer.load("mContactingNeighbourIds", mContactingNeighbourIds); diff --git a/applications/DEMApplication/custom_python/DEM_python_application.cpp b/applications/DEMApplication/custom_python/DEM_python_application.cpp index 89decaf01681..06f8a37aac23 100644 --- a/applications/DEMApplication/custom_python/DEM_python_application.cpp +++ b/applications/DEMApplication/custom_python/DEM_python_application.cpp @@ -77,6 +77,7 @@ PYBIND11_MODULE(KratosDEMApplication,m) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, LOCAL_COORDINATION_NUMBER_OPTION) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, GLOBAL_COORDINATION_NUMBER_OPTION) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, CLEAN_INDENT_OPTION) + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, CLEAN_INDENT_V2_OPTION) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, TRIHEDRON_OPTION) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, ROLLING_FRICTION_OPTION) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, POISSON_EFFECT_OPTION) @@ -100,6 +101,7 @@ PYBIND11_MODULE(KratosDEMApplication,m) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, DEM_DRAG_CONSTANT_X) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, DEM_DRAG_CONSTANT_Y) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, DEM_DRAG_CONSTANT_Z) + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, ENERGY_CALCULATION_OPTION) //KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, CLUSTER_INFORMATION) //Dangerous. Requires adding ClusterInformation to python KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, CLUSTER_FILE_NAME) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, INJECTOR_ELEMENT_TYPE) @@ -211,6 +213,9 @@ PYBIND11_MODULE(KratosDEMApplication,m) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, DEBUG_PRINTING_ID_1) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, DEBUG_PRINTING_ID_2) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, FRACTURE_ENERGY) + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, FRACTURE_ENERGY_NORMAL) + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, FRACTURE_ENERGY_TANGENTIAL) + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, FRACTURE_ENERGY_EXPONENT) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, SIGMA_SLOPE_CHANGE_THRESHOLD) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, INTERNAL_FRICTION_AFTER_THRESHOLD) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, SEARCH_RADIUS_INCREMENT_FOR_BONDS_CREATION) @@ -273,6 +278,12 @@ PYBIND11_MODULE(KratosDEMApplication,m) KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS(m, LOCAL_AUX_ANGULAR_VELOCITY) // ******************* Quaternion Integration END ******************* + // ****************Radius expansion method BEGIN******************* + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, IS_RADIUS_EXPANSION) + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, RADIUS_EXPANSION_RATE) + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, RADIUS_MULTIPLIER_MAX) + // *****************Radius expansion method END******************** + // FORCE AND MOMENTUM KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS(m, PARTICLE_MOMENT) @@ -296,6 +307,7 @@ PYBIND11_MODULE(KratosDEMApplication,m) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PARTICLE_INELASTIC_VISCODAMPING_ENERGY) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PARTICLE_INELASTIC_FRICTIONAL_ENERGY) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PARTICLE_INELASTIC_ROLLING_RESISTANCE_ENERGY) + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PARTICLE_MAX_NORMAL_BALL_TO_BALL_FORCE_TIMES_RADIUS) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, COMPUTE_ENERGY_OPTION) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, GLOBAL_DAMPING) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, GLOBAL_VISCOUS_DAMPING) diff --git a/applications/DEMApplication/custom_python/add_custom_constitutive_laws_to_python.cpp b/applications/DEMApplication/custom_python/add_custom_constitutive_laws_to_python.cpp index dc4cff84edfb..a18abf488aeb 100644 --- a/applications/DEMApplication/custom_python/add_custom_constitutive_laws_to_python.cpp +++ b/applications/DEMApplication/custom_python/add_custom_constitutive_laws_to_python.cpp @@ -62,6 +62,8 @@ #include "custom_constitutive/DEM_parallel_bond_CL.h" #include "custom_constitutive/DEM_smooth_joint_CL.h" #include "custom_constitutive/DEM_parallel_bond_for_membrane_CL.h" +#include "custom_constitutive/DEM_parallel_bond_bilinear_damage_CL.h" +#include "custom_constitutive/DEM_parallel_bond_bilinear_damage_mixed_CL.h" #include "custom_constitutive/DEM_rolling_friction_model_constant_torque.h" #include "custom_constitutive/DEM_rolling_friction_model_viscous_torque.h" #include "custom_constitutive/DEM_rolling_friction_model_bounded.h" @@ -271,14 +273,22 @@ void AddCustomConstitutiveLawsToPython(pybind11::module& m) { ; py::class_(m, "DEM_parallel_bond") - .def(py::init<>()) + .def(py::init<>()) ; py::class_(m, "DEM_smooth_joint") .def(py::init<>()) ; - py::class_(m, "DEM_parallel_bond_for_membrane") + py::class_(m, "DEM_parallel_bond_for_membrane") + .def(py::init<>()) + ; + + py::class_(m, "DEM_parallel_bond_bilinear_damage") + .def(py::init<>()) + ; + + py::class_(m, "DEM_parallel_bond_bilinear_damage_mixed") .def(py::init<>()) ; @@ -294,6 +304,30 @@ void AddCustomConstitutiveLawsToPython(pybind11::module& m) { .def(py::init<>()) ; + py::class_, DEM_compound_constitutive_law_for_PBM::Pointer, DEM_parallel_bond_bilinear_damage>(m, "DEM_parallel_bond_bilinear_damage_Hertz") + .def(py::init<>()) + ; + + py::class_, DEM_compound_constitutive_law_for_PBM::Pointer, DEM_parallel_bond_bilinear_damage>(m, "DEM_parallel_bond_bilinear_damage_Linear") + .def(py::init<>()) + ; + + py::class_, DEM_compound_constitutive_law_for_PBM::Pointer, DEM_parallel_bond_bilinear_damage>(m, "DEM_parallel_bond_bilinear_damage_Quadratic") + .def(py::init<>()) + ; + + py::class_, DEM_compound_constitutive_law_for_PBM::Pointer, DEM_parallel_bond_bilinear_damage_mixed>(m, "DEM_parallel_bond_bilinear_damage_mixed_Hertz") + .def(py::init<>()) + ; + + py::class_, DEM_compound_constitutive_law_for_PBM::Pointer, DEM_parallel_bond_bilinear_damage_mixed>(m, "DEM_parallel_bond_bilinear_damage_mixed_Linear") + .def(py::init<>()) + ; + + py::class_, DEM_compound_constitutive_law_for_PBM::Pointer, DEM_parallel_bond_bilinear_damage_mixed>(m, "DEM_parallel_bond_bilinear_damage_mixed_Quadratic") + .def(py::init<>()) + ; + //for compound constitutive law diff --git a/applications/DEMApplication/custom_python/add_custom_strategies_to_python.cpp b/applications/DEMApplication/custom_python/add_custom_strategies_to_python.cpp index ef2bc26c7bf1..8a268d5f11b0 100755 --- a/applications/DEMApplication/custom_python/add_custom_strategies_to_python.cpp +++ b/applications/DEMApplication/custom_python/add_custom_strategies_to_python.cpp @@ -161,6 +161,7 @@ void AddCustomStrategiesToPython(pybind11::module& m) .def("SearchFemNeighbours", &ExplicitSolverStrategy::SearchFEMOperations) .def("SetSearchRadiiOnAllParticles", &ExplicitSolverStrategy::SetSearchRadiiOnAllParticles) .def("SetNormalRadiiOnAllParticles", &ExplicitSolverStrategy::SetNormalRadiiOnAllParticles) + .def("SetNormalRadiiOnAllParticlesBeforeInitilization", &ExplicitSolverStrategy::SetNormalRadiiOnAllParticlesBeforeInitilization) .def("SetSearchRadiiWithFemOnAllParticles", &ExplicitSolverStrategy::SetSearchRadiiWithFemOnAllParticles) .def("RebuildListOfDiscontinuumSphericParticles", &ExplicitSolverStrategy::RebuildListOfDiscontinuumSphericParticles) .def("PrepareElementsForPrinting", &ExplicitSolverStrategy::PrepareElementsForPrinting) diff --git a/applications/DEMApplication/custom_python/add_custom_utilities_to_python.cpp b/applications/DEMApplication/custom_python/add_custom_utilities_to_python.cpp index 8f271acbf4a9..e9a5931550a8 100755 --- a/applications/DEMApplication/custom_python/add_custom_utilities_to_python.cpp +++ b/applications/DEMApplication/custom_python/add_custom_utilities_to_python.cpp @@ -236,6 +236,7 @@ void AddCustomUtilitiesToPython(pybind11::module& m) { .def("CalculateTotalMomentum", &SphericElementGlobalPhysicsCalculator::CalculateTotalMomentum) .def("CalulateTotalAngularMomentum", &SphericElementGlobalPhysicsCalculator::CalulateTotalAngularMomentum) .def("CalculateSumOfInternalForces", &SphericElementGlobalPhysicsCalculator::CalculateSumOfInternalForces) + .def("CalculateParticleNumberTimesMaxNormalBallToBallForceTimesRadius", &SphericElementGlobalPhysicsCalculator::CalculateParticleNumberTimesMaxNormalBallToBallForceTimesRadius) ; void (DemSearchUtilities::*SearchNodeNeigboursDistancesMM)(ModelPart&,ModelPart&,const double&,const Variable&) = &DemSearchUtilities::SearchNodeNeigboursDistances >; diff --git a/applications/DEMApplication/custom_strategies/strategies/explicit_solver_continuum.cpp b/applications/DEMApplication/custom_strategies/strategies/explicit_solver_continuum.cpp index 0ad1fe59c514..02e2e547aa9e 100644 --- a/applications/DEMApplication/custom_strategies/strategies/explicit_solver_continuum.cpp +++ b/applications/DEMApplication/custom_strategies/strategies/explicit_solver_continuum.cpp @@ -97,6 +97,11 @@ namespace Kratos { //RebuildListsOfPointersOfEachParticle(); //Serialized pointers are lost, so we rebuild them using Id's } + // Finding overlapping of initial configurations + if (r_process_info[CLEAN_INDENT_OPTION]) { + for (int i = 0; i < 10; i++) CalculateInitialMaxIndentations(r_process_info); + } + if (fem_model_part.Nodes().size() > 0) { SetSearchRadiiWithFemOnAllParticles(r_model_part, mpDem_model_part->GetProcessInfo()[SEARCH_RADIUS_INCREMENT_FOR_WALLS], 1.0); SearchRigidFaceNeighbours(); @@ -126,11 +131,6 @@ namespace Kratos { ComputeNewNeighboursHistoricalData(); } - - - - - AttachSpheresToStickyWalls(); if (r_process_info[CONTACT_MESH_OPTION] == 1) { diff --git a/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.cpp b/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.cpp index bdd51609204e..cc45c8c5ceb7 100644 --- a/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.cpp +++ b/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.cpp @@ -631,7 +631,7 @@ namespace Kratos { RebuildListOfSphericParticles(r_model_part.GetCommunicator().LocalMesh().Elements(), mListOfSphericParticles); - SetNormalRadiiOnAllParticles(*mpDem_model_part); + SetNormalRadiiOnAllParticlesBeforeInitilization(*mpDem_model_part); #pragma omp parallel { @@ -1307,7 +1307,7 @@ namespace Kratos { KRATOS_CATCH("") } - void ExplicitSolverStrategy::SetNormalRadiiOnAllParticles(ModelPart& r_model_part) { + void ExplicitSolverStrategy::SetNormalRadiiOnAllParticlesBeforeInitilization(ModelPart& r_model_part) { KRATOS_TRY int number_of_elements = r_model_part.GetCommunicator().LocalMesh().ElementsArray().end() - r_model_part.GetCommunicator().LocalMesh().ElementsArray().begin(); @@ -1317,6 +1317,30 @@ namespace Kratos { KRATOS_CATCH("") } + + void ExplicitSolverStrategy::SetNormalRadiiOnAllParticles(ModelPart& r_model_part) { + KRATOS_TRY + int number_of_elements = r_model_part.GetCommunicator().LocalMesh().ElementsArray().end() - r_model_part.GetCommunicator().LocalMesh().ElementsArray().begin(); + + ProcessInfo& r_process_info = GetModelPart().GetProcessInfo(); + bool is_radius_expansion = r_process_info[IS_RADIUS_EXPANSION]; + double radius_expansion_rate = r_process_info[RADIUS_EXPANSION_RATE]; + double radius_multiplier_max = r_process_info[RADIUS_MULTIPLIER_MAX]; + const double time = r_process_info[TIME]; + const double delta_time = r_process_info[DELTA_TIME]; + double radius_multiplier = 1.0 + time * radius_expansion_rate; + double radius_multiplier_old = 1.0 + (time - delta_time) * radius_expansion_rate; + + if (radius_multiplier > radius_multiplier_max) { + is_radius_expansion = false; + } + + IndexPartition(number_of_elements).for_each([&](unsigned int i){ + mListOfSphericParticles[i]->SetRadius(is_radius_expansion, radius_expansion_rate, radius_multiplier_max, radius_multiplier, radius_multiplier_old); + }); + + KRATOS_CATCH("") + } void ExplicitSolverStrategy::SetSearchRadiiWithFemOnAllParticles(ModelPart& r_model_part, const double added_search_distance, const double amplification) { KRATOS_TRY diff --git a/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.h b/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.h index d960eab40d83..de193c3b5bf1 100644 --- a/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.h +++ b/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.h @@ -233,6 +233,7 @@ namespace Kratos { void ApplyPrescribedBoundaryConditions(); void ApplyInitialConditions(); virtual void SetSearchRadiiOnAllParticles(ModelPart& r_model_part, const double added_search_distance = 0.0, const double amplification = 1.0); + void SetNormalRadiiOnAllParticlesBeforeInitilization(ModelPart& r_model_part); void SetNormalRadiiOnAllParticles(ModelPart& r_model_part); virtual void SetSearchRadiiWithFemOnAllParticles(ModelPart& r_model_part, const double added_search_distance = 0.0, const double amplification = 1.0); virtual void SearchNeighbours(); diff --git a/applications/DEMApplication/custom_utilities/GeometryFunctions.h b/applications/DEMApplication/custom_utilities/GeometryFunctions.h index 160d685ff027..ff5e79f053e5 100644 --- a/applications/DEMApplication/custom_utilities/GeometryFunctions.h +++ b/applications/DEMApplication/custom_utilities/GeometryFunctions.h @@ -150,7 +150,9 @@ namespace Kratos { for (int i=0; i<3; i++) { LocalVector[i] = 0.0; for (int j=0; j<3; j++) { - LocalVector[i]+=LocalCoordSystem[i][j]*GlobalVector[j]; + long double temp_1 = LocalCoordSystem[i][j]; + long double temp_2 = GlobalVector[j]; + LocalVector[i] += temp_1 * temp_2; } } } @@ -160,7 +162,9 @@ namespace Kratos { for (int i=0; i<3; i++) { LocalVector[i] = 0.0; for (int j=0; j<3; j++) { - LocalVector[i]+=LocalCoordSystem[i][j]*GlobalVector[j]; + long double temp_1 = LocalCoordSystem[i][j]; + long double temp_2 = GlobalVector[j]; + LocalVector[i] += temp_1 * temp_2; } } } @@ -170,7 +174,9 @@ namespace Kratos { for (int i=0; i<3; i++) { LocalVector[i] = 0.0; for (int j=0; j<3; j++) { - LocalVector[i]+=LocalCoordSystem[i][j]*GlobalVector[j]; + long double temp_1 = LocalCoordSystem[i][j]; + long double temp_2 = GlobalVector[j]; + LocalVector[i] += temp_1 * temp_2; } } } @@ -180,7 +186,9 @@ namespace Kratos { for (int i=0; i<3; i++) { GlobalVector[i] = 0.0; for (int j=0; j<3; j++) { - GlobalVector[i]+=LocalCoordSystem[j][i]*LocalVector[j]; + long double temp_1 = LocalCoordSystem[j][i]; + long double temp_2 = LocalVector[j]; + GlobalVector[i] += temp_1 * temp_2; } } } @@ -190,7 +198,9 @@ namespace Kratos { for (int i=0; i<3; i++) { GlobalVector[i] = 0.0; for (int j=0; j<3; j++) { - GlobalVector[i]+=LocalCoordSystem[j][i]*LocalVector[j]; + long double temp_1 = LocalCoordSystem[j][i]; + long double temp_2 = LocalVector[j]; + GlobalVector[i] += temp_1 * temp_2; } } } @@ -200,7 +210,9 @@ namespace Kratos { for (int i=0; i<3; i++) { GlobalVector[i] = 0.0; for (int j=0; j<3; j++) { - GlobalVector[i]+=LocalCoordSystem[j][i]*LocalVector[j]; + long double temp_1 = LocalCoordSystem[j][i]; + long double temp_2 = LocalVector[j]; + GlobalVector[i] += temp_1 * temp_2; } } } @@ -210,7 +222,9 @@ namespace Kratos { for (int i=0; i<3; i++) { GlobalVector[i] = 0.0; for (int j=0; j<3; j++) { - GlobalVector[i]+=LocalCoordSystem[j][i]*LocalVector[j]; + long double temp_1 = LocalCoordSystem[j][i]; + long double temp_2 = LocalVector[j]; + GlobalVector[i] += temp_1 * temp_2; } } } @@ -936,6 +950,7 @@ namespace Kratos { }*/ + /* -------------Not verified function static inline void RotateCoordToDirection(const double OLdCoordSystem[3][3], double Vector[3], double NewCoordSystem[3][3]) { double x_axis[3] = {0.0}; @@ -972,6 +987,83 @@ namespace Kratos { NewCoordSystem[2][1] = rotation_matrix[1][0] * x_axis[1] + rotation_matrix[1][1] * y_axis[1] + rotation_matrix[1][2] * z_axis[1]; NewCoordSystem[2][2] = rotation_matrix[1][0] * x_axis[2] + rotation_matrix[1][1] * y_axis[2] + rotation_matrix[1][2] * z_axis[2]; } + */ + + static inline void RotateVectorToVector(double InitialVector[3], double TargetVector[3], double NewVector[3]) + { + if (InitialVector[0] == 0.0 && InitialVector[1] == 0.0 && InitialVector[2] == 0.0) { + NewVector[0] = 0.0; + NewVector[1] = 0.0; + NewVector[2] = 0.0; + } else { + double InitialVector_BeforeNormalize[3] = {0.0}; + InitialVector_BeforeNormalize[0] = InitialVector[0]; + InitialVector_BeforeNormalize[1] = InitialVector[1]; + InitialVector_BeforeNormalize[2] = InitialVector[2]; + + normalize(InitialVector); + normalize(TargetVector); + + double InitialVector_AfterNormalize[3] = {0.0}; + InitialVector_AfterNormalize[0] = InitialVector[0]; + InitialVector_AfterNormalize[1] = InitialVector[1]; + InitialVector_AfterNormalize[2] = InitialVector[2]; + + double TargetVector_AfterNormalize[3] = {0.0}; + TargetVector_AfterNormalize[0] = TargetVector[0]; + TargetVector_AfterNormalize[1] = TargetVector[1]; + TargetVector_AfterNormalize[2] = TargetVector[2]; + + //Pay attention here! This function will rotate the InitialVector to the direction of TargetVector if the angle between them is smaller than 90 degree. Otherwise, rotating to the opposite direction. + if (DotProduct(InitialVector_AfterNormalize, TargetVector_AfterNormalize) < 0.0){ + TargetVector_AfterNormalize[0] = -1 * TargetVector_AfterNormalize[0]; + TargetVector_AfterNormalize[1] = -1 * TargetVector_AfterNormalize[1]; + TargetVector_AfterNormalize[2] = -1 * TargetVector_AfterNormalize[2]; + } + + double rotation_angle = acos(DotProduct(InitialVector_AfterNormalize, TargetVector_AfterNormalize)); + + // Calculate the cross product + std::vector crossProduct; + crossProduct.push_back(InitialVector_AfterNormalize[1] * TargetVector_AfterNormalize[2] - InitialVector_AfterNormalize[2] * TargetVector_AfterNormalize[1]); + crossProduct.push_back(InitialVector_AfterNormalize[2] * TargetVector_AfterNormalize[0] - InitialVector_AfterNormalize[0] * TargetVector_AfterNormalize[2]); + crossProduct.push_back(InitialVector_AfterNormalize[0] * TargetVector_AfterNormalize[1] - InitialVector_AfterNormalize[1] * TargetVector_AfterNormalize[0]); + + double norm_cross_product = std::sqrt(crossProduct[0] * crossProduct[0] + crossProduct[1] * crossProduct[1] + crossProduct[2] * crossProduct[2]); + crossProduct[0] = crossProduct[0] / norm_cross_product; + crossProduct[1] = crossProduct[1] / norm_cross_product; + crossProduct[2] = crossProduct[2] / norm_cross_product; + + + // Calculate the rotation matrix + double c = cos(rotation_angle); + double s = sin(rotation_angle); + double t = 1 - c; + double x = crossProduct[0]; + double y = crossProduct[1]; + double z = crossProduct[2]; + + // Apply the rotation matrix to the vector + std::vector rotatedVector; + rotatedVector.push_back(t * x * x + c); + rotatedVector.push_back(t * x * y - s * z); + rotatedVector.push_back(t * x * z + s * y); + + rotatedVector.push_back(t * x * y + s * z); + rotatedVector.push_back(t * y * y + c); + rotatedVector.push_back(t * y * z - s * x); + + rotatedVector.push_back(t * x * z - s * y); + rotatedVector.push_back(t * y * z + s * x); + rotatedVector.push_back(t * z * z + c); + + // Apply the rotation to the vector + NewVector[0] = rotatedVector[0] * InitialVector_BeforeNormalize[0] + rotatedVector[1] * InitialVector_BeforeNormalize[1] + rotatedVector[2] * InitialVector_BeforeNormalize[2]; + NewVector[1] = rotatedVector[3] * InitialVector_BeforeNormalize[0] + rotatedVector[4] * InitialVector_BeforeNormalize[1] + rotatedVector[5] * InitialVector_BeforeNormalize[2]; + NewVector[2] = rotatedVector[6] * InitialVector_BeforeNormalize[0] + rotatedVector[7] * InitialVector_BeforeNormalize[1] + rotatedVector[8] * InitialVector_BeforeNormalize[2]; + } + + } static inline bool InsideOutside(const array_1d& Coord1, const array_1d& Coord2, diff --git a/applications/DEMApplication/custom_utilities/calculate_global_physical_properties.h b/applications/DEMApplication/custom_utilities/calculate_global_physical_properties.h index 675c664a38ea..b10be275cd4c 100644 --- a/applications/DEMApplication/custom_utilities/calculate_global_physical_properties.h +++ b/applications/DEMApplication/custom_utilities/calculate_global_physical_properties.h @@ -425,6 +425,32 @@ class SphericElementGlobalPhysicsCalculator return rollingresistance_energy; } + double CalculateParticleNumberTimesMaxNormalBallToBallForceTimesRadius(ModelPart& r_model_part) + { + OpenMPUtils::CreatePartition(ParallelUtilities::GetNumThreads(), r_model_part.GetCommunicator().LocalMesh().Elements().size(), mElementsPartition); + + double particle_max_normal_ball_to_ball_force_times_radius = 0.0; + + #pragma omp parallel for reduction(+ : particle_max_normal_ball_to_ball_force_times_radius) + for (int k = 0; k < ParallelUtilities::GetNumThreads(); k++){ + + for (ElementsArrayType::iterator it = GetElementPartitionBegin(r_model_part, k); it != GetElementPartitionEnd(r_model_part, k); ++it){ + if ((it)->IsNot(DEMFlags::BELONGS_TO_A_CLUSTER)) { + double particle_normal_ball_to_ball_force_times_radius = 0.0; + + (it)->Calculate(PARTICLE_MAX_NORMAL_BALL_TO_BALL_FORCE_TIMES_RADIUS, particle_normal_ball_to_ball_force_times_radius, r_model_part.GetProcessInfo()); + + if (particle_max_normal_ball_to_ball_force_times_radius < particle_normal_ball_to_ball_force_times_radius){ + particle_max_normal_ball_to_ball_force_times_radius = particle_normal_ball_to_ball_force_times_radius; + } + } + } + + } + + return particle_max_normal_ball_to_ball_force_times_radius * r_model_part.GetCommunicator().LocalMesh().Elements().size(); + } + //*************************************************************************************************************** //*************************************************************************************************************** diff --git a/applications/DEMApplication/python_scripts/DEM_analysis_stage.py b/applications/DEMApplication/python_scripts/DEM_analysis_stage.py index 5d7686762e72..792421cbf22b 100644 --- a/applications/DEMApplication/python_scripts/DEM_analysis_stage.py +++ b/applications/DEMApplication/python_scripts/DEM_analysis_stage.py @@ -302,6 +302,13 @@ def Initialize(self): if self.DEM_parameters["output_configuration"]["print_number_of_neighbours_histogram"].GetBool(): self.PreUtilities.PrintNumberOfNeighboursHistogram(self.spheres_model_part, os.path.join(self.graphs_path, "number_of_neighbours_histogram.txt")) + self.BoundingBoxMinX_update = self.DEM_parameters["BoundingBoxMinX"].GetDouble() + self.BoundingBoxMinY_update = self.DEM_parameters["BoundingBoxMinY"].GetDouble() + self.BoundingBoxMinZ_update = self.DEM_parameters["BoundingBoxMinZ"].GetDouble() + self.BoundingBoxMaxX_update = self.DEM_parameters["BoundingBoxMaxX"].GetDouble() + self.BoundingBoxMaxY_update = self.DEM_parameters["BoundingBoxMaxY"].GetDouble() + self.BoundingBoxMaxZ_update = self.DEM_parameters["BoundingBoxMaxZ"].GetDouble() + def SetMaterials(self): self.ReadMaterialsFile() @@ -531,21 +538,22 @@ def InitializeSolutionStep(self): def UpdateSearchStartegyAndCPlusPlusStrategy(self): + delta_time = self.spheres_model_part.ProcessInfo.GetValue(DELTA_TIME) move_velocity = self.DEM_parameters["BoundingBoxMoveVelocity"].GetDouble() - BoundingBoxMinX_update = self.DEM_parameters["BoundingBoxMinX"].GetDouble() + self.time * move_velocity - BoundingBoxMinY_update = self.DEM_parameters["BoundingBoxMinY"].GetDouble() + self.time * move_velocity - BoundingBoxMinZ_update = self.DEM_parameters["BoundingBoxMinZ"].GetDouble() + self.time * move_velocity - BoundingBoxMaxX_update = self.DEM_parameters["BoundingBoxMaxX"].GetDouble() - self.time * move_velocity - BoundingBoxMaxY_update = self.DEM_parameters["BoundingBoxMaxY"].GetDouble() - self.time * move_velocity - BoundingBoxMaxZ_update = self.DEM_parameters["BoundingBoxMaxZ"].GetDouble() - self.time * move_velocity - - self._GetSolver().search_strategy = OMP_DEMSearch(BoundingBoxMinX_update, - BoundingBoxMinY_update, - BoundingBoxMinZ_update, - BoundingBoxMaxX_update, - BoundingBoxMaxY_update, - BoundingBoxMaxZ_update) + self.BoundingBoxMinX_update += delta_time * move_velocity + self.BoundingBoxMinY_update += delta_time * move_velocity + self.BoundingBoxMinZ_update += delta_time * move_velocity + self.BoundingBoxMaxX_update -= delta_time * move_velocity + self.BoundingBoxMaxY_update -= delta_time * move_velocity + self.BoundingBoxMaxZ_update -= delta_time * move_velocity + + self._GetSolver().search_strategy = OMP_DEMSearch(self.BoundingBoxMinX_update, + self.BoundingBoxMinY_update, + self.BoundingBoxMinZ_update, + self.BoundingBoxMaxX_update, + self.BoundingBoxMaxY_update, + self.BoundingBoxMaxZ_update) self._GetSolver().UpdateCPlusPlusStrategy() def UpdateIsTimeToPrintInModelParts(self, is_time_to_print): @@ -733,12 +741,12 @@ def MeasureSphereForGettingPackingProperties(self, radius, center_x, center_y, c ''' if type == "porosity": - measure_sphere_volume = 4/3 * math.pi * radius * radius * radius + measure_sphere_volume = 4.0 / 3.0 * math.pi * radius * radius * radius sphere_volume_inside_range = 0.0 measured_porosity = 0.0 for node in self.spheres_model_part.Nodes: - + r = node.GetSolutionStepValue(RADIUS) x = node.X y = node.Y @@ -750,7 +758,7 @@ def MeasureSphereForGettingPackingProperties(self, radius, center_x, center_y, c sphere_volume_inside_range += 4/3 * math.pi * r * r * r - elif center_to_sphere_distance < (radius + r): + elif center_to_sphere_distance <= (radius + r): other_part_d = radius - (radius * radius + center_to_sphere_distance * center_to_sphere_distance - r * r) / (center_to_sphere_distance * 2) @@ -759,8 +767,8 @@ def MeasureSphereForGettingPackingProperties(self, radius, center_x, center_y, c cross_volume = math.pi * other_part_d * other_part_d * (radius - 1/3 * other_part_d) + math.pi * my_part_d * my_part_d * (r - 1/3 * my_part_d) sphere_volume_inside_range += cross_volume - - measured_porosity = 1 - (sphere_volume_inside_range / measure_sphere_volume) + + measured_porosity = 1.0 - (sphere_volume_inside_range / measure_sphere_volume) return measured_porosity diff --git a/applications/DEMApplication/python_scripts/DEM_procedures.py b/applications/DEMApplication/python_scripts/DEM_procedures.py index 7b3702b570f7..10d22c219aae 100644 --- a/applications/DEMApplication/python_scripts/DEM_procedures.py +++ b/applications/DEMApplication/python_scripts/DEM_procedures.py @@ -306,6 +306,17 @@ def CalculateEnergyAndPlot(self, time): self.energy_graph_counter += 1 + def CalculateNormalizedKinematicEnergy(self): + + particle_number_times_max_normal_ball_to_ball_force_times_radius = self.SpheresEnergyUtil.CalculateParticleNumberTimesMaxNormalBallToBallForceTimesRadius(self.SpheresModelPart) + + if particle_number_times_max_normal_ball_to_ball_force_times_radius == 0.0: + normalized_kinematic_energy = 0.0 + else: + normalized_kinematic_energy = self.kinematic_energy / particle_number_times_max_normal_ball_to_ball_force_times_radius + + return normalized_kinematic_energy + def CalculateEnergy(self): self.translational_kinematic_energy = self.SpheresEnergyUtil.CalculateTranslationalKinematicEnergy(self.SpheresModelPart) + self.ClusterEnergyUtil.CalculateTranslationalKinematicEnergy(self.ClusterModelPart) diff --git a/applications/DEMApplication/python_scripts/dem_default_input_parameters.py b/applications/DEMApplication/python_scripts/dem_default_input_parameters.py index d8ee80e3b625..ce32f1543f12 100644 --- a/applications/DEMApplication/python_scripts/dem_default_input_parameters.py +++ b/applications/DEMApplication/python_scripts/dem_default_input_parameters.py @@ -39,7 +39,11 @@ def GetDefaultInputParameters(): "VelocityTrapMinZ" : 0.0, "RotationOption" : true, "CleanIndentationsOption" : false, + "CleanIndentationsV2Option" : false, "RemoveBallsInEmbeddedOption" : false, + "RadiusExpansionOption" : false, + "RadiusExpansionRate" : 0.0, + "RadiusMultiplierMax" : 0.0, "solver_settings" : {}, "creator_destructor_settings" : {}, "echo_level" : 1, diff --git a/applications/DEMApplication/python_scripts/sphere_strategy.py b/applications/DEMApplication/python_scripts/sphere_strategy.py index fb0d46fa777a..63079b748436 100644 --- a/applications/DEMApplication/python_scripts/sphere_strategy.py +++ b/applications/DEMApplication/python_scripts/sphere_strategy.py @@ -66,11 +66,17 @@ def __init__(self, all_model_parts, creator_destructor, dem_fem_search, DEM_para #self.time_integration_scheme.SetRotationOption(self.rotation_option) self.clean_init_indentation_option = DEM_parameters["CleanIndentationsOption"].GetBool() + self.clean_init_indentation_v2_option = DEM_parameters["CleanIndentationsV2Option"].GetBool() if self.clean_init_indentation_option and self._GetInputType() == 'rest': Logger.PrintWarning("DEM", '\nWARNING!: \'clean_indentations_option\' is set to true in a restarted simulation. The particles\' radii could be modified before the first time step.\n' * 50) - self.contact_mesh_option = 0 + if self.clean_init_indentation_option and self.clean_init_indentation_v2_option: + Logger.PrintWarning("DEM", '\nWARNING!: "CleanIndentationsOption" and "CleanIndentationsV2Option" can not both be "True" at the same time.\n') + from sys import exit + exit(0) + + self.contact_mesh_option = 0 if "ContactMeshOption" in DEM_parameters.keys(): self.contact_mesh_option = DEM_parameters["ContactMeshOption"].GetBool() self.automatic_bounding_box_option = DEM_parameters["AutomaticBoundingBoxOption"].GetBool() @@ -177,6 +183,18 @@ def __init__(self, all_model_parts, creator_destructor, dem_fem_search, DEM_para else: self.global_viscous_damping = DEM_parameters["GlobalViscousDamping"].GetDouble() + if "RadiusExpansionOption" in DEM_parameters.keys(): + self.radius_expansion_option = DEM_parameters["RadiusExpansionOption"].GetBool() + + if "RadiusExpansionRate" in DEM_parameters.keys(): + self.radius_expansion_rate = DEM_parameters["RadiusExpansionRate"].GetDouble() + + if "RadiusMultiplierMax" in DEM_parameters.keys(): + self.radius_multiplier_max = DEM_parameters["RadiusMultiplierMax"].GetDouble() + + if "EnergyCalculationOption" in DEM_parameters.keys(): + self.energy_calculation_option = DEM_parameters["EnergyCalculationOption"].GetBool() + # PRINTING VARIABLES self.print_export_id = DEM_parameters["PostExportId"].GetBool() self.print_export_skin_sphere = 0 @@ -251,6 +269,7 @@ def SetVariablesAndOptions(self): self.spheres_model_part.ProcessInfo.SetValue(FIX_VELOCITIES_FLAG, self.fix_velocities_flag) self.spheres_model_part.ProcessInfo.SetValue(NEIGH_INITIALIZED, 0) self.spheres_model_part.ProcessInfo.SetValue(CLEAN_INDENT_OPTION, self.clean_init_indentation_option) + self.spheres_model_part.ProcessInfo.SetValue(CLEAN_INDENT_V2_OPTION, self.clean_init_indentation_v2_option) self.spheres_model_part.ProcessInfo.SetValue(BOUNDING_BOX_START_TIME, self.bounding_box_start_time) self.spheres_model_part.ProcessInfo.SetValue(BOUNDING_BOX_STOP_TIME, self.bounding_box_stop_time) self.spheres_model_part.ProcessInfo.SetValue(COMPUTE_STRESS_TENSOR_OPTION, self.compute_stress_tensor_option) @@ -273,6 +292,13 @@ def SetVariablesAndOptions(self): self.spheres_model_part.ProcessInfo.SetValue(GLOBAL_DAMPING, self.global_damping) self.spheres_model_part.ProcessInfo.SetValue(GLOBAL_VISCOUS_DAMPING, self.global_viscous_damping) + #Radius expansion method + self.spheres_model_part.ProcessInfo.SetValue(IS_RADIUS_EXPANSION, self.radius_expansion_option) + self.spheres_model_part.ProcessInfo.SetValue(RADIUS_EXPANSION_RATE, self.radius_expansion_rate) + self.spheres_model_part.ProcessInfo.SetValue(RADIUS_MULTIPLIER_MAX, self.radius_multiplier_max) + + self.spheres_model_part.ProcessInfo.SetValue(ENERGY_CALCULATION_OPTION, self.energy_calculation_option) + # SEARCH-RELATED self.spheres_model_part.ProcessInfo.SetValue(SEARCH_RADIUS_INCREMENT, self.search_increment) self.spheres_model_part.ProcessInfo.SetValue(SEARCH_RADIUS_INCREMENT_FOR_WALLS, self.search_increment_for_walls) @@ -456,7 +482,7 @@ def InitializeSolutionStep(self): self.cplusplus_strategy.InitializeSolutionStep() def SetNormalRadiiOnAllParticles(self): - self.cplusplus_strategy.SetNormalRadiiOnAllParticles(self.spheres_model_part) + self.cplusplus_strategy.SetNormalRadiiOnAllParticlesBeforeInitilization(self.spheres_model_part) def SetSearchRadiiOnAllParticles(self): diff --git a/applications/DEMApplication/python_scripts/verlet_continuum_sphere_strategy.py b/applications/DEMApplication/python_scripts/verlet_continuum_sphere_strategy.py index 31f7fcb3ff21..7174edae67a1 100644 --- a/applications/DEMApplication/python_scripts/verlet_continuum_sphere_strategy.py +++ b/applications/DEMApplication/python_scripts/verlet_continuum_sphere_strategy.py @@ -34,6 +34,7 @@ def CreateCPlusPlusStrategy(self): self.model_part.ProcessInfo.SetValue(FIX_VELOCITIES_FLAG, self.fix_velocities_flag) self.model_part.ProcessInfo.SetValue(NEIGH_INITIALIZED, 0) self.model_part.ProcessInfo.SetValue(CLEAN_INDENT_OPTION, self.clean_init_indentation_option) + self.model_part.ProcessInfo.SetValue(CLEAN_INDENT_V2_OPTION, self.clean_init_indentation_v2_option) self.model_part.ProcessInfo.SetValue(COMPUTE_STRESS_TENSOR_OPTION, self.compute_stress_tensor_option) # GLOBAL PHISICAL ASPECTS diff --git a/applications/DEMApplication/tests/test_DEM_3D_continuum.py b/applications/DEMApplication/tests/test_DEM_3D_continuum.py index 5d3b98a66aec..10f90d47a2f6 100644 --- a/applications/DEMApplication/tests/test_DEM_3D_continuum.py +++ b/applications/DEMApplication/tests/test_DEM_3D_continuum.py @@ -28,16 +28,16 @@ def GetProblemNameWithPath(self): def CheckValues(self, x_vel, z_vel, x_force, z_force, dem_pressure, z_elastic, shear, x_tangential): tol = 1.0e-8 # DEM reference values - x_vel_ref = 0.028907825348927448 - z_vel_ref = -0.8757276957864403 - x_force_ref = -26919.437972831598 - z_force_ref = 1970950.3578554934 + x_vel_ref = 0.02677117811730761 + z_vel_ref = -0.8785947012104277 + x_force_ref = -28004.11984689704 + z_force_ref = 1976786.9183720595 #FEM reference values - dem_pressure_ref = 21566.85065708402 - z_elastic_ref = -273575.41245014494 - shear_ref = 362.391011482587 - x_tangential_ref = 6039.850191376444 + dem_pressure_ref = 21615.368565335117 + z_elastic_ref = -274076.08804810344 + shear_ref = 368.8978792516498 + x_tangential_ref = 6148.297987527496 self.assertAlmostEqual(x_vel, x_vel_ref, delta=tol) self.assertAlmostEqual(z_vel, z_vel_ref, delta=tol) diff --git a/applications/DEMApplication/tests/test_DEM_3D_continuum_vs_discontinuum.py b/applications/DEMApplication/tests/test_DEM_3D_continuum_vs_discontinuum.py index f60810b6d48e..ab57001ebcc7 100644 --- a/applications/DEMApplication/tests/test_DEM_3D_continuum_vs_discontinuum.py +++ b/applications/DEMApplication/tests/test_DEM_3D_continuum_vs_discontinuum.py @@ -16,6 +16,10 @@ def GetFilePath(fileName): class DEM3D_ContinuumTestVsDiscontinuumSolution(KratosMultiphysics.DEMApplication.DEM_analysis_stage.DEMAnalysisStage, KratosUnittest.TestCase): + def Initialize(self): + super().Initialize() + self._GetSolver().cplusplus_strategy.BreakAllBonds() + def InitializeSolutionStep(self): super().InitializeSolutionStep() diff --git a/applications/DEMApplication/tests/test_DEM_schemes.py b/applications/DEMApplication/tests/test_DEM_schemes.py index f31bb9cc80cb..fdfca4602e9a 100644 --- a/applications/DEMApplication/tests/test_DEM_schemes.py +++ b/applications/DEMApplication/tests/test_DEM_schemes.py @@ -36,10 +36,10 @@ def GetProblemNameWithPath(self): def CheckValues(self, x_vel, dem_pressure): tol = 1.0e-10 - x_vel_ref = 0.028907825348927448 + x_vel_ref = 0.02677117811730761 self.assertAlmostEqual(x_vel, x_vel_ref, delta=tol) - dem_pressure_ref = 21566.85065708402 + dem_pressure_ref = 21615.368565335117 self.assertAlmostEqual(dem_pressure, dem_pressure_ref, delta=tol) def Finalize(self): @@ -105,29 +105,29 @@ class DEM3D_TaylorTestSolution(DEM3D_ForwardEulerTestSolution): def CheckValues(self, x_vel, dem_pressure): tol = 1.0e-10 - x_vel_ref = 0.028709756132288513 + x_vel_ref = 0.026581603224742834 self.assertAlmostEqual(x_vel, x_vel_ref, delta=tol) - dem_pressure_ref = 21550.45232404601 + dem_pressure_ref = 21598.828562997718 self.assertAlmostEqual(dem_pressure, dem_pressure_ref, delta=tol) class DEM3D_SymplecticTestSolution(DEM3D_ForwardEulerTestSolution): def CheckValues(self, x_vel, dem_pressure): tol = 1.0e-10 - x_vel_ref = 0.028515905722678703 + x_vel_ref = 0.026396099600870276 self.assertAlmostEqual(x_vel, x_vel_ref, delta=tol) - dem_pressure_ref = 21534.129347072263 + dem_pressure_ref = 21582.363307819433 self.assertAlmostEqual(dem_pressure, dem_pressure_ref, delta=tol) class DEM3D_VerletTestSolution(DEM3D_ForwardEulerTestSolution): def CheckValues(self,x_vel, dem_pressure): tol = 1.0e-10 - x_vel_ref = 0.028603162986524718 + x_vel_ref = 0.026480762680319958 self.assertAlmostEqual(x_vel, x_vel_ref, delta=tol) - dem_pressure_ref = 21560.85829962128 + dem_pressure_ref = 21609.179544197963 self.assertAlmostEqual(dem_pressure, dem_pressure_ref, delta=tol) class TestDEMSchemes(KratosUnittest.TestCase): diff --git a/applications/DEMApplication/tests/test_dem_3d_parallel_bond_model.py b/applications/DEMApplication/tests/test_dem_3d_parallel_bond_model.py index b7303af84466..9f60be9f0a37 100644 --- a/applications/DEMApplication/tests/test_dem_3d_parallel_bond_model.py +++ b/applications/DEMApplication/tests/test_dem_3d_parallel_bond_model.py @@ -23,21 +23,21 @@ def GetProblemNameWithPath(self): def FinalizeSolutionStep(self): super().FinalizeSolutionStep() - tolerance = 1e-7 + tolerance = 1e-6 for node in self.spheres_model_part.Nodes: velocity = node.GetSolutionStepValue(Kratos.VELOCITY) angular_velocity = node.GetSolutionStepValue(Kratos.ANGULAR_VELOCITY) if node.Id == 1: if self.time > 0.0098 and self.time < 0.00981: - expected_value = -3.365862016784862e-05 + expected_value = -3.387935641635233e-05 self.CheckValues(velocity, 0, expected_value, tolerance) expected_value = 0.0 self.CheckValues(velocity, 1, expected_value, tolerance) - expected_value = 3.345233292428133e-05 + expected_value = 3.20846623155654e-05 self.CheckValues(velocity, 2, expected_value, tolerance) expected_value = 0.0 self.CheckValues(angular_velocity, 0, expected_value, tolerance) - expected_value = 4.011134246117183e-05 + expected_value = 7.160580099884259e-05 self.CheckValues(angular_velocity, 1, expected_value, tolerance) expected_value = 0.0 self.CheckValues(angular_velocity, 2, expected_value, tolerance) diff --git a/applications/DEMApplication/tests/test_dem_3d_smooth_joint_model.py b/applications/DEMApplication/tests/test_dem_3d_smooth_joint_model.py index d1130d88ad9b..e4c488d6e0a9 100644 --- a/applications/DEMApplication/tests/test_dem_3d_smooth_joint_model.py +++ b/applications/DEMApplication/tests/test_dem_3d_smooth_joint_model.py @@ -29,7 +29,7 @@ def FinalizeSolutionStep(self): angular_velocity = node.GetSolutionStepValue(Kratos.ANGULAR_VELOCITY) if node.Id == 1: if self.time > 0.00098 and self.time < 0.000982: - expected_value = -0.0007822842347468785 + expected_value = 0.0 self.CheckValues(velocity, 0, expected_value, tolerance) expected_value = 0.0 self.CheckValues(velocity, 1, expected_value, tolerance)