Discrete Element Method in Lethe

Launching a DEM simulation in Lethe requires an executable of the required solver (dem_2d or dem_3d for two-dimensional and three-dimensional simulations, respectively), and a parameters file (with extension .prm if it uses the parameter file format, or .json if it uses the JSON file format). The parameter file controls all aspects of the simulations. This section aims at describing the various parameters available within Lethe-DEM. The readers are recommended to visit here before reading the following.

Parameters file (.prm)

The parameter file is composed of different subsections. In the following, the principal subsections of a DEM parameter template file are explained.

💡 In the parameter file format, the parameters are established one by one using the following syntax example: set time end = 10 sets the end time to 10 seconds.

💡 The arguments can be either double, integer, string, or a choice between a predefined number of variables. In the parameter files, comments are preceded by the sharp symbol (e.g. '#comment').

Simulation Control
Timer
Test
Model Parameters
Physical Properties
Insertion Info
Floating walls
Mesh
Boundary Conditions
Boundary Forces
Triangulation Motion

Simulation Control

This subsection contains the general information of the simulation, including time step, end time, output path, log frequency, and output frequency. It is the most commonly modified section for a simulation. time step in DEM simulations is generally in the range of 1e-7 to 1e-5 seconds. With this small time step, the DEM simulation can capture a single collision in several iterations, which leads to the high accuracy of the simulation.

❗ Lethe-DEM compares the selected time step to Rayleigh time-step and prints a warning if the selected time step is too small or too large. Users should modify the time step if they see this warning.

❗ It should be mentioned that if the output path is specified in the parameter handler file, a folder with the specified name should be created before starting simulation in the working directory.

subsection simulation control
  # DEM time-step 
  set time step                 			 = 1e-5

  # Simulation end time
  set time end       					 = 0.2

  # File output prefix
  set output path                                        = ./

  # Log frequency
  set log frequency				         = 1000

  # Output frequency
  set output frequency            			 = 10000
end

Timer

This subsection deals with simulation timing information. If it is set to end, the simulation timing information will be printed at the end of the simulation.

subsection timer
  set type    					         = end
end

Test

This subsection deals with simulation testing. If it is set to true, the sorted (using particle id) information of all particles will be printed in xyz format at the end of the simulation. Then this information can be compared with predefined values to ensure the reproducibility of the simulations.

subsection test
  # DEM test
  # Choices are true|false
  set type    					         = false
end

Model Parameters

In this subsection, DEM simulation parameters are defined. These parameters include contact detection method and its subsequent information (dynamic contact search size coefficient or contact detection frequency for dynamic or constant contact detection method), repartition method (once, frequent or dynamic) for load balancing (in parallel simulations), particle weight which defines the particle weight compared to a default cell weight of 1000, neighborhood threshold (which defines the contact neighbor list size: neighborhood threshold * particle diameter), particle particle contact force method, particle wall contact force method and integration method.

💡 In constant contact search, Lethe-DEM searches for particle-particle and particle-wall collisions at a constant frequency (every contact detection frequency iterations). Normally, contact detection frequency should be selected in the range of 5-50. Small values of contact detection frequency lead to long simulation times, while large values of contact detection frequency may lead to late detection of collisions. Late detection of collisions can result in very large particles velocities (popcorn jump of particles in a simulation) or particles leaving the simulation domain.

💡 In dynamic contact search, Lethe-DEM searches for particle-particle and particle-wall collisions automatically. In dynamic contact search, Lethe-DEM predicts (using the velocities of particles at each iteration) and stores (adds up) the displacements of each particle in the simulation. If the maximum displacement of particles exceeds the smallest contact search criterion (explained in the following), then the iteration is recognized as a contact detection iteration.

❕ smallest contact search criterion is defined as:

$\phi=\min({d_c^{min}r_p^{max},\epsilon(\alpha-1)r_p^{max}})$

where $\phi$ , $d_c^{min}$ , $r_p^{max}$ , $\epsilon$ , $\alpha$ denote smallest contact search criterion, minimum cell size (in the triangulation), maximum particle radius (in polydisperse simulations), dynamic contact search size coefficient, and neighborhood threshold.

❗ dynamic contact search size coefficient, as illustrated in the equation above, is a safety factor to ensure the late detection of particles will not happen in the simulations with dynamic contact search; and its value should be defined generally in the range of 0.5-1. 0.5 is a rather conservative value for dynamic contact search size coefficient.

💡 To optimize the computational performance of Lethe-DEM, the contact search is performed in two steps (broad and fine contact searches). In broad contact search, neighbor particles in adjacent cells are stored in a list. Then this list of adjacent particles is refined once more in the fine search. More information about the contact search algorithm in Lethe-DEM can be found here. In the fine search, the particles (from the output list of broad search) which are located in a spherical region around each particle are stored in another list and sent to force calculation class. The diameter of this spherical region ( around each particle is defined as:

$D={\alpha}d_p^{max}$

where $d_p^{max}$ denotes the maximum particle size in the simulation.

❗ According to the definition of the spherical region in the fine search and contact force (equation above), the value of the neighborhood threshold ( $\alpha$ ) must be larger than 1. It is generally defined in the range of 1.3-1.5.

💡 Load-balancing updates the distribution of the subdomains between the processes in parallel simulation to achieve better computational performance (less simulation time). Three load-balancing methods are available in Lethe-DEM: once, frequent, or dynamic. Read here for more information about different load-balancing methods and their performances in various types of DEM simulations. The total weight of each cell with particles in load-balancing is defined as:

W=1000+W_pn_p

where W_p is the particle weight and n_p is the number of particles in the cell. 1000 is the default weight assigned to one cell.

❗ Selecting repartition method = once, requires defining the step at which the code calls load balancing (load balance step). Dynamic repartition method requires defining load balance frequency, and in dynamic repartition method, we should define load balance threshold and dynamic load balance check frequency. In dynamic load balancing, the code checks the distribution of particles among the processors, every dynamic load balance check frequency steps, and if $L_{max}-L_{min}>{\beta}\bar{L}$ it calls load-balancing. and $\beta$ denote computational load on a process and load balance threshold, respectively.

💡 Four particle particle contact models are available in Lethe-DEM (hertz_mindlin_limit_overlap, hertz_mindlin_limit_force, hertz, and linear). hertz_mindlin_limit_overlap and hertz_mindlin_limit_force are non-linear contact models in which the stiffness and damping forces in both normal and tangential directions are considered. The only difference between these models is in their limiting method of the tangential force during gross sliding (where the tangential force exceeds the coulomb's limit). In hertz_mindlin_limit_overlap model, Lethe-DEM limits the tangential overlap and with limiting the overlap, the tangential force is limited; while in hertz_mindlin_limit_force model, the tangential force is limited directly without limiting the tangential overlap. hertz is another non-linear model in which the damping force is not considered in the tangential direction and the tangential force is limited in gross sliding. Lethe-DEM also has a linear contact model (the stiffness and damping forces are linear functions of overlap and relative velocity, respectively).

💡 Lethe-DEM has two linear and non-linear particle-wall models.

💡 euler (1st order) and velocity-verlet (2nd order) are the available integration methods in Lethe-DEM.

:bulb" Three rolling resistance models are available in Lethe-DEM: no_resistance, constant_resistance, viscous_resistance.

subsection model parameters
  # Contact detection method
  # Choices are constant|dynamic
  set contact detection method   			 = dynamic

  # Depending on the contact detection method, contact search size coefficient (safety factor multiplier for dynamic contact search) or contact detection frequency should be defined for dynamic and constant contact search methods, respectively.
  set dynamic contact search size coefficient	 	 = 0.9
  set contact detection frequency                 	 = 20

  # Load balancing method and its subsequent information
  set load balance method				 = frequent
  set load balance frequency				 = 200000
  set load balance particle weight                       = 10000

  # Particle-particle contact neighborhood size
  set neighborhood threshold				 = 1.6

  # Particle-particle contact force model
  # Choices are linear|hertz_mindlin_limit_overlap|hertz_mindlin_limit_force|hertz
  set particle particle contact force method             = hertz_mindlin_limit_overlap

  # Particle-wall contact force model
  # Choices are linear|nonlinear
  set particle wall contact force method                 = nonlinear

  # Integration method
  # Choices are euler|velocity_verlet
  set integration method				 = velocity_verlet

  # Rolling resistance method
  # choices are no_resistance|constant_resistance|viscous_resistance
  set rolling resistance torque method			 = no_resistance
end

Lagrangian Physical Properties

In this subsection, the physical properties of the simulation are defined. These properties include gravitational acceleration (gx, gy and gz), number of particle types, and for each type, particle diameter, particle density, Young's modulus of particle and wall, Poisson ratio of particle and wall, restitution coefficient of particle and wall, friction coefficient of particle and wall and rolling friction coefficient of particle and wall.

💡 If the particles in the simulation are monodispersed, the number of particle types should be equal to zero. For polydispersed systems, the number of particle types is selected equal to the number of particles types in the simulation. For each particle type, a separate subsection particle type n should be defined (n starts from zero to number of particle types - 1) which contains all the physical properties related to that particle type.

❗ For each particle type, two size distribution types can be defined: uniform and normal. In uniform size distribution, the diameter of the particles is constant, while in normal size distribution, the particle diameters are sampled from a normal distribution with an average of average diameter and standard deviation of standard deviation.

subsection model parameters
  # Gravitational acceleration in x direction
  set gx            		 			 = 0.0

  # Gravitational acceleration in y direction
  set gy            		 			 = 0.0

  # Gravitational acceleration in z direction
  set gz            		 			 = -9.81

  # Number of particle types
  set number of particle types	                         = 1

  # Entering particle type 0
    	subsection particle type 0

  # Size distribution of particle type 0
		set size distribution type		= uniform

  # Particle diameter
                set diameter            	 	= 0.005

  # Number of particles in type 0
		set number				= 132300

  # Particle density
                set density particles           	= 2000

  # Young's modulus of particle
                set young modulus particles         	= 1000000

  # Poisson ratio of particle
                set poisson ratio particles            	= 0.3

  # Coefficient of restitution of particle
                set restitution coefficient particles   = 0.95

  # Coefficient of friction of particle
                set friction coefficient particles      = 0.05

  # Coefficient of rolling friction of particle
                set rolling friction particles          = 0.1
	end

  # Young's modulus of wall
  set young modulus wall            			 = 1000000

  # Poisson ratio of wall
  set poisson ratio wall            			 = 0.3

  # Coefficient of restitution of wall
  set restitution coefficient wall           		 = 0.95

  # Coefficient of friction of wall
  set friction coefficient wall         		 = 0.05

  # Coefficient of rolling friction of wall
  set rolling friction wall         	      	  	 = 0.1
end

Insertion Info

Particle insertion information is defined in this section. This information includes insertion method, inserted number of particles at each insertion step, insertion frequency, dimensions of the insertion box, approximate initial distance between particles at insertion steps, and information about the randomness of initial positions of particles. If the insertion method is uniform, the particles are uniformly (without randomness in initial location) inserted. On the other hand, in non_uniform insertion, initial locations of particles in x and y directions are chosen randomly. Using non_uniform insertion, insertion random number range and insertion random number seed are also required in the parameter handler file.

💡 Insertion in Lethe-DEM starts inserting particles from type 0 and proceeds to the next type when all the particles from the previous type are inserted.

💡 If the insertion box is not large enough to insert inserted number of particles at each time step, Lethe-DEM prints a warning to inform the user about the number of particles that can fit in the insertion box.

💡 insertion distance threshold determines the initial distance between the particles in the insertion phase. As a result, it must be larger than 1 to avoid any initial collision between the inserted particles. It should also be compatible with the random number range; especially if the random number range is large, a large value should be defined for insertion distance threshold. Generally, we recommend users to use a value in the range of 1.3-2 (depending on the value of random number range) for insertion distance threshold.

subsection insertion info
  # Insertion method
  # Choices are uniform|non_uniform
  set insertion method				        = non_uniform

  # Number of inserted particles at each insertion step. This value may change automatically if the insertion box is not adequately large to handle all the inserted particles
  set inserted number of particles at each time step    = 100

  # Insertion frequency
  set insertion frequency            		 	= 20000

  # Insertion box dimensions (xmin, xmax, ymin, ymax, zmin, zmax)
  set insertion box minimum x            	 	= -0.05
  set insertion box minimum y            	        = -0.05
  set insertion box minimum z            	        = -0.03
  set insertion box maximum x            	        = 0.05
  set insertion box maximum y           	 	= 0.05
  set insertion box maximum z            	        = 0.07

  # This value controls the distance between the center of inserted particles (the distance is: [distance threshold] * [diameter of particles]). The distance is modified by a random number if non_uniform insertion is chosen
  set insertion distance threshold			= 2

  # Random number range and seed for non_uniform insertion
  set insertion random number range			= 0.75
  set insertion random number seed			= 19
end

Floating Walls

In this subsection, the information on floating walls is defined. First of all, the total number of floating walls is specified. Then for each floating wall, we should specify its normal vector, a point on the wall, start and end times.

subsection floating walls
  # Total number of floating walls
    set number of floating walls	                = 1
    	subsection wall 0

  # Defining a point on the floating wall
		subsection point on wall
			set x				= 0
			set y				= 0
			set z				= -0.06
		end

  # Defining normal vector of the floating wall
		subsection normal vector
			set nx				= 0	
			set ny				= 0.5
			set nz				= 1
		end

  # Starting time of the floating wall
		set start time				= 0

  # Ending time of the floating wall
		set end time				= 0.2
    	end
end

Mesh

This subsection provides information on the simulation geometry and meshing. The simulation geometry shape and size, the number of refinements and other required information can be provided here. It should be mentioned that meshes from gmsh can also be defined in this section by setting the type = gmsh.

subsection mesh
  # Mesh type
  # Choices are dealii|gmsh
  set type                  				= dealii

  # Grid type in dealii
  set grid type      	     				= hyper_ball

  # Dealii grid arguments
  set grid arguments       				= 0.0, 0.0, 0.0 : 0.1 : true

  # Grid initial refinement
  set initial refinement   				= 3
end

Boundary Conditions

In this subsection, the boundary conditions of the DEM simulation are defined. First of all, the number of boundary conditions is specified. Then for each boundary condition, its information is defined. There are four boundary types: fixed_wall, outlet, rotational (around the center), and translational. For rotational motion, rotational speed and rotational vector are required, while for translational motion, the speed should be defined in each direction.

💡 The default boundary condition type is fixed_wall. If the outlet condition is chosen for a boundary, particles can leave the simulation domain via this outlet.

subsection DEM boundary conditions
  # Total number of boundary motion
  set number of boundary conditions         		= 1

    # For each motion, we need a separate subsection
    subsection boundary condition 0

        # ID of boundary
	set boundary id					= 0

        # Boundary type
        # Choices are fixed_wall|outlet|rotational|translational
        set type              				= rotational

        # Rotational speed magnitude
	set rotational speed				= 2.5

        # Rotational vector
	set rotational vector x				= 1
	set rotational vector y				= 0
	set rotational vector z				= 0
    end

  # OR for translational motion
    subsection boundary condition 0

  # ID of moving boundary
	set boundary id	 				= 4

  # Motion type
        set type              				= translational

  # Speed in each direction
	set speed x					= 0.15
	set speed y					= 0
	set speed z					= 0
    end
end

Boundary Forces

In this subsection, boundary forces are calculated and thanks to that, specific motion is possible.

Calculation makes the calculation possible.
verbosity shows the forces (Force and Torque) in the screen with the output frequency
filename is the filename where the outputs will be written
center of mass corresponds to the initial position of the triangulation's center of mass
moment of inertia is the same but it won't change over time
triangulation mass corresponds to as the name suggests to the triangulation mass

triangulation motion

In this subsection, you make the triangulation motion thanks to the particles forces possible!

Enable is the parameter that allows it
Initial translational and rotational velocity sets the triangulation's initial velocity.
case permits the triangulation to move with specific conditions. There is just the cylinder case now.
slope's inclination is only used for the case cylinder right now and defines as the name suggests, a slope's inclination.

subsection boundary forces
    set calculation               	                                        = true
    set verbosity                                                               = quiet
    set filename                                            	     		= force_and_torque
    set output frequency                                                        = 1
    set center of mass                                                          = 0,0,0
    set moment of inertia                                                       = 0.00239,0,0
    set triangulation mass                                                      = 0.545

    subsection triangulation motion
        set enable                                                                  = true
        set initial translational velocity                                          = 0,0,0
        set initial rotational velocity                                             = 0,0,0
        set case                                                                    = cylinder
        set slope's inclination                                                     = 4
    end
end

Provide feedback

Saved searches

Use saved searches to filter your results more quickly