A modified discrete element method for concave granular materials based on energy-conserving contact model

Ting Qiao, Ji Li, Shunying Ji

State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, China

KeyWords:Discrete element method Concave shapes Energy conservation Contact volume-based contact model Volumetric mesh representation

ABSTRACT The development of a general discrete element method for irregularly shaped particles is the core issue of the simulation of the dynamic behavior of granular materials.The general energy-conserving contact theory is used to establish a universal discrete element method suitable for particle contact of arbitrary shape.In this study, three dimentional (3D) modeling and scanning techniques are used to obtain a trian- gular mesh representation of the true particles containing typical concave particles.The contact volume- based energy-conserving model is used to realize the contact detection between irregularly shaped par- ticles, and the contact force model is refined and modified to describe the contact under real conditions.The inelastic collision processes between the particles and boundaries are simulated to verify the ro- bustness of the modified contact force model and its applicability to the multi-point contact mode.In addition, the packing process and the flow process of a large number of irregular particles are simu- lated with the modified discrete element method (DEM) to illustrate the applicability of the method of complex problems.

The discrete element method (DEM) takes each particle or each particle cluster of granular material as an independent discrete el- ement, which is one of the typical methods to describe the macro- scopic mechanical behavior of granular material in the complex in- dustrial environment [1-3] .The discrete elements are connected through contact, so the mechanical modeling of the contact be- tween particles is the core issue of particle dynamics.In the de- velopment and application of the discrete element method, how to establish a suitable mathematical model for complex particle shapes and at the same time reasonably describe the contact be- tween complex-shaped particles has become one of the hot issues in computational granular mechanics [4] .

In the current representation methods for particles with com- plex shapes, single particle methods mainly include polygonal and polyhedral elements [5] , superquadric elements [ 6 , 7 ] and heart- shaped elements [8] represented by continuous functions, spheri- cal harmonic elements [9] represented by a discrete function, and complex particles elements represented by level sets [10] .The composite particle method mainly includes a combined geometric element model and bonding model [6] .These various geometrical expressions enable the real particle morphology in nature and in- dustry to be described relatively accurately.

In recent years, with the continuous development of complex particles, the contact model of non-spherical particles has gradu- ally been enriched and improved.The mainstream non-spherical particle contact model is based on the traditional nonlinear visco- elastic contact model of the sphere [4] , and further considers the local radius of curvature or overlap volume at the contact point [11] .Different contact theories require different contact geome- try solution goals.Therefore, based on different contact models, corresponding contact detection algorithms can be established for particles of different shapes.The common feature of mainstream non-spherical particle contact models is to decouple and simplify physical and geometric information.The geometric solution targets of contact detection are often: contact point, contact normal, and contact overlap [12] .However, it’s difficult to solve the contact be- tween concave particles of the problem of multi-point contact in such methods.For particle elements of any shape, the diversity of contact modes not only increases the complexity of the calculation procedure but also increases the instability of the solution system.Recently, Feng developed the universal contact theory based on the conservation of energy [ 13 , 14 ].This theory considers the physical and geometric information of the contact at the same time by con- structing a reasonable contact energy function to determine the contact force and moment in any contact mode.

In this letter, the discrete element method for arbitrary irreg- ularly shaped particles based on the energy conservation method is introduced.The contact force model is refined and modified and is used to describe the multi-point contact between concave parti- cles.Granular flows with concave particles are simulated with the modified DEM.

The energy conservation contact theory is an efficient contact detection method for concave particles.In the general energy- conserving contact theory [13] , if the particleA(x,θ)in con- tact with particleB(x＇,θ＇), there exists a scalar potential func- tionω(x,θ,x＇,θ＇)which defines a global function in the contact areaΩAas shown in Fig 1 .And the(F,M)is the negative gradient of the functionω(x,θ,x＇,θ＇), where(F,M)need to be a poten- tial vector field.The normal contact force and the corresponding normal contact moment acting on particleAfrom particleBare determined as:

The contact normalncan be deduced as the unit vector ofFn, and can written as:

The contact pointxc, at where the contact moment equal to zero, is given as:

whereλis a free parameter.

If a reasonable form of the contact energy functionω(x,θ,x＇,θ＇)has been chosen, the full normal contact features can be determined.Therefore, the form of the contact energy functionω(x,θ,x＇,θ＇)become the key of the general energy-conserving contact theory.To reasonably give a specified contact energy function form, the contact volume-based energy-conserving model is proposed by the contact energy function as [15] :

whereVc is the volume of the contact regionΩAbetween two con- tact particles.Then, the normal force is determined by:

whereGnis a vector related to the geometric characteristics of the contact area, which can be understood as a static moment in the form of a vector about the contact surface.In this model, only three geometric quantitiesSn,GnandVc(orω’(Vc)) need to be calculated.When using the linear contact energy function of the contact volume, contact energy function is determined by

whereknis the penalty coefficient.For such contact energy func- tion form, only two geometric quantitiesSnandGnneed to be calculated.An easy-to-implement scheme was developed by Feng [15] to efficiently calculate geometric feature quantities, which can be expressed as:

For complex particle morphology, it is still difficult to solve the geometric feature quantity using the above-mentioned form.A simplified idea is to discretize the complex particle morphol- ogy by using a triangular mesh [ 16 , 17 ].For two contacting triangu- lar mesh discretized particles, their intersecting boundary changes from a closed curve to a closed polygonal line as shown in Fig 2 .The line segments that make up the closed polygonal line can be obtained by intersecting the triangles.And it is easy to get themcorner pointsxi(i= 1, 2, …, m) of the closed polygon line.There- fore, the integral expressions for calculating two geometric feature quantities can be simplified into the summation form, and can be written as:

whereΔxi=xi+1 −xi.Due to the additivity of the calculation for- mat of geometric features, the multi-point contact of concave par- ticles is realized.

The normal contact force solved by the energy conservation contact theory is only the elastic part of the contact.The selection scheme of the penalty coefficientknin the contact energy function given by Eq.(9) needs to be further determined.In order to de- scribe the contact under real conditions, the viscous and frictional effects of the contact need to be additionally considered.Therefore, it is necessary to refine and modify the contact force model based on the energy conservation contact theory.

Firstly, consider the contact energy of the contact area as the elastic strain energy of the contact area,knis the elastic strain en- ergy density, which is expressed as:

whereEis Young’s modulus of the material, andεis the char- acteristic strain that depicts the relative deformation represented by the overlap of the discrete element method based on the soft sphere model.

By comparing with the classical Hertz contact model, the spe- cific form of characteristic strainεcan be obtained and its value range can be reasonably determined.Take the case where the sphere element with radius R is in contact with the plane as an ex- ample, when the sphere element and the plane are in contact with the overlapδn.The directions of the contact force and the con- tact overlap are along the normal direction of the plane.Then the amount of normal elastic contact force determined by the Hertz contact model on a spherical particle is:

and the amount of normal elastic contact force determined by the contact volume-based energy-conserving model (when using the linear contact energy function of the contact volume) is

Comparing Eqs.(16) and (17) , we have

Noteα=δn/R, and,βdescribes the relationship between relative deformationαand characteristic strainε.The form ofβshows that the characteristic strains ex- pressed when the mass points in the overlapping area of the par- ticles are in contact at different overlapping depths are different.This is because the pressure of the contact surface increases with the increase of contact deformation when the actual elastic body is in contact.Considering that the contact surface is very small com- pared to the particle surface, it is assumed that during the contact process, the change ofβwith the relative deformationαis con- stant.Therefore, it is necessary to select an appropriate value ofβwithin a reasonable range.It is assumed that when the con- tact stiffness in the three-dimensional space is 106～108N/m, the maximum overlap is 0.1% ～0.5% of the particle radius, soαmax= 0.1% ～0.5% ,βcan take a value in the range of 0.005 ～0.01, and the corresponding characteristic strainεcan take a value between 0.1 ～0.2 when Poission’s ratioν= 0.3 .Hereεis selected as 0.1.

Then considering a parallel spring-damping contact model, the contact point is not affected by viscous damping, thus viscous damping is taken into account in the contact force model directly.When the normal elastic contact force is determined by the gen-

eral energy-conserving contact theory as:

the normal linear damping forceis considered as:

where,Cn,vn,ijandmi(i= 1, 2) respectively represent the nor- mal damping coefficient, the normal relative velocity between two contacting objects, and the mass of the contacting objecti.

In addition, the tangential forceFs between two contacting ob- jects includes elastic forceand viscous force, respectively ex- pressed as:

whereμs,Cs ,vs,ij, andδsrespectively represent the sliding friction coefficient, tangential viscosity coefficient, tangential relative ve- locity between elements, and tangential overlap considered the Coulomb law of friction.The rolling friction moment(Mr)is used to hinder the relative rotation between the elements, which is ex- pressed as:

whereμris the rolling friction coefficient.is obtained by=ωij/|ωij| , whereωijis the relative rotational velocity.

According to the introduced contact detection scheme and modified contact force model, the discrete element method for concave particles only leaves the step of how to reasonably rep- resent the shape of the particles in the form of triangular meshes.Using scanning technology and three dimentional (3D) modeling technology, it is easy to get very rich irregularly shaped particles.For instance, Fig.3 a and 3 b respectively show a ball and a con- cave block reconstructed by 3D modeling technology, Fig.3 c and 3 d show two real ballast particle shapes reconstructed by scan- ning technology.It is easy to observe that particle shapes (b)-(d) are non-convex particles, so the multi-point contact mode is one of the main contact modes.For the boundary also divided by the tri- angular mesh, the boundary-particle and particle-particle contact detection schemes are the same.

Fig. 1. Schematic diagram of the contact geometry when two particles are in con- tact (The area indicated by the red dashed line represents the contact areaΩ A).

Fig. 2. The calculation scheme of vector areaS n.(a) The boundary intersectionΓin the form of a closed curve; (b) The boundary intersectionΓin the form of a closed polygonal line.

For further verifying the concave particle discrete element method, take the boundary-particle collision process as an exam- ple to illustrate the rationality of the contact model and to show the stability of multi-point contact more clearly.Thus, these four shaped elements are used to simulate the inelastic collision pro- cess between particle and boundary.Particles of different shapes with the same equivalent diameterd= 0.1 m are placed at a dis- tance of 0.1 m from the rigid wall.Among them, the simulation of spherical particle Fig.3 a is used as the reference group.The rele- vant discrete element calculation parameters are listed in Table 1 .

Table 1 DEM simulation parameters.

Fig. 3. Four shaped elements and their triangular mesh representation.

Fig. 4. Inelastic collisions between particles and boundary in a gravity environment.

Figure 4 shows the inelastic collision process of four kinds of particles with boundaries in a gravity environment.The particles finally stably fall on the boundaries due to the existence of gravity and energy dissipation.Figure 5 shows the energy evolution dur- ing the inelastic collision of four kinds of particles with boundaries.Among them, the black line represents the total kinetic energy of the particle, including translational kinetic energy and rotational kinetic energy, the red line represents the gravitational potential energy of the particle (the height of the plane boundary is the zero potential energy surface), and the blue line represents the total energy of the particle which is the sum of the particle’s kinetic energy and gravitational potential energy..The results show that the kinetic energy of the four types of particles eventually tends to zero, and the total energy tends to a stable value, therefore, this contact model is stable.It can be observed that the concave ge- ometry (Fig.4 b) of the typical multi-point contact mode steadily finally stays on the boundary, and there is no sudden change in the energy evolution process (Fig.4 b).

For further evaluating the reliability of the above method, the packing processes of four particle shapes shown in Fig.3 are sim- ulated.800 particles with random orientations and positions are packed with gravity in each simulation.The length and width of the rectangular container are 0.8 m, and the height is 1.5 m.The remaining calculation parameters are listed in Table 1 .The con- tact model considering damping and friction which is adopted in the simulations of inelastic collision between particles and bound- ary is also used here.Figure 6 shows the falling and packing pro- cess of four shapes of particles.The color of the particles indicates the distance between the particles and the bottom of the container in the initial state.Finally, all the particles remain stationary and form stable packing systems.According to Fig.7 , irregular ballast- shaped particles Figs.3 c and 3 d have a higher packing fraction, while concave geometry Fig.3 b has a lower packing fraction.This is because the interlocking between concave geometries causes the local arching structure and more voids.The model of the granular column collapse is also established as shown in Fig.8 .The granu- lar column is filled with particles of 30 shapes, and the same color in Fig.8 represents the same particle shape.Figure 9 shows the flow process of granular accumulation with 30 particle shapes ob- tained by scanning.The color in Fig.9 represents the velocity of the particles.The remaining variables are consistent with Table 1 .Figures.6-9 show the applicability and rationality of the discrete element method for irregular particles based on the energy conser- vation method for describing the process of particle accumulation and particle flow.

Fig. 5. Energy evolutions in the process of inelastic collision between particles and boundary.

Fig. 6. Falling and packing process of four shapes of particles.

Fig. 7. Comparison of packing fractions of different particle shapes.

Fig. 8. The model of the granular column collapse.

Fig. 9. Flow process of granular accumulation with multiple particle shapes.

This letter introduced the application of energy-conserving con- tact theory in the discrete element method for irregularly shaped particles and modified the contact force of the contact volume- based energy-conserving model.The rationality of the modified contact force model and the stability of multi-point contact were illustrated by the simulations of the boundary-particle collision process.In addition, the simulation of the packing and flow process of a large number of particles shows that the method is effective and applicable in describing the contact between a large number of particles and characterizing the behavior of particle clusters.All the simulations are implemented under the GPU parallel frame- work, and the GPU processor of NVIDIA Quadro GV100 is used for calculations.The computational efficiency is significantly affected by the triangular mesh density of irregular particles.For irregu- lar particle shapes with a triangular mesh number of about 400, the calculation of the accumulation process of 800 particles takes about 3.5 hours.Therefore, while ensuring the calculation accuracy, the calculation efficiency of this method is also acceptable.

Declaration of Competing Interest

We declare that we have no known competing financial inter- ests or personal relationships that could have appeared to influ- ence the work reported in this paper.

Acknowledgments

This study is financially supported by the National Key Re- search and Development Program of China (2018YFA0605902) and the National Natural Science Foundation of China (42176241 and 11872136).We also appreciate the discussion on the energy- conserving contact theory with Professor Yun-tian Feng of Swansea University, UK.