Discrete element method modification for the transition to a linearly elastic body model

The article addresses the issue of description of lateral pressure in the discrete element method. The authors formulate a one-dimensional problem on elastic deformation of a particle under the action of the uncompensated normal force and torque. The elastic change in the profile of the particle is calculated and is then included in problem solving by DEM and a specifically developed algorithm. The calculations of compression applied to a regular packing of discrete particles with regard to lateral pressure are exemplified. The developed procedure is the DEM modification enabling transition to modeling a linearly elastic body.


Introduction
Currently, the calculations of deformation and failure of solids use mostly the discrete element method (DEM), or the methods based on continuum models (finite elements, boundary integral equations, gridless methods, etc). In the discrete element method, a real medium is represented by the packing of discrete particles. Then, the laws of interaction between particles are assigned, and the problem of deformation of a body is replaced by the problem of deformation of the set packing of discrete particles. Such approach is faced with no fundamental difficulties in description of physical nonlinearity or even any large deformation. Moreover, it readily takes into account dynamic effects. On the other hand, you know, when you win something, you lose something. In the case discussed, we lose in the degree of adequacy of the obtained results. Indeed, after setting a synthesis algorithm for a stochastic particle packing, as well as properties of the particles and their interaction laws, and then having solved the problem of deformation of the packing, it will be difficult to answer which macroproperties of the packing the solution of the problem corresponds to [1][2][3][4]. Furthermore, a problem connected with the consistency arises. The classical approach to deformation problems involves the models of elasticity, plasticity, creep, etc. Naturally, when one and the same object is analyzed using different methods, the results should be consistent. For example, transition to a linearly elastic body should be fulfilled in the limit case. The method of discrete elements [5][6][7][8] dissatisfies this condition since the forces at the contact of particles are connected with the relative overlaps and displacements of particles only at that contact. Therefore, it is simple to set such a regular packing of particles that no lateral pressure arises in the packing under the uniaxial compression. It is noteworthy that the attempts to modify DEM with a view to describing the lateral pressure phenomenon are currently undertaken [9,10].
This article describes a DEM modification that allows avoiding such paradoxes. To this effect, it is required to set more complex conditions at the contacts of particles. Now, the forces on a contact of particles depend on the overlaps at the given contact and at the other contacts of these particles.

Particle deformation under uncompensated forces
Let us consider dynamics of a particle packing. It is intended to take into account deformation of particles. In principle, this problem can be formulated sufficiently strictly. It is required to write dynamic equations inside each particle and to take into account contact conditions between the particles and between the particles and a loading device. At the same time, it is necessary to include conditions at the free boundary of particles and all initial conditions. Evidently, such formulation makes the problem immeasurable both in terms of solving and analyzing results even if reachablethere is too much information.
Thus, it is of no sense to deal with a strict problem formulation. Some simplifications should be made. For this reason, we neglect wave dynamics inside particles and account for the dynamics of a particle as a whole. This is a major premiss.
Let us discuss the simplest and sufficiently adequate variant. In the discrete element method, as a rule, we deal with an unbalanced system of forces applied to a particle. The particle accelerates, interacts with the neighbor particles and, finally, attains equilibrium. Based on that, first, we analyze a problem on movement of a particle under an unbalanced force. The particle is assumed linearly elastic and rounded.
We limit ourselves to planar deformation. It is convenient to demonstrate the key ideas using a particle shaped as a rod with a length L (figure 1): Here, x x u ,  -stresses and displacements, respectively;  -density, E -elasticity modulus. Let the right-hand end L x  of the rod be free of stress and the left-hand end 0  x be subject to a constant compression stress: To put it simpler, the rod is pushed rightwards in the line of the Ox axis. Let us also assume that the wave process quickly attenuates in the rod (the causes of which are beyond the scope of the discussion). At the steady-state stage, the rod has a constant acceleration a given by: This follows from the equation of dynamics: S -cross section area of the rod, and m -mass of the rod. Consequently, x C x C -arbitrary functions. Then, the first equation in the system (1) with the constant acceleration a , 0 or , with regard to the boundary conditions (2), yields a solution: Regarding the second equation in the system (1): its solution, considering (4), is given by: Based on the problem discussed, some conclusions are drawn: 1. The initially dynamic problem can be reduced to a static problem with the mass forces equal to the inertia forces with an opposite sign (Eq. (5)) -d'Alembert principle.
2. The force application point, i.e. , 0  x moves by the law (10). 3. Including lateral distortion in the solution of (10), we determine the shape of the particle after its deformation.
In this case, the numerical solution of a dynamic problem on motion and deformation of an elastic particle under an uncompensated force follows the steps below: 1. Presetting the value of the uncompensated force, using the formula (3) and finding the acceleration a.
2. Setting the boundary condition which "automatically" provides the solution with the stress ) 0 ( Regarding the issue of moments and rotations, we start from the simplest problem as well. Let a particle be shaped as a cylinder and rotate around an immobile center. The angle of rotation is denoted by ) (t  (figure 2). It is obvious that the element of the medium has the coordinates given by: whence it follows that: Considering that M -uncompensated torque applied at the boundary of a particle and J -inertia moment, it can be written that: . , The forces (14) contribute to the particle deformation. These constructs allow formulating a general algorithm of constitutive equations for an elastic particle in the method of discrete elements: 1. At a current time t the position and configuration of a particle are known. The velocity of the particle center of inertia and the rotational speed are known, too. These values are the initial conditions of the problem.
2. At the time t at the known boundary points, the particle is in contact with the other particles. By the values of overlaps, the forces applied to the particle by the neighbor particles are determined:  3. The static problem solution is connected with the initial position of the particle and acts as the constitutive relations later on. As a result, we find: Moreover, the shape of the particle is determined:

The time increment  is selected. From the equations of motion
and the initial conditions, new average values of the displacements and rotation are determined: The resultant displacements of points at the particle contact are determined: Accordingly, the new position of the particle at the time is calculated. Then, the new system of forces is determined, etc.

Example of calculation of an elastic particle profile
Below in this article an example of calculating a plane static problem with the preset bulk forces is given. Let a round particle at a boundary point A be subjected to an unbalanced surface force F and an unbalanced torque M ( figure 4). In accord with the above algorithm, we replace the surface forces by the bulk forces, i.e. analyze the equilibrium equation (16). The boundary conditions are taken as: Such conditions ensure the required preset load at the point .
A Figure 4. Formulation of the static problem.
The solution of (16), (21) was implemented in the elastic formulation using the finite element method. The elastic characteristics of the particle were chosen as the average values of the elasticity modulus Regarding DEM, the normal forces and, accordingly, the shearing forces due to the particle interactions are applied through a finite size area (governed by the overlap of the particles). For this reason, in the numerical calculation, the point A was replaced by the finite area with a size characterized by the polar angle  ( figure 4). The calculations assumed The problem solution was the displacement profile of the particle boundary. Figure 5 figure 6a and figure 6b, respectively). It is seen in the figures that the particle exposed to the unbalanced load undergoes the change in the profile of its boundary. Compression of the particle induces lateral pressure (refer to figure 5b). As a result, the obtained field can be used in the method of discrete element by the algorithm described in Section 2. The profile of an elastic particle during deformation is only calculated once. Later on, the calculated profile is a parameter of the problem solving by DEM.

Application of the modified discrete element method
In the framework of DEM, we present an example of calculating compression of a flat regular packing of particles centers of which lie at the points of a square lattice in the initial state ( figure 7). It is obvious that in the problem of compression of a solid elastic body with a square cross section, the coefficient of lateral pressure will coincide with Poisson's ratio.
In the classical method of discrete elements [5][6][7][8], movement of particles is calculated using equations of translational and rotational motion (18). The computation of forces and moments of forces at the contacts of particles and at the contacts of particles and loading plates involves Kelvin-Voigt viscoelastic model with regard to the normal and shear components. The normal force is found using the Hertz law from the contact mechanics [11], and the shear force is determined based on the Mindlin-Deresiewicz concept [12,13]. In order to include lateral pressure induced by elastic interaction between particles, we take the algorithm from Section 2 and the calculated elastic change in the particle profile from Section 3. It appears sufficient to consider the profile change only under the action of the uncompensated normal force without regard to the torque: Let us vary elastic properties of particles by increasing compressive stiffness and ductility orthogonally to the compression orientation. Assume  . Given such parameters, an elastic problem is again calculated using the finite element method (Section 3) and a new profile of a particle is determined. Inclusion of the latter in the problem In this way, by varying elastic characteristics of a particle and sizes of contact areas between particles, it is possible to obtain sufficiency realistic values of lateral pressure. Consequently, the presented procedure of numerical calculation is the discrete element method modification which permits transition to linear elasticity.

Conclusions
The presented approach allows taking into account elastic properties of particles and enables description of lateral pressure using the method of discrete elements.
The proposed procedure is the modification of the discrete element method for the transition to modeling a linearly elastic body.
This study was supported by the Russian Science Foundation, Project No. 16-17-10121.