Add absolute value to coulomb criterion

chaos-polymtl · Jul 31, 2024 · 9f1ad69 · 9f1ad69
1 parent 63f64ab
commit 9f1ad69
Showing 1 changed file with 23 additions and 23 deletions.
diff --git a/include/dem/particle_particle_contact_force.h b/include/dem/particle_particle_contact_force.h
@@ -384,10 +384,12 @@ class ParticleParticleContactForce
 
     double tangential_damping_constant =
       normal_damping_constant * 0.6324555320336759; // sqrt(0.4)
-    // Calculation of the normal force
-    normal_force = (normal_spring_constant * normal_overlap +
-                    normal_damping_constant * normal_relative_velocity_value) *
-                   normal_unit_vector;
+
+    // Calculation of normal force
+    const double normal_force_norm =
+      normal_spring_constant * normal_overlap +
+      normal_damping_constant * normal_relative_velocity_value;
+    normal_force = normal_force_norm * normal_unit_vector;
 
     // Calculation of tangential force. Since we need damping tangential force
     // in the gross sliding again, we define it as a separate variable
@@ -401,8 +403,7 @@ class ParticleParticleContactForce
     double coulomb_threshold =
       this->effective_coefficient_of_friction[vec_particle_type_index(
         particle_one_type, particle_two_type)] *
-      normal_force.norm();
-
+      std::fabs(normal_force_norm);
 
     // Check for gross sliding
     if (tangential_force.norm() > coulomb_threshold)
@@ -556,7 +557,7 @@ class ParticleParticleContactForce
     double coulomb_threshold =
       this->effective_coefficient_of_friction[vec_particle_type_index(
         particle_one_type, particle_two_type)] *
-      normal_force_norm;
+      std::fabs(normal_force_norm);
 
     // Check for gross sliding
     const double tangential_force_norm = tangential_force.norm();
@@ -693,11 +694,11 @@ class ParticleParticleContactForce
     double tangential_damping_constant =
       normal_damping_constant * sqrt(model_parameter_st / model_parameter_sn);
 
-    // Calculation of normal force using spring and dashpot normal forces
-    normal_force =
-      ((normal_spring_constant * normal_overlap) * normal_unit_vector) +
-      ((normal_damping_constant * normal_relative_velocity_value) *
-       normal_unit_vector);
+    // Calculation of normal force
+    const double normal_force_norm =
+      normal_spring_constant * normal_overlap +
+      normal_damping_constant * normal_relative_velocity_value;
+    normal_force = normal_force_norm * normal_unit_vector;
 
     // Calculation of tangential force using spring and dashpot tangential
     // forces. Since we need dashpot tangential force in the gross sliding
@@ -709,7 +710,7 @@ class ParticleParticleContactForce
     double coulomb_threshold =
       this->effective_coefficient_of_friction[vec_particle_type_index(
         particle_one_type, particle_two_type)] *
-      normal_force.norm();
+      std::fabs(normal_force_norm);
 
     // Check for gross sliding
     if (tangential_force.norm() > coulomb_threshold)
@@ -834,11 +835,10 @@ class ParticleParticleContactForce
 
 
     // Calculation of normal force using spring and dashpot normal forces
-    normal_force =
-      ((normal_spring_constant * normal_overlap) * normal_unit_vector) +
-      ((normal_damping_constant * normal_relative_velocity_value) *
-       normal_unit_vector);
-
+    const double normal_force_norm =
+      normal_spring_constant * normal_overlap +
+      normal_damping_constant * normal_relative_velocity_value;
+    normal_force = normal_force_norm * normal_unit_vector;
     // Calculation of tangential force using spring and dashpot tangential
     // forces. Since we need dashpot tangential force in the gross sliding
     // again, we define it as a separate variable
@@ -848,7 +848,7 @@ class ParticleParticleContactForce
     double coulomb_threshold =
       this->effective_coefficient_of_friction[vec_particle_type_index(
         particle_one_type, particle_two_type)] *
-      normal_force.norm();
+      std::fabs(normal_force_norm);
 
     // Check for gross sliding
     if (tangential_force.norm() > coulomb_threshold)
@@ -1031,9 +1031,9 @@ class ParticleParticleContactForce
         particle_one_type, particle_two_type)] *
       this->effective_radius;
     const double modified_coulomb_threshold =
-      (normal_force_coefficient + two_pull_off_force) *
-      this->effective_coefficient_of_friction[vec_particle_type_index(
-        particle_one_type, particle_two_type)];
+      std::fabs((normal_force_coefficient + two_pull_off_force) *
+                this->effective_coefficient_of_friction[vec_particle_type_index(
+                  particle_one_type, particle_two_type)]);
 
     if (tangential_force.norm() > modified_coulomb_threshold)
       {
@@ -1180,7 +1180,7 @@ class ParticleParticleContactForce
     double coulomb_threshold =
       this->effective_coefficient_of_friction[vec_particle_type_index(
         particle_one_type, particle_two_type)] *
-      normal_force_norm;
+      std::fabs(normal_force_norm);
 
     // Check for gross sliding
     const double tangential_force_norm = tangential_force.norm();