diff --git a/applications/DEMApplication/CMakeLists.txt b/applications/DEMApplication/CMakeLists.txt index ac0ce65062ee..bd2b642978e6 100755 --- a/applications/DEMApplication/CMakeLists.txt +++ b/applications/DEMApplication/CMakeLists.txt @@ -68,6 +68,7 @@ set(KRATOS_DEM_APPLICATION_CORE ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_KDEM_CamClay_CL.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/dem_kdem_fissured_rock_cl.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_ExponentialHC_CL.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_D_Linear_Simple_Coulomb_CL.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_D_Linear_viscous_Coulomb_CL.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_D_Linear_HighStiffness_CL.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_D_Linear_HighStiffness_2D_CL.cpp diff --git a/applications/DEMApplication/DEM_application.cpp b/applications/DEMApplication/DEM_application.cpp index 6a05ec219e23..e68c916d7657 100644 --- a/applications/DEMApplication/DEM_application.cpp +++ b/applications/DEMApplication/DEM_application.cpp @@ -25,6 +25,7 @@ #include "custom_constitutive/DEM_Dempack_dev_CL.h" #include "custom_constitutive/dem_kdem_2d_cl.h" #include "custom_constitutive/dem_kdem_fabric_2d_cl.h" +#include "custom_constitutive/DEM_D_Linear_Simple_Coulomb_CL.h" #include "custom_constitutive/DEM_D_Linear_viscous_Coulomb_CL.h" #include "custom_constitutive/DEM_D_Linear_viscous_Coulomb_2D_CL.h" #include "custom_constitutive/DEM_D_Hertz_viscous_Coulomb_CL.h" @@ -1034,28 +1035,21 @@ void KratosDEMApplication::Register() { KRATOS_REGISTER_CONDITION("RigidEdge2D2N", mRigidEdge2D2N) KRATOS_REGISTER_CONDITION("RigidEdge2D1N", mRigidEdge2D1N) - // SERIALIZER Serializer::Register("PropertiesProxy", PropertiesProxy()); - Serializer::Register( - "DEM_D_Linear_viscous_Coulomb", DEM_D_Linear_viscous_Coulomb()); - Serializer::Register( - "DEM_D_Linear_viscous_Coulomb2D", DEM_D_Linear_viscous_Coulomb2D()); - Serializer::Register( - "DEM_D_Hertz_viscous_Coulomb", DEM_D_Hertz_viscous_Coulomb()); - Serializer::Register( - "DEM_D_Hertz_viscous_Coulomb2D", DEM_D_Hertz_viscous_Coulomb2D()); + Serializer::Register("DEM_D_Linear_Simple_Coulomb", DEM_D_Linear_Simple_Coulomb()); + Serializer::Register("DEM_D_Linear_viscous_Coulomb", DEM_D_Linear_viscous_Coulomb()); + Serializer::Register("DEM_D_Linear_viscous_Coulomb2D", DEM_D_Linear_viscous_Coulomb2D()); + Serializer::Register("DEM_D_Hertz_viscous_Coulomb", DEM_D_Hertz_viscous_Coulomb()); + Serializer::Register("DEM_D_Hertz_viscous_Coulomb2D", DEM_D_Hertz_viscous_Coulomb2D()); Serializer::Register("DEM_D_JKR_Cohesive_Law", DEM_D_JKR_Cohesive_Law()); Serializer::Register("DEM_D_Bentonite_Colloid", DEM_D_Bentonite_Colloid()); Serializer::Register("DEM_D_DMT_Cohesive_Law", DEM_D_DMT_Cohesive_Law()); Serializer::Register("DEM_D_Stress_Dependent_Cohesive", DEM_D_Stress_Dependent_Cohesive()); - Serializer::Register( - "DEM_D_Linear_Custom_Constants", DEM_D_Linear_Custom_Constants()); - Serializer::Register( - "DEM_D_Conical_damage", DEM_D_Conical_damage()); - Serializer::Register("DEM_D_Hertz_viscous_Coulomb_Nestle", - DEM_D_Hertz_viscous_Coulomb_Nestle()); + Serializer::Register("DEM_D_Linear_Custom_Constants", DEM_D_Linear_Custom_Constants()); + Serializer::Register("DEM_D_Conical_damage", DEM_D_Conical_damage()); + Serializer::Register("DEM_D_Hertz_viscous_Coulomb_Nestle", DEM_D_Hertz_viscous_Coulomb_Nestle()); Serializer::Register("DEM_D_Quadratic", DEM_D_Quadratic()); Serializer::Register("DEM_D_Linear_classic", DEM_D_Linear_classic()); Serializer::Register("DEM_D_void", DEM_D_void()); diff --git a/applications/DEMApplication/custom_constitutive/DEM_D_Linear_Simple_Coulomb_CL.cpp b/applications/DEMApplication/custom_constitutive/DEM_D_Linear_Simple_Coulomb_CL.cpp new file mode 100644 index 000000000000..e43a59922bb5 --- /dev/null +++ b/applications/DEMApplication/custom_constitutive/DEM_D_Linear_Simple_Coulomb_CL.cpp @@ -0,0 +1,103 @@ +// Kratos Multi-Physics - DEM Application +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Rafael Rangel (rrangel@cimne.upc.edu) +// +// References using this model: +// R.L. Rangel et al. (2024). Multiscale data-driven modeling of the thermomechanical behavior of granular media with thermal expansion effects. Computers and Geotechnics, 176:106789. +// R.L. Rangel et al. (2024). A continuum-discrete multiscale methodology using machine learning for thermal analysis of granular media. Computers and Geotechnics, 168:106118. +// Guo & Zhao (2014). coupled FEM/DEM approach for hierarchical multiscale modelling of granular media. IJNME, 99(11):789-818. +// + +#include "DEM_D_Linear_Simple_Coulomb_CL.h" +#include "custom_elements/spheric_particle.h" + +namespace Kratos { + + DEMDiscontinuumConstitutiveLaw::Pointer DEM_D_Linear_Simple_Coulomb::Clone() const { + DEMDiscontinuumConstitutiveLaw::Pointer p_clone(new DEM_D_Linear_Simple_Coulomb(*this)); + return p_clone; + } + + std::unique_ptr DEM_D_Linear_Simple_Coulomb::CloneUnique() { + return Kratos::make_unique(); + } + + std::string DEM_D_Linear_Simple_Coulomb::GetTypeOfLaw() { + std::string type_of_law = "Linear"; + return type_of_law; + } + + void DEM_D_Linear_Simple_Coulomb::Check(Properties::Pointer pProp) const { + if (!pProp->Has(STATIC_FRICTION)) { + if (!pProp->Has(FRICTION)) { //deprecated since April 6th, 2020 + KRATOS_WARNING("DEM") << std::endl; + KRATOS_WARNING("DEM") << "WARNING: Variable STATIC_FRICTION or FRICTION should be present in the properties when using DEMDiscontinuumConstitutiveLaw. 0.0 value assigned by default." << std::endl; + KRATOS_WARNING("DEM") << std::endl; + pProp->GetValue(STATIC_FRICTION) = 0.0; + } + else { + pProp->GetValue(STATIC_FRICTION) = pProp->GetValue(FRICTION); + } + } + } + + void DEM_D_Linear_Simple_Coulomb::CalculateForces(const ProcessInfo& r_process_info, + const double OldLocalElasticContactForce[3], + double LocalElasticContactForce[3], + double LocalDeltDisp[3], + double LocalRelVel[3], + double indentation, + double previous_indentation, + double ViscoDampingLocalContactForce[3], + double& cohesive_force, + SphericParticle* element1, + SphericParticle* element2, + bool& sliding, + double LocalCoordSystem[3][3]) { + // Compute stiffness coefficients + const double my_radius = element1->GetRadius(); + const double other_radius = element2->GetRadius(); + const double equiv_radius = 2.0 * my_radius * other_radius / (my_radius + other_radius); + const double my_young = element1->GetYoung(); + const double my_poisson = element1->GetPoisson(); + + mKn = my_young * equiv_radius; // normal + mKt = my_poisson * mKn; // tangent + + // Compute normal and tangent forces + const double normal_contact_force = mKn * indentation; + LocalElasticContactForce[2] = normal_contact_force; + + CalculateTangentialForceWithNeighbour(normal_contact_force, OldLocalElasticContactForce, LocalElasticContactForce, LocalDeltDisp, sliding, element1, element2); + } + + void DEM_D_Linear_Simple_Coulomb::CalculateForcesWithFEM(const ProcessInfo& r_process_info, + const double OldLocalElasticContactForce[3], + double LocalElasticContactForce[3], + double LocalDeltDisp[3], + double LocalRelVel[3], + double indentation, + double previous_indentation, + double ViscoDampingLocalContactForce[3], + double& cohesive_force, + SphericParticle* const element, + Condition* const wall, + bool& sliding) { + // Compute stiffness coefficients + const double my_radius = element->GetRadius(); + const double my_young = element->GetYoung(); + const double my_poisson = element->GetPoisson(); + + mKn = my_young * my_radius; // normal + mKt = my_poisson * mKn; // tangent + + // Compute normal and tangent forces + const double normal_contact_force = mKn * indentation; + LocalElasticContactForce[2] = normal_contact_force; + + CalculateTangentialForceWithNeighbour(normal_contact_force, OldLocalElasticContactForce, LocalElasticContactForce, LocalDeltDisp, sliding, element, wall); + } +} // namespace Kratos diff --git a/applications/DEMApplication/custom_constitutive/DEM_D_Linear_Simple_Coulomb_CL.h b/applications/DEMApplication/custom_constitutive/DEM_D_Linear_Simple_Coulomb_CL.h new file mode 100644 index 000000000000..35cae548cb30 --- /dev/null +++ b/applications/DEMApplication/custom_constitutive/DEM_D_Linear_Simple_Coulomb_CL.h @@ -0,0 +1,106 @@ +// Kratos Multi-Physics - DEM Application +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Rafael Rangel (rrangel@cimne.upc.edu) +// +// References using this model: +// R.L. Rangel et al. (2024). Multiscale data-driven modeling of the thermomechanical behavior of granular media with thermal expansion effects. Computers and Geotechnics, 176:106789. +// R.L. Rangel et al. (2024). A continuum-discrete multiscale methodology using machine learning for thermal analysis of granular media. Computers and Geotechnics, 168:106118. +// Guo & Zhao (2014). coupled FEM/DEM approach for hierarchical multiscale modelling of granular media. IJNME, 99(11):789-818. +// + +#pragma once + +#include +#include "DEM_discontinuum_constitutive_law.h" + +namespace Kratos { + + class SphericParticle; + class KRATOS_API(DEM_APPLICATION) DEM_D_Linear_Simple_Coulomb : public DEMDiscontinuumConstitutiveLaw { + + public: + + using DEMDiscontinuumConstitutiveLaw::CalculateNormalForce; + + KRATOS_CLASS_POINTER_DEFINITION(DEM_D_Linear_Simple_Coulomb); + + DEM_D_Linear_Simple_Coulomb() {} + ~DEM_D_Linear_Simple_Coulomb() {} + + DEMDiscontinuumConstitutiveLaw::Pointer Clone() const override; + std::unique_ptr CloneUnique() override; + + std::string GetTypeOfLaw() override; + virtual void Check(Properties::Pointer pProp) const override; + + void CalculateForces(const ProcessInfo& r_process_info, + const double OldLocalElasticContactForce[3], + double LocalElasticContactForce[3], + double LocalDeltDisp[3], + double LocalRelVel[3], + double indentation, + double previous_indentation, + double ViscoDampingLocalContactForce[3], + double& cohesive_force, + SphericParticle* element1, + SphericParticle* element2, + bool& sliding, + double LocalCoordSystem[3][3]) override; + + void CalculateForcesWithFEM(const ProcessInfo& r_process_info, + const double OldLocalElasticContactForce[3], + double LocalElasticContactForce[3], + double LocalDeltDisp[3], + double LocalRelVel[3], + double indentation, + double previous_indentation, + double ViscoDampingLocalContactForce[3], + double& cohesive_force, + SphericParticle* const element, + Condition* const wall, + bool& sliding) override; + + template + void CalculateTangentialForceWithNeighbour(const double normal_contact_force, + const double OldLocalElasticContactForce[3], + double LocalElasticContactForce[3], + const double LocalDeltDisp[3], + bool& sliding, + SphericParticle* const element, + NeighbourClassType* const neighbour) { + // Compute shear force + LocalElasticContactForce[0] = OldLocalElasticContactForce[0] - mKt * LocalDeltDisp[0]; + LocalElasticContactForce[1] = OldLocalElasticContactForce[1] - mKt * LocalDeltDisp[1]; + const double tangent_contact_force = sqrt(LocalElasticContactForce[0] * LocalElasticContactForce[0] + LocalElasticContactForce[1] * LocalElasticContactForce[1]); + + // Compute maximum admisible shear force + Properties& properties_of_this_contact = element->GetProperties().GetSubProperties(neighbour->GetProperties().Id()); + const double friction_angle_tg = std::tan(properties_of_this_contact[STATIC_FRICTION]); + const double MaximumAdmisibleShearForce = normal_contact_force * friction_angle_tg; + + // Check for sliding: apply Coulomb friction condition + if (tangent_contact_force > MaximumAdmisibleShearForce) { + sliding = true; + const double fraction = MaximumAdmisibleShearForce / tangent_contact_force; + LocalElasticContactForce[0] *= fraction; + LocalElasticContactForce[1] *= fraction; + } + } + + private: + + friend class Serializer; + + virtual void save(Serializer& rSerializer) const override { + KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, DEMDiscontinuumConstitutiveLaw) + } + + virtual void load(Serializer& rSerializer) override { + KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, DEMDiscontinuumConstitutiveLaw) + } + + }; +} 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 a18abf488aeb..5567e78773b6 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 @@ -14,6 +14,7 @@ #include "custom_constitutive/DEM_compound_constitutive_law.h" #include "custom_constitutive/DEM_compound_constitutive_law_for_PBM.h" +#include "custom_constitutive/DEM_D_Linear_Simple_Coulomb_CL.h" #include "custom_constitutive/DEM_D_Linear_viscous_Coulomb_CL.h" #include "custom_constitutive/DEM_D_Hertz_viscous_Coulomb_CL.h" #include "custom_constitutive/DEM_D_Hertz_viscous_Coulomb_Nestle_CL.h" @@ -97,6 +98,10 @@ void AddCustomConstitutiveLawsToPython(pybind11::module& m) { .def(py::init<>()) ; + py::class_(m, "DEM_D_Linear_Simple_Coulomb") + .def(py::init<>()) + ; + py::class_(m, "DEM_D_Linear_viscous_Coulomb2D") .def(py::init<>()) ;