[DEMApplication] Fixed some bugs and added a smooth-joint model

applications/DEMApplication/CMakeLists.txt

@@ -84,10 +84,12 @@ set(KRATOS_DEM_APPLICATION_CORE
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/dem_d_linear_custom_constants_cl.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_D_Conical_damage_CL.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_D_Quadratic_CL.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_D_void_CL.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/dem_kdem_2d_cl.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/dem_kdem_fabric_2d_cl.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_D_Hertz_viscous_Coulomb_Nestle_CL.cpp
     ${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_utilities/inlet.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/force_based_inlet.cpp

applications/DEMApplication/DEM_application.cpp

@@ -36,6 +36,7 @@
 #include "custom_constitutive/dem_d_linear_custom_constants_cl.h"
 #include "custom_constitutive/DEM_D_Conical_damage_CL.h"
 #include "custom_constitutive/DEM_D_Quadratic_CL.h"
+#include "custom_constitutive/DEM_D_void_CL.h"
 #include "custom_constitutive/DEM_D_Linear_classic_CL.h"
 #include "custom_constitutive/DEM_KDEM_fabric_CL.h"
 #include "custom_constitutive/DEM_beam_constitutive_law.h"
@@ -46,6 +47,7 @@
 #include "custom_constitutive/DEM_compound_constitutive_law.h"
 #include "custom_constitutive/DEM_compound_constitutive_law_for_PBM.h"
 #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_rolling_friction_model.h"
 #include "custom_constitutive/DEM_rolling_friction_model_constant_torque.h"
@@ -266,6 +268,12 @@ KRATOS_CREATE_VARIABLE(double, BOND_INTERNAL_FRICC)
 KRATOS_CREATE_VARIABLE(double, BOND_ROTATIONAL_MOMENT_COEFFICIENT_NORMAL)
 KRATOS_CREATE_VARIABLE(double, BOND_ROTATIONAL_MOMENT_COEFFICIENT_TANGENTIAL)
 KRATOS_CREATE_VARIABLE(double, BOND_RADIUS_FACTOR)
+KRATOS_CREATE_VARIABLE(double, JOINT_NORMAL_STIFFNESS)
+KRATOS_CREATE_VARIABLE(double, JOINT_TANGENTIAL_STIFFNESS)
+KRATOS_CREATE_VARIABLE(double, JOINT_NORMAL_DIRECTION_X)
+KRATOS_CREATE_VARIABLE(double, JOINT_NORMAL_DIRECTION_Y)
+KRATOS_CREATE_VARIABLE(double, JOINT_NORMAL_DIRECTION_Z)
+KRATOS_CREATE_VARIABLE(double, JOINT_FRICTION_COEFF)
 KRATOS_CREATE_VARIABLE(double, K_ALPHA) // for DEM_D_Quadratic_LAW
 
 // *************** Continuum only END *************
@@ -275,6 +283,7 @@ KRATOS_CREATE_VARIABLE(double, K_ALPHA) // for DEM_D_Quadratic_LAW
 KRATOS_CREATE_VARIABLE(double, LOCAL_CONTACT_AREA_HIGH)
 KRATOS_CREATE_VARIABLE(double, LOCAL_CONTACT_AREA_LOW)
 KRATOS_CREATE_VARIABLE(double, MEAN_CONTACT_AREA)
+KRATOS_CREATE_VARIABLE(double, CONTACT_RADIUS)
 KRATOS_CREATE_VARIABLE(double, REPRESENTATIVE_VOLUME)
 KRATOS_CREATE_VARIABLE(DenseVector<int>, NEIGHBOUR_IDS)
 KRATOS_CREATE_VARIABLE(DenseVector<double>, NEIGHBOURS_CONTACT_AREAS)
@@ -340,6 +349,7 @@ KRATOS_CREATE_VARIABLE(double, PARTICLE_INELASTIC_FRICTIONAL_ENERGY)
 KRATOS_CREATE_VARIABLE(double, PARTICLE_INELASTIC_ROLLING_RESISTANCE_ENERGY)
 KRATOS_CREATE_VARIABLE(int, COMPUTE_ENERGY_OPTION)
 KRATOS_CREATE_VARIABLE(double, GLOBAL_DAMPING)
+KRATOS_CREATE_VARIABLE(double, GLOBAL_VISCOUS_DAMPING)
 KRATOS_CREATE_VARIABLE(double, NORMAL_IMPACT_VELOCITY)
 KRATOS_CREATE_VARIABLE(double, TANGENTIAL_IMPACT_VELOCITY)
 KRATOS_CREATE_VARIABLE(double, FACE_NORMAL_IMPACT_VELOCITY)
@@ -734,6 +744,12 @@ void KratosDEMApplication::Register() {
     KRATOS_REGISTER_VARIABLE(BOND_ROTATIONAL_MOMENT_COEFFICIENT_NORMAL)
     KRATOS_REGISTER_VARIABLE(BOND_ROTATIONAL_MOMENT_COEFFICIENT_TANGENTIAL)
     KRATOS_REGISTER_VARIABLE(BOND_RADIUS_FACTOR)
+    KRATOS_REGISTER_VARIABLE(JOINT_NORMAL_STIFFNESS)
+    KRATOS_REGISTER_VARIABLE(JOINT_TANGENTIAL_STIFFNESS)
+    KRATOS_REGISTER_VARIABLE(JOINT_NORMAL_DIRECTION_X)
+    KRATOS_REGISTER_VARIABLE(JOINT_NORMAL_DIRECTION_Y)
+    KRATOS_REGISTER_VARIABLE(JOINT_NORMAL_DIRECTION_Z)
+    KRATOS_REGISTER_VARIABLE(JOINT_FRICTION_COEFF)
 
     // *************** Continuum only END *************
 
@@ -743,6 +759,7 @@ void KratosDEMApplication::Register() {
     KRATOS_REGISTER_VARIABLE(LOCAL_CONTACT_AREA_HIGH)
     KRATOS_REGISTER_VARIABLE(LOCAL_CONTACT_AREA_LOW)
     KRATOS_REGISTER_VARIABLE(MEAN_CONTACT_AREA)
+    KRATOS_REGISTER_VARIABLE(CONTACT_RADIUS)
     KRATOS_REGISTER_VARIABLE(REPRESENTATIVE_VOLUME)
     KRATOS_REGISTER_VARIABLE(NEIGHBOUR_IDS)
     KRATOS_REGISTER_VARIABLE(NEIGHBOURS_CONTACT_AREAS)
@@ -815,6 +832,7 @@ void KratosDEMApplication::Register() {
     KRATOS_REGISTER_VARIABLE(PARTICLE_INELASTIC_ROLLING_RESISTANCE_ENERGY)
     KRATOS_REGISTER_VARIABLE(COMPUTE_ENERGY_OPTION)
     KRATOS_REGISTER_VARIABLE(GLOBAL_DAMPING)
+    KRATOS_REGISTER_VARIABLE(GLOBAL_VISCOUS_DAMPING)
     KRATOS_REGISTER_VARIABLE(NORMAL_IMPACT_VELOCITY)
     KRATOS_REGISTER_VARIABLE(TANGENTIAL_IMPACT_VELOCITY)
     KRATOS_REGISTER_VARIABLE(FACE_NORMAL_IMPACT_VELOCITY)
@@ -992,6 +1010,7 @@ void KratosDEMApplication::Register() {
         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());
 
     Serializer::Register("DEM_Dempack", DEM_Dempack());
     Serializer::Register("DEM_Dempack2D", DEM_Dempack2D());
@@ -1005,6 +1024,7 @@ void KratosDEMApplication::Register() {
     Serializer::Register("DEM_KDEM2D", DEM_KDEM2D());
     Serializer::Register("DEM_ExponentialHC", DEM_ExponentialHC());
     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("DEMRollingFrictionModelConstantTorque", DEMRollingFrictionModelConstantTorque());
     Serializer::Register("DEMRollingFrictionModelViscousTorque", DEMRollingFrictionModelViscousTorque());

applications/DEMApplication/DEM_application_variables.h

@@ -236,6 +236,15 @@ namespace Kratos
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, BOND_RADIUS_FACTOR)
   // ****** for parallel bond model END *****
 
+  // ****** for smooth joint model BEGIN ***
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, JOINT_NORMAL_STIFFNESS)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, JOINT_TANGENTIAL_STIFFNESS)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, JOINT_NORMAL_DIRECTION_X)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, JOINT_NORMAL_DIRECTION_Y)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, JOINT_NORMAL_DIRECTION_Z)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, JOINT_FRICTION_COEFF)
+  // ****** for smooth joint model BEGIN ***
+
   // *************** Continuum only END *************
 
   // GEOMETRIC PARAMETERS
@@ -244,6 +253,7 @@ namespace Kratos
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, LOCAL_CONTACT_AREA_HIGH)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, LOCAL_CONTACT_AREA_LOW)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, MEAN_CONTACT_AREA)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, CONTACT_RADIUS)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, REPRESENTATIVE_VOLUME)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, DenseVector<int>,  NEIGHBOUR_IDS)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, Vector,  NEIGHBOURS_CONTACT_AREAS)
@@ -316,6 +326,7 @@ namespace Kratos
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, PARTICLE_INELASTIC_ROLLING_RESISTANCE_ENERGY)
   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)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, NORMAL_IMPACT_VELOCITY)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, TANGENTIAL_IMPACT_VELOCITY)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, FACE_NORMAL_IMPACT_VELOCITY)

applications/DEMApplication/custom_constitutive/DEM_D_Quadratic_CL.h

@@ -33,7 +33,7 @@ namespace Kratos{
 
         std::string GetTypeOfLaw() override;
 
-        virtual void Check(Properties::Pointer pProp) const override;
+        void Check(Properties::Pointer pProp) const override;
 
         DEMDiscontinuumConstitutiveLaw::Pointer Clone() const override;
 

applications/DEMApplication/custom_constitutive/DEM_D_void_CL.cpp

@@ -0,0 +1,80 @@
+/////////////////////////////////////////////////
+// Author: Chengshun Shang (CIMNE)
+// Email: [email protected], [email protected]
+// Date: June 2023
+/////////////////////////////////////////////////
+
+#include "DEM_D_void_CL.h"
+#include "custom_elements/spheric_particle.h"
+
+namespace Kratos {
+
+    DEMDiscontinuumConstitutiveLaw::Pointer DEM_D_void::Clone() const {
+        DEMDiscontinuumConstitutiveLaw::Pointer p_clone(new DEM_D_void(*this));
+        return p_clone;
+    }
+
+    std::unique_ptr<DEMDiscontinuumConstitutiveLaw> DEM_D_void::CloneUnique() {
+        return Kratos::make_unique<DEM_D_void>();
+    }
+
+    std::string DEM_D_void::GetTypeOfLaw() {
+        std::string type_of_law = "Void"; 
+        return type_of_law;
+    }
+
+    void DEM_D_void::Check(Properties::Pointer pProp) const {}
+
+    /////////////////////////
+    // DEM-DEM INTERACTION //
+    /////////////////////////
+
+    void DEM_D_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]) {
+        KRATOS_TRY
+        LocalElasticContactForce[0]  = 0.0;
+        LocalElasticContactForce[1]  = 0.0;
+        LocalElasticContactForce[2]  = 0.0;
+        ViscoDampingLocalContactForce[0] = 0.0;
+        ViscoDampingLocalContactForce[1] = 0.0;
+        ViscoDampingLocalContactForce[2] = 0.0;
+        KRATOS_CATCH("")
+    }
+
+    /////////////////////////
+    // DEM-FEM INTERACTION //
+    /////////////////////////
+
+    void DEM_D_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) {
+        KRATOS_TRY
+        LocalElasticContactForce[0]  = 0.0;
+        LocalElasticContactForce[1]  = 0.0;
+        LocalElasticContactForce[2]  = 0.0;
+        ViscoDampingLocalContactForce[0] = 0.0;
+        ViscoDampingLocalContactForce[1] = 0.0;
+        ViscoDampingLocalContactForce[2] = 0.0;
+        KRATOS_CATCH("")
+    }
+
+} // namespace Kratos

applications/DEMApplication/custom_constitutive/DEM_D_void_CL.h

@@ -0,0 +1,79 @@
+/////////////////////////////////////////////////
+// Author: Chengshun Shang (CIMNE)
+// Email: [email protected], [email protected]
+// Date: June 2023
+/////////////////////////////////////////////////
+
+#if !defined(DEM_D_VOID_CL_H_INCLUDED)
+#define DEM_D_VOID_CL_H_INCLUDED
+
+#include <string>
+#include <iostream>
+#include "DEM_discontinuum_constitutive_law.h"
+
+namespace Kratos {
+
+    class SphericParticle;
+
+    class KRATOS_API(DEM_APPLICATION) DEM_D_void : public DEMDiscontinuumConstitutiveLaw {
+
+    public:
+
+        KRATOS_CLASS_POINTER_DEFINITION(DEM_D_void);
+
+        DEM_D_void() {}
+
+        ~DEM_D_void() {}
+
+        std::string GetTypeOfLaw() override;
+
+        void Check(Properties::Pointer pProp) const override;
+
+        DEMDiscontinuumConstitutiveLaw::Pointer Clone() const override;
+
+        std::unique_ptr<DEMDiscontinuumConstitutiveLaw> CloneUnique() 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;
+
+    private:
+
+        friend class Serializer;
+
+        virtual void save(Serializer& rSerializer) const override {
+            KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, DEMDiscontinuumConstitutiveLaw)
+                    //rSerializer.save("MyMemberName",myMember);
+        }
+
+        virtual void load(Serializer& rSerializer) override {
+            KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, DEMDiscontinuumConstitutiveLaw)
+                    //rSerializer.load("MyMemberName",myMember);
+        }
+
+    }; //class DEM_D_void
+
+} /* namespace Kratos.*/
+#endif /* DEM_D_VOID_CL_H_INCLUDED  defined */

applications/DEMApplication/custom_constitutive/DEM_compound_constitutive_law_for_PBM.h

@@ -34,6 +34,13 @@ class KRATOS_API(DEM_APPLICATION) DEM_compound_constitutive_law_for_PBM : public
         return mCCL.InitializeContact(element1, element2, indentation);
     }
 
+    void CalculateUnbondedViscoDampingForce(double LocalRelVel[3],
+                                                double UnbondedViscoDampingLocalContactForce[3],
+                                                SphericParticle* const element1,
+                                                SphericParticle* const element2) override{
+        return mCCL.CalculateViscoDampingForce(LocalRelVel, UnbondedViscoDampingLocalContactForce, element1, element2);                                            
+    }
+
 private:
 
     UnbondCL mCCL;

applications/DEMApplication/custom_constitutive/DEM_continuum_constitutive_law.cpp

@@ -10,7 +10,7 @@ namespace Kratos {
     DEMContinuumConstitutiveLaw::~DEMContinuumConstitutiveLaw() {}
 
     std::string DEMContinuumConstitutiveLaw::GetTypeOfLaw() {
-        KRATOS_ERROR << "This function (DEMContinuumConstitutiveLaw::GetTypeOfLaw) shouldn't be accessed, use derived class instead"<<std::endl;
+        //KRATOS_ERROR << "This function (DEMContinuumConstitutiveLaw::GetTypeOfLaw) shouldn't be accessed, use derived class instead"<<std::endl;
         std::string type_of_law = "";
         return type_of_law;
     }