-
Notifications
You must be signed in to change notification settings - Fork 60
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.
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 Motion
This subsection contains the general information of the simulation, including time-step, end time, output path, log and output frequency. It is the most commonly modified section for a simulation. 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
This subsection deals with simulation timing information. If it is set to end, the simulation timing information will be printed at the end of simulation.
subsection timer
set type = end
end
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 simulation. Then these information can be compared with predefined values to ensure reproducibility of the simulations.
subsection test
# DEM test
# Choices are true|false
set type = false
end
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 mode), 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. Choosing once repartition method, 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 its non-uniformity exceeds load balance threshold, it calls load-balancing.
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 pp_linear|pp_nonlinear
set particle particle contact force method = pp_nonlinear
# Particle-wall contact force model
# Choices are pw_linear|pw_nonlinear
set particle wall contact force method = pw_nonlinear
# Integration method
# Choices are euler|velocity_verlet
set integration method = velocity_verlet
end
In this subsection, 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.
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 = 2000
# Young's modulus of particle
set young modulus particle = 1000000
# Poisson ratio of particle
set poisson ratio particle = 0.3
# Coefficient of restitution of particle
set restitution coefficient particle = 0.95
# Coefficient of friction of particle
set friction coefficient particle = 0.05
# Coefficient of rolling friction of particle
set rolling friction particle = 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
Particle insertion information are defined in this section. These information include 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 randomness of initial positions of particles. If 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.
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 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
In this subsection, the information of floating walls are defined. First of all, 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
This subsection provides information of the simulation geometry and meshing. The simulation geometry shape and size, 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
In this subsection, motion of boundaries are defined. First of all, number of boundary motion is specified. Then for each boundary motion, information of motion are defined. There are two boundary motion types: 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.
subsection mesh
# Total number of boundary motion
set number of boundary motion = 1
# For each motion, we need a separate subsection
subsection moving boundary 0
# ID of moving boundary
set boundary id = 4
# Motion type
# Choices are 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 moving boundary 0
# ID of moving boundary
set boundary id = 4
# Motion type
# Choices are rotational|translational
set type = translational
# Speed in each direction
set speed x = 0.15
set speed y = 0
set speed z = 0
end
end
In this subsection, boundary forces are calculated and thanks to that, specific motion is possible.
- Calculation make 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 overtime
- triangulation mass corresponds to as the name suggests to the triangulation mass
In this sub-subection, you make the triangulation motion thanks to the particles forces possible!
- Enable is the parameters that allows it
- Initial translational and rotational velocity sets the triangulation's initial velocity.
- case permits to the triangulation to moves with specifics conditions. There is just 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