Fix the Coulomb's criterion in linear, Hertz, and Hertz-Mindlin with …

…limit force in DEM (#1216) Description The Coulomb's criterion was wrong in the particle-particle contact force for Hertz-Mindlin with limit force, Hertz and Linear. The normal force norm was explicitly positive when doing normal_force.norm(), making the Coulomb's criterion always positive even if the particles are in repulsion, so in slidling. Solution Apply the the same way the Hertz-Mindlin with limit overlap is calculated. Testing Tests with low restitution coefficient may give negative coulomb criterion, this is why some tests have changed. Co-authored-by: Bruno Blais <[email protected]> Co-authored-by: OGaboriault <[email protected]> Former-commit-id: 38900e1
chaos-polymtl · Jul 31, 2024 · c6bd61d · c6bd61d
1 parent bd08437
commit c6bd61d
Show file tree

Hide file tree

Showing 9 changed files with 132 additions and 129 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -5,17 +5,20 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/).
 
 ## [Master] - 2024-07-31
 
+### Fixed
+
+- MINOR The Coulomb's criterion was wrong in the particle-particle contact force for Hertz-Mindlin with limit force, Hertz and Linear in DEM. The normal force norm was explicitly positive when doing normal_force.norm(), making the Coulomb's criterion always positive even if the particles are in repulsion, so in slidling. This has been fixed using the same method as for Hertz-Mindlin with limit overlap is calculated. [#1216](https://github.com/chaos-polymtl/lethe/pull/1216)
+
 ### Added
 
 - MINOR P-multigrid was added to the gcmg preconditioner of the lethe-fluid-matrix-free application. It supports three different coarsening strategies to define the degree p of the different levels. It also allows to use hybrid hp- and ph-multigrid strategies. [#1209](https://github.com/chaos-polymtl/lethe/pull/1209)
 
-
 ## [Master] - 2024-07-24
 
 ### Changed
 
 - MINOR Forces a contact search at the last DEM iteration of a CFD iteration for more robustness related to the update of the reference location of the particles prior the void fraction calculation [#1205](https://github.com/chaos-polymtl/lethe/pull/1205)
-- 
+
 ## [Master] - 2024-07-23
 
 ### Changed

diff --git a/applications_tests/lethe-particles/ball_with_floating_walls.output b/applications_tests/lethe-particles/ball_with_floating_walls.output
@@ -12,13 +12,13 @@ This feature should only be activated in geometries with concave boundaries. (Fo
 10 particles of type 0 were inserted, 0 particles of type 0 remaining
 **********************************************************************
 id, type, dp, x, y, z 
-0 0 0.00500 0.0051 -0.0247 -0.0915
-1 0 0.00500 -0.0044 -0.0148 -0.0936
-2 0 0.00500 -0.0008 -0.0094 -0.0953
-3 0 0.00500 0.0089 -0.0215 -0.0914
-4 0 0.00500 0.0003 -0.0162 -0.0941
-5 0 0.00500 0.0045 -0.0187 -0.0928
-6 0 0.00500 -0.0003 -0.0034 -0.0966
-7 0 0.00500 0.0046 -0.0038 -0.0957
-8 0 0.00500 0.0075 -0.0077 -0.0944
-9 0 0.00500 0.0054 -0.0125 -0.0939
+0 0 0.00500 0.0062 -0.0244 -0.0914
+1 0 0.00500 -0.0034 -0.0142 -0.0939
+2 0 0.00500 -0.0008 -0.0093 -0.0954
+3 0 0.00500 0.0103 -0.0215 -0.0911
+4 0 0.00500 0.0012 -0.0161 -0.0940
+5 0 0.00500 0.0038 -0.0202 -0.0927
+6 0 0.00500 -0.0004 -0.0034 -0.0966
+7 0 0.00500 0.0045 -0.0025 -0.0959
+8 0 0.00500 0.0071 -0.0083 -0.0943
+9 0 0.00500 0.0052 -0.0131 -0.0938
diff --git a/applications_tests/lethe-particles/box_grid_rotation.output b/applications_tests/lethe-particles/box_grid_rotation.output
@@ -12,31 +12,31 @@ This feature is useful in geometries with concave boundaries.
 100 particles of type 0 were inserted, 100 particles of type 0 remaining
 *************************************************************************
 id, type, dp, x, y, z 
-0 0 0.00500 -0.0304 -0.0312 -0.0408
-1 0 0.00500 -0.0246 -0.0325 -0.0409
-2 0 0.00500 -0.0159 -0.0315 -0.0408
-3 0 0.00500 -0.0069 -0.0314 -0.0408
-4 0 0.00500 -0.0021 -0.0327 -0.0410
-5 0 0.00500 0.0065 -0.0336 -0.0411
-6 0 0.00500 0.0120 -0.0302 -0.0407
-7 0 0.00500 0.0231 -0.0328 -0.0410
-8 0 0.00500 0.0307 -0.0336 -0.0411
-9 0 0.00500 -0.0302 -0.0248 -0.0402
-10 0 0.00500 -0.0222 -0.0259 -0.0403
-11 0 0.00500 -0.0158 -0.0265 -0.0403
-12 0 0.00500 -0.0115 -0.0238 -0.0401
-13 0 0.00500 0.0016 -0.0254 -0.0402
+0 0 0.00500 -0.0304 -0.0313 -0.0408
+1 0 0.00500 -0.0245 -0.0324 -0.0409
+2 0 0.00500 -0.0160 -0.0317 -0.0408
+3 0 0.00500 -0.0068 -0.0314 -0.0408
+4 0 0.00500 -0.0020 -0.0326 -0.0410
+5 0 0.00500 0.0066 -0.0336 -0.0411
+6 0 0.00500 0.0121 -0.0300 -0.0407
+7 0 0.00500 0.0233 -0.0329 -0.0410
+8 0 0.00500 0.0309 -0.0336 -0.0411
+9 0 0.00500 -0.0302 -0.0249 -0.0402
+10 0 0.00500 -0.0220 -0.0262 -0.0403
+11 0 0.00500 -0.0157 -0.0268 -0.0404
+12 0 0.00500 -0.0116 -0.0238 -0.0401
+13 0 0.00500 0.0017 -0.0255 -0.0403
 14 0 0.00500 0.0073 -0.0230 -0.0400
-15 0 0.00500 0.0122 -0.0251 -0.0402
-16 0 0.00500 0.0205 -0.0251 -0.0402
-17 0 0.00500 0.0300 -0.0236 -0.0401
-18 0 0.00500 -0.0323 -0.0176 -0.0395
+15 0 0.00500 0.0124 -0.0250 -0.0402
+16 0 0.00500 0.0205 -0.0250 -0.0402
+17 0 0.00500 0.0300 -0.0237 -0.0401
+18 0 0.00500 -0.0324 -0.0177 -0.0395
 19 0 0.00500 -0.0229 -0.0167 -0.0394
-20 0 0.00500 -0.0156 -0.0159 -0.0393
+20 0 0.00500 -0.0156 -0.0158 -0.0393
 21 0 0.00500 -0.0080 -0.0170 -0.0394
 22 0 0.00500 -0.0003 -0.0172 -0.0394
 23 0 0.00500 0.0064 -0.0169 -0.0394
-24 0 0.00500 0.0142 -0.0161 -0.0393
+24 0 0.00500 0.0142 -0.0162 -0.0393
 25 0 0.00500 0.0222 -0.0173 -0.0394
 26 0 0.00500 0.0294 -0.0175 -0.0394
 27 0 0.00500 -0.0304 -0.0092 -0.0386
@@ -93,22 +93,22 @@ id, type, dp, x, y, z
 78 0 0.00500 0.0143 0.0280 -0.0349
 79 0 0.00500 0.0215 0.0278 -0.0349
 80 0 0.00500 0.0297 0.0281 -0.0349
-81 0 0.00500 -0.0310 -0.0340 -0.0367
-82 0 0.00500 -0.0204 -0.0305 -0.0390
-83 0 0.00500 -0.0166 -0.0340 -0.0366
-84 0 0.00500 -0.0113 -0.0336 -0.0406
-85 0 0.00500 0.0021 -0.0304 -0.0394
-86 0 0.00500 0.0068 -0.0282 -0.0405
-87 0 0.00500 0.0166 -0.0336 -0.0411
-88 0 0.00500 0.0192 -0.0297 -0.0391
-89 0 0.00500 0.0264 -0.0285 -0.0400
+81 0 0.00500 -0.0309 -0.0340 -0.0367
+82 0 0.00500 -0.0203 -0.0307 -0.0385
+83 0 0.00500 -0.0166 -0.0341 -0.0364
+84 0 0.00500 -0.0113 -0.0336 -0.0410
+85 0 0.00500 0.0021 -0.0304 -0.0391
+86 0 0.00500 0.0067 -0.0281 -0.0405
+87 0 0.00500 0.0163 -0.0336 -0.0410
+88 0 0.00500 0.0193 -0.0298 -0.0397
+89 0 0.00500 0.0265 -0.0285 -0.0400
 90 0 0.00500 -0.0308 -0.0222 -0.0360
-91 0 0.00500 -0.0263 -0.0207 -0.0393
-92 0 0.00500 -0.0165 -0.0215 -0.0397
-93 0 0.00500 -0.0055 -0.0260 -0.0399
-94 0 0.00500 -0.0021 -0.0230 -0.0379
-95 0 0.00500 0.0068 -0.0262 -0.0359
-96 0 0.00500 0.0170 -0.0217 -0.0395
-97 0 0.00500 0.0250 -0.0231 -0.0395
-98 0 0.00500 0.0293 -0.0268 -0.0362
-99 0 0.00500 -0.0275 -0.0150 -0.0388
+91 0 0.00500 -0.0262 -0.0208 -0.0393
+92 0 0.00500 -0.0165 -0.0213 -0.0395
+93 0 0.00500 -0.0056 -0.0261 -0.0399
+94 0 0.00500 -0.0021 -0.0229 -0.0383
+95 0 0.00500 0.0068 -0.0261 -0.0360
+96 0 0.00500 0.0169 -0.0216 -0.0394
+97 0 0.00500 0.0250 -0.0230 -0.0396
+98 0 0.00500 0.0293 -0.0267 -0.0361
+99 0 0.00500 -0.0276 -0.0150 -0.0392
diff --git a/applications_tests/lethe-particles/circle_restart.mpirun=1.output b/applications_tests/lethe-particles/circle_restart.mpirun=1.output
@@ -14,21 +14,21 @@ This feature should only be activated in geometries with concave boundaries. (Fo
 id, type, dp, x, y, z 
 0 0 0.00500 -0.0191 -0.0956
 1 0 0.00500 -0.0141 -0.0961
-2 0 0.00500 -0.0092 -0.0966
-3 0 0.00500 -0.0042 -0.0971
-4 0 0.00500 0.0011 -0.0974
-5 0 0.00500 0.0062 -0.0969
-6 0 0.00500 0.0112 -0.0964
-7 0 0.00500 0.0161 -0.0959
-8 0 0.00500 0.0210 -0.0950
-9 0 0.00500 -0.0287 -0.0927
-10 0 0.00500 -0.0239 -0.0941
+2 0 0.00500 -0.0091 -0.0966
+3 0 0.00500 -0.0041 -0.0971
+4 0 0.00500 0.0012 -0.0974
+5 0 0.00500 0.0064 -0.0969
+6 0 0.00500 0.0114 -0.0964
+7 0 0.00500 0.0164 -0.0959
+8 0 0.00500 0.0213 -0.0949
+9 0 0.00500 -0.0286 -0.0927
+10 0 0.00500 -0.0239 -0.0942
 11 0 0.00500 -0.0112 -0.0920
-12 0 0.00500 -0.0063 -0.0925
+12 0 0.00500 -0.0062 -0.0925
 13 0 0.00500 -0.0013 -0.0930
-14 0 0.00500 0.0037 -0.0926
-15 0 0.00500 0.0086 -0.0921
-16 0 0.00500 0.0258 -0.0936
-17 0 0.00500 0.0306 -0.0921
-18 0 0.00500 0.0354 -0.0907
-19 0 0.00500 -0.0162 -0.0915
+14 0 0.00500 0.0037 -0.0927
+15 0 0.00500 0.0087 -0.0922
+16 0 0.00500 0.0261 -0.0935
+17 0 0.00500 0.0308 -0.0920
+18 0 0.00500 0.0356 -0.0906
+19 0 0.00500 -0.0162 -0.0916
diff --git a/applications_tests/lethe-particles/circle_with_floating_walls.mpirun=1.output b/applications_tests/lethe-particles/circle_with_floating_walls.mpirun=1.output
@@ -12,21 +12,21 @@ This feature should only be activated in geometries with concave boundaries. (Fo
 18 particles of type 0 were inserted, 0 particles of type 0 remaining
 **********************************************************************
 id, type, dp, x, y, z 
-0 0 0.00500 -0.0121 -0.0963
-1 0 0.00500 -0.0070 -0.0968
+0 0 0.00500 -0.0122 -0.0963
+1 0 0.00500 -0.0069 -0.0968
 2 0 0.00500 -0.0017 -0.0973
 3 0 0.00500 0.0034 -0.0972
 4 0 0.00500 0.0084 -0.0967
 5 0 0.00500 0.0135 -0.0962
-6 0 0.00500 -0.0366 -0.0335
-7 0 0.00500 -0.0314 -0.0335
-8 0 0.00500 -0.0153 -0.0335
-9 0 0.00500 0.0022 -0.0335
-10 0 0.00500 0.0236 -0.0335
-11 0 0.00500 0.0294 -0.0335
-12 0 0.00500 -0.0573 -0.0335
-13 0 0.00500 -0.0211 -0.0335
-14 0 0.00500 -0.0101 -0.0335
-15 0 0.00500 -0.0046 -0.0335
-16 0 0.00500 0.0076 -0.0335
-17 0 0.00500 0.0429 -0.0335
+6 0 0.00500 -0.0376 -0.0335
+7 0 0.00500 -0.0312 -0.0335
+8 0 0.00500 -0.0157 -0.0335
+9 0 0.00500 0.0027 -0.0335
+10 0 0.00500 0.0226 -0.0335
+11 0 0.00500 0.0288 -0.0335
+12 0 0.00500 -0.0570 -0.0335
+13 0 0.00500 -0.0214 -0.0335
+14 0 0.00500 -0.0107 -0.0335
+15 0 0.00500 -0.0050 -0.0335
+16 0 0.00500 0.0081 -0.0335
+17 0 0.00500 0.0431 -0.0335
diff --git a/applications_tests/lethe-particles/packing_in_ball_dynamic_contact.output b/applications_tests/lethe-particles/packing_in_ball_dynamic_contact.output
@@ -23,7 +23,7 @@ id, type, dp, x, y, z
 8 0 0.00500 0.0106 -0.0360 -0.0809
 9 0 0.00500 0.0163 -0.0357 -0.0797
 10 0 0.00500 0.0211 -0.0353 -0.0790
-11 0 0.00500 0.0271 -0.0357 -0.0783
+11 0 0.00500 0.0271 -0.0357 -0.0784
 12 0 0.00500 0.0327 -0.0352 -0.0769
 13 0 0.00500 -0.0358 -0.0283 -0.0781
 14 0 0.00500 -0.0301 -0.0288 -0.0800

diff --git a/applications_tests/lethe-particles/packing_in_ball_dynamic_load_balance.mpirun=2.output b/applications_tests/lethe-particles/packing_in_ball_dynamic_load_balance.mpirun=2.output
@@ -26,12 +26,12 @@ id, type, dp, x, y, z
 3 0 0.00500 -0.0186 -0.0353 -0.0790
 4 0 0.00500 -0.0132 -0.0364 -0.0809
 5 0 0.00500 -0.0068 -0.0360 -0.0808
-6 0 0.00500 -0.0003 -0.0367 -0.0817
-7 0 0.00500 0.0051 -0.0360 -0.0811
+6 0 0.00500 -0.0004 -0.0367 -0.0817
+7 0 0.00500 0.0051 -0.0359 -0.0811
 8 0 0.00500 0.0106 -0.0360 -0.0809
 9 0 0.00500 0.0163 -0.0357 -0.0797
 10 0 0.00500 0.0211 -0.0353 -0.0790
-11 0 0.00500 0.0271 -0.0357 -0.0783
+11 0 0.00500 0.0271 -0.0357 -0.0784
 12 0 0.00500 0.0327 -0.0352 -0.0769
 13 0 0.00500 -0.0358 -0.0283 -0.0781
 14 0 0.00500 -0.0301 -0.0288 -0.0800

diff --git a/applications_tests/lethe-particles/two_types_packing_in_circle.output b/applications_tests/lethe-particles/two_types_packing_in_circle.output
@@ -15,23 +15,23 @@ This feature should only be activated in geometries with concave boundaries. (Fo
 10 particles of type 1 were inserted, 0 particles of type 1 remaining
 **********************************************************************
 id, type, dp, x, y, z 
-0 0 0.00300 -0.0056 -0.0980
-1 0 0.00300 -0.0022 -0.0983
-2 0 0.00300 0.0122 -0.0943
-3 0 0.00300 0.0037 -0.0981
-4 0 0.00300 0.0008 -0.0984
-5 0 0.00300 0.0066 -0.0974
-6 0 0.00300 0.0096 -0.0975
-7 0 0.00300 0.0126 -0.0973
-8 0 0.00300 0.0156 -0.0970
-9 0 0.00300 0.0157 -0.0940
-10 1 0.00500 -0.0144 -0.0961
-11 1 0.00500 -0.0094 -0.0966
-12 1 0.00500 -0.0082 -0.0917
-13 1 0.00500 -0.0036 -0.0946
-14 1 0.00500 0.0019 -0.0946
-15 1 0.00500 -0.0127 -0.0703
-16 1 0.00500 0.0075 -0.0932
-17 1 0.00500 -0.0118 -0.0785
-18 1 0.00500 0.0280 -0.0745
-19 1 0.00500 0.0194 -0.0955
+0 0 0.00300 -0.0063 -0.0979
+1 0 0.00300 -0.0002 -0.0955
+2 0 0.00300 0.0036 -0.0954
+3 0 0.00300 -0.0033 -0.0982
+4 0 0.00300 0.0015 -0.0930
+5 0 0.00300 -0.0003 -0.0985
+6 0 0.00300 0.0026 -0.0982
+7 0 0.00300 0.0056 -0.0979
+8 0 0.00300 0.0209 -0.0930
+9 0 0.00300 0.0087 -0.0976
+10 1 0.00500 -0.0202 -0.0952
+11 1 0.00500 -0.0153 -0.0960
+12 1 0.00500 -0.0086 -0.0917
+13 1 0.00500 -0.0043 -0.0942
+14 1 0.00500 -0.0101 -0.0965
+15 1 0.00500 0.0280 -0.0929
+16 1 0.00500 0.0071 -0.0597
+17 1 0.00500 0.0083 -0.0861
+18 1 0.00500 0.0125 -0.0963
+19 1 0.00500 0.0333 -0.0913
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_value =
+      normal_spring_constant * normal_overlap +
+      normal_damping_constant * normal_relative_velocity_value;
+    normal_force = normal_force_value * 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
@@ -398,11 +400,10 @@ class ParticleParticleContactForce
       (tangential_spring_constant * contact_info.tangential_overlap) +
       damping_tangential_force;
 
-    double coulomb_threshold =
+    const double coulomb_threshold =
       this->effective_coefficient_of_friction[vec_particle_type_index(
         particle_one_type, particle_two_type)] *
-      normal_force.norm();
-
+      normal_force_value;
 
     // Check for gross sliding
     if (tangential_force.norm() > coulomb_threshold)
@@ -540,10 +541,10 @@ class ParticleParticleContactForce
       normal_damping_constant * sqrt(model_parameter_st / model_parameter_sn);
 
     // Calculation of normal force
-    const double normal_force_norm =
+    const double normal_force_value =
       normal_spring_constant * normal_overlap +
       normal_damping_constant * normal_relative_velocity_value;
-    normal_force = normal_force_norm * normal_unit_vector;
+    normal_force = normal_force_value * 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
@@ -553,10 +554,10 @@ class ParticleParticleContactForce
       (tangential_spring_constant * contact_info.tangential_overlap) +
       damping_tangential_force;
 
-    double coulomb_threshold =
+    const double coulomb_threshold =
       this->effective_coefficient_of_friction[vec_particle_type_index(
         particle_one_type, particle_two_type)] *
-      normal_force_norm;
+      normal_force_value;
 
     // 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_value =
+      normal_spring_constant * normal_overlap +
+      normal_damping_constant * normal_relative_velocity_value;
+    normal_force = normal_force_value * normal_unit_vector;
 
     // Calculation of tangential force using spring and dashpot tangential
     // forces. Since we need dashpot tangential force in the gross sliding
@@ -706,10 +707,10 @@ class ParticleParticleContactForce
       tangential_spring_constant * contact_info.tangential_overlap +
       tangential_damping_constant * tangential_relative_velocity;
 
-    double coulomb_threshold =
+    const double coulomb_threshold =
       this->effective_coefficient_of_friction[vec_particle_type_index(
         particle_one_type, particle_two_type)] *
-      normal_force.norm();
+      normal_force_value;
 
     // Check for gross sliding
     if (tangential_force.norm() > coulomb_threshold)
@@ -834,21 +835,20 @@ 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_value =
+      normal_spring_constant * normal_overlap +
+      normal_damping_constant * normal_relative_velocity_value;
+    normal_force = normal_force_value * 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
     tangential_force =
       tangential_spring_constant * contact_info.tangential_overlap;
 
-    double coulomb_threshold =
+    const double coulomb_threshold =
       this->effective_coefficient_of_friction[vec_particle_type_index(
         particle_one_type, particle_two_type)] *
-      normal_force.norm();
+      normal_force_value;
 
     // Check for gross sliding
     if (tangential_force.norm() > coulomb_threshold)