Dem dmt force model (#1045)

New feature
New DMT contact force model.
How Has This Been Tested?
Application test need to be added. I think I'm just going to use the same test as the JKR model. Coming soon
Documentation
Theory has been refactor.
Co-authored-by: Audrey Collard-Daigneault <[email protected]>
Co-authored-by: Victor Oliveira Ferreira <[email protected]>

applications_tests/lethe-particles/pp_dmt_equilibrium.output

@@ -0,0 +1,27 @@
+
+*********************
+Running on 1 rank(s)
+*********************
+DEM time-step is 0.011331% of Rayleigh time step
+Warning: It is recommended to increase the time-step
+Reading triangulation
+
+Finished reading triangulation
+Warning: expansion of particle-wall contact list is disabled.
+This feature is useful in geometries with concave boundaries.
+*********************************************************************
+2 particles of type 0 were inserted, 0 particles of type 0 remaining
+*********************************************************************
+
+*****************************************************************
+Transient iteration: 175000   Time: 1.75     Time step: 1e-05
+*****************************************************************
+| Variable                     | Min        | Max         | Average    | Total      |
+| Contact list generation      | 0.0000e+00 | 1.8000e+01  | 1.8000e+01 | 1.8000e+01 |
+| Velocity magnitude           | 4.6564e-04 | 5.3436e-04  | 5.0000e-04 | 1.0000e-03 |
+| Angular velocity magnitude   | 0.0000e+00 | 2.2251e-308 | 0.0000e+00 | 0.0000e+00 |
+| Translational kinetic energy | 1.1807e-06 | 1.5549e-06  | 1.3678e-06 | 2.7356e-06 |
+| Rotational kinetic energy    | 0.0000e+00 | 2.2251e-308 | 0.0000e+00 | 0.0000e+00 |
+id, type, dp, x, y, z
+0 0 0.20000 1.0086 1.0000 1.0000
+1 0 0.20000 1.1932 1.0000 1.0000

applications_tests/lethe-particles/pp_dmt_equilibrium.prm

@@ -0,0 +1,100 @@
+# Listing of Parameters
+#----------------------
+set dimension = 3
+#---------------------------------------------------
+# Test
+#---------------------------------------------------
+subsection test
+  set enable = true
+end
+#---------------------------------------------------
+# Simulation Control
+#---------------------------------------------------
+subsection simulation control
+  set time step         = 1e-5
+  set time end          = 1.75
+  set log frequency     = 175000
+  set output frequency  = -1
+  set output path       = ./out/
+  set output boundaries = true
+end
+#---------------------------------------------------
+# Model parameters
+#---------------------------------------------------
+subsection model parameters
+  subsection contact detection
+    set contact detection method                = dynamic
+    set dynamic contact search size coefficient = 0.9
+    set neighborhood threshold                  = 20
+  end
+  subsection load balancing
+    set load balance method = none
+  end
+  set particle particle contact force method = DMT
+  set particle wall contact force method     = DMT
+  set integration method                     = velocity_verlet
+end
+#---------------------------------------------------
+# Physical Properties
+#---------------------------------------------------s
+subsection lagrangian physical properties
+  set gx                       = 0.0
+  set gy                       = 0.0
+  set gz                       = 0.0
+  set number of particle types = 1
+  subsection particle type 0
+    set size distribution type            = uniform
+    set diameter                          = 0.2
+    set number of particles               = 2
+    set density particles                 = 2600
+    set young modulus particles           = 1e5
+    set poisson ratio particles           = 0.3
+    set restitution coefficient particles = 0.5
+    set friction coefficient particles    = 0.0
+    set surface energy particles          = 50.
+  end
+  set young modulus wall           = 1000000000000
+  set poisson ratio wall           = 0.30
+  set restitution coefficient wall = 0.8
+  set friction coefficient wall    = 0.
+  set surface energy wall          = 50.
+end
+#---------------------------------------------------
+# Insertion Info
+#---------------------------------------------------
+subsection insertion info
+  set insertion method    = list
+  set insertion frequency = 10000
+  set list x              = 1.0 , 1.2
+  set list y              = 1.0 , 1.0
+  set list z              = 1.0 , 1.0
+  set list velocity x     = 0.001 , 0.0
+  set list velocity y     = 0.0 , 0.0
+  set list velocity z     = 0.0 , 0.0
+  set list omega x        = 0.0 , 0.0
+  set list omega y        = 0.0 , 0.0
+  set list omega z        = 0.0 , 0.0
+  set list diameters      = 0.2 , 0.2
+end
+#---------------------------------------------------
+# Mesh
+#---------------------------------------------------
+subsection mesh
+  set type                                = dealii
+  set grid type                           = hyper_cube
+  set grid arguments                      = 0 : 2 : false
+  set initial refinement                  = 4
+  set expand particle-wall contact search = false
+end
+#---------------------------------------------------
+# Boundary Condition
+#---------------------------------------------------
+subsection DEM boundary conditions
+  set number of boundary conditions = 1
+  subsection boundary condition 0
+    set boundary id       = 0
+    set type              = fixed_wall
+    set rotational speed  = 0
+    set rotational vector = 1,0,0
+  end
+end

applications_tests/lethe-particles/pp_jkr_equilibrium.output

@@ -2,7 +2,7 @@
 *********************
 Running on 1 rank(s)
 *********************
-DEM time-step is 0.00011331% of Rayleigh time step
+DEM time-step is 0.011331% of Rayleigh time step
 Warning: It is recommended to increase the time-step
 Reading triangulation
 
@@ -12,23 +12,15 @@ This feature is useful in geometries with concave boundaries.
 *********************************************************************
 2 particles of type 0 were inserted, 0 particles of type 0 remaining
 *********************************************************************
--->Repartitionning triangulation
-Load balance finished
-Average, minimum and maximum number of particles on the processors are 2 , 2 and 2
-Minimum and maximum number of cells owned by the processors are 4096 and 4096
--->Repartitionning triangulation
-Load balance finished
-Average, minimum and maximum number of particles on the processors are 2 , 2 and 2
-Minimum and maximum number of cells owned by the processors are 4096 and 4096
 
-********************************************************************
-Transient iteration: 17500000 Time: 1.75     Time step: 9.98119e-08
-********************************************************************
+*****************************************************************
+Transient iteration: 175000   Time: 1.75     Time step: 1e-05
+*****************************************************************
 | Variable                     | Min        | Max         | Average    | Total      |
-| Contact list generation      | 0.0000e+00 | 1.7510e+03  | 1.7510e+03 | 1.7510e+03 |
-| Velocity magnitude           | 4.9976e-04 | 5.0024e-04  | 5.0000e-04 | 1.0000e-03 |
+| Contact list generation      | 0.0000e+00 | 1.8000e+01  | 1.8000e+01 | 1.8000e+01 |
+| Velocity magnitude           | 4.1322e-04 | 5.8678e-04  | 5.0000e-04 | 1.0000e-03 |
 | Angular velocity magnitude   | 0.0000e+00 | 2.2251e-308 | 0.0000e+00 | 0.0000e+00 |
-| Translational kinetic energy | 1.3601e-06 | 1.3626e-06  | 1.3614e-06 | 2.7227e-06 |
+| Translational kinetic energy | 9.2979e-07 | 1.8749e-06  | 1.4024e-06 | 2.8047e-06 |
 | Rotational kinetic energy    | 0.0000e+00 | 2.2251e-308 | 0.0000e+00 | 0.0000e+00 |
 id, type, dp, x, y, z
 0 0 0.20000 1.0062 1.0000 1.0000

applications_tests/lethe-particles/pp_jkr_equilibrium.prm

@@ -11,9 +11,9 @@ end
 # Simulation Control
 #---------------------------------------------------
 subsection simulation control
-  set time step         = 1e-7
+  set time step         = 1e-5
   set time end          = 1.75
-  set log frequency     = 17500000
+  set log frequency     = 175000
   set output frequency  = -1
   set output path       = ./out/
   set output boundaries = true
@@ -28,7 +28,7 @@ subsection model parameters
     set neighborhood threshold                  = 20
   end
   subsection load balancing
-    set load balance method = once
+    set load balance method = none
   end
   set particle particle contact force method = hertz_JKR
   set particle wall contact force method     = linear
@@ -97,4 +97,4 @@ subsection DEM boundary conditions
     set rotational speed  = 0
     set rotational vector = 1,0,0
   end
-end
+end

applications_tests/lethe-particles/pw_dmt_equilibrium.output

@@ -0,0 +1,26 @@
+
+*********************
+Running on 1 rank(s)
+*********************
+DEM time-step is 0.11331% of Rayleigh time step
+Warning: It is recommended to increase the time-step
+Reading triangulation
+
+Finished reading triangulation
+Warning: expansion of particle-wall contact list is disabled.
+This feature is useful in geometries with concave boundaries.
+*********************************************************************
+1 particles of type 0 were inserted, 0 particles of type 0 remaining
+*********************************************************************
+
+*****************************************************************
+Transient iteration: 10000    Time: 10       Time step: 0.001
+*****************************************************************
+| Variable                     | Min        | Max         | Average    | Total      |
+| Contact list generation      | 0.0000e+00 | 1.0000e+00  | 1.0000e+00 | 1.0000e+00 |
+| Velocity magnitude           | 7.6385e-04 | 7.6385e-04  | 7.6385e-04 | 7.6385e-04 |
+| Angular velocity magnitude   | 0.0000e+00 | 2.2251e-308 | 0.0000e+00 | 0.0000e+00 |
+| Translational kinetic energy | 3.1772e-03 | 3.1772e-03  | 3.1772e-03 | 3.1772e-03 |
+| Rotational kinetic energy    | 0.0000e+00 | 2.2251e-308 | 0.0000e+00 | 0.0000e+00 |
+id, type, dp, x, y, z
+0 0 2.00000 0.9586 1.0000 1.0000

applications_tests/lethe-particles/pw_dmt_equilibrium.prm

@@ -0,0 +1,99 @@
+# Listing of Parameters
+#----------------------
+set dimension = 3
+#---------------------------------------------------
+# Test
+#---------------------------------------------------
+subsection test
+  set enable = true
+end
+#---------------------------------------------------
+# Simulation Control
+#---------------------------------------------------
+subsection simulation control
+  set time step         = 1e-3
+  set time end          = 10.
+  set log frequency     = 10000
+  set output frequency  = -1
+  set output boundaries = true
+end
+#---------------------------------------------------
+# Model parameters
+#---------------------------------------------------
+subsection model parameters
+  subsection contact detection
+    set contact detection method                = dynamic
+    set dynamic contact search size coefficient = 0.9
+    set neighborhood threshold                  = 20
+  end
+  subsection load balancing
+    set load balance method = none
+  end
+  set particle particle contact force method = DMT
+  set particle wall contact force method     = DMT
+  set integration method                     = velocity_verlet
+end
+#---------------------------------------------------
+# Physical Properties
+#---------------------------------------------------s
+subsection lagrangian physical properties
+  set gx                       = 0.0
+  set gy                       = 0.0
+  set gz                       = 0.0
+  set number of particle types = 1
+  subsection particle type 0
+    set size distribution type            = uniform
+    set diameter                          = 2.
+    set number of particles               = 1
+    set density particles                 = 2600
+    set young modulus particles           = 1e5
+    set poisson ratio particles           = 0.3
+    set restitution coefficient particles = 0.3
+    set friction coefficient particles    = 0.0
+    set surface energy particles          = 50.
+  end
+  set young modulus wall           = 1e5
+  set poisson ratio wall           = 0.30
+  set restitution coefficient wall = 0.4
+  set friction coefficient wall    = 0.
+  set surface energy wall          = 50.
+end
+#---------------------------------------------------
+# Insertion Info
+#---------------------------------------------------
+subsection insertion info
+  set insertion method    = list
+  set insertion frequency = 10000
+  set list x              = 1.000001
+  set list y              = 1.0
+  set list z              = 1.0
+  set list velocity x     = -0.001
+  set list velocity y     = 0.0
+  set list velocity z     = 0.0
+  set list omega x        = 0.0
+  set list omega y        = 0.0
+  set list omega z        = 0.0
+  set list diameters      = 2.0
+end
+#---------------------------------------------------
+# Mesh
+#---------------------------------------------------
+subsection mesh
+  set type                                = dealii
+  set grid type                           = hyper_cube
+  set grid arguments                      = 0 : 4 : false
+  set initial refinement                  = 0
+  set expand particle-wall contact search = false
+end
+#---------------------------------------------------
+# Boundary Condition
+#---------------------------------------------------
+subsection DEM boundary conditions
+  set number of boundary conditions = 1
+  subsection boundary condition 0
+    set boundary id       = 0
+    set type              = fixed_wall
+    set rotational speed  = 0
+    set rotational vector = 1,0,0
+  end
+end

applications_tests/lethe-particles/pw_jkr_equilibrium.output

@@ -2,7 +2,7 @@
 *********************
 Running on 1 rank(s)
 *********************
-DEM time-step is 0.00011331% of Rayleigh time step
+DEM time-step is 0.11331% of Rayleigh time step
 Warning: It is recommended to increase the time-step
 Reading triangulation
 
@@ -12,23 +12,15 @@ This feature is useful in geometries with concave boundaries.
 *********************************************************************
 1 particles of type 0 were inserted, 0 particles of type 0 remaining
 *********************************************************************
--->Repartitionning triangulation
-Load balance finished
-Average, minimum and maximum number of particles on the processors are 1 , 1 and 1
-Minimum and maximum number of cells owned by the processors are 1 and 1
--->Repartitionning triangulation
-Load balance finished
-Average, minimum and maximum number of particles on the processors are 1 , 1 and 1
-Minimum and maximum number of cells owned by the processors are 1 and 1
 
 *****************************************************************
-Transient iteration: 20000000 Time: 20       Time step: 1e-06
+Transient iteration: 10000    Time: 10       Time step: 0.001
 *****************************************************************
 | Variable                     | Min        | Max         | Average    | Total      |
-| Contact list generation      | 0.0000e+00 | 2.0010e+03  | 2.0010e+03 | 2.0010e+03 |
-| Velocity magnitude           | 3.0049e-05 | 3.0049e-05  | 3.0049e-05 | 3.0049e-05 |
+| Contact list generation      | 0.0000e+00 | 1.0000e+00  | 1.0000e+00 | 1.0000e+00 |
+| Velocity magnitude           | 1.5446e-04 | 1.5446e-04  | 1.5446e-04 | 1.5446e-04 |
 | Angular velocity magnitude   | 0.0000e+00 | 2.2251e-308 | 0.0000e+00 | 0.0000e+00 |
-| Translational kinetic energy | 4.9169e-06 | 4.9169e-06  | 4.9169e-06 | 4.9169e-06 |
+| Translational kinetic energy | 1.2992e-04 | 1.2992e-04  | 1.2992e-04 | 1.2992e-04 |
 | Rotational kinetic energy    | 0.0000e+00 | 2.2251e-308 | 0.0000e+00 | 0.0000e+00 |
 id, type, dp, x, y, z
-0 0 2.00000 0.9710 1.0000 1.0000
+0 0 2.00000 0.9716 1.0000 1.0000