En memoria del meu amic Pablo

## Acknowledgements

In the very last lesson of Estructuras III of my undergraduate degree in Civil Engineering, Professor Oñate presented several animations of simulations carried out in CIMNE. As a student, to see that all those bunch of equations of finite elements could be applied in real engineering problems greatly amazed me. Surprisingly, when I expressed to Oñate my interest in the field, he offered me a position to start implementing the Discrete Element Method in the new plataform Kratos. I combined it with the development of my final degree monograph and I moved on undertaking the Master in Numerical Methods in Engineering. Afterwards, I got a scholarship from the programme Doctorats Industrials de la Generalitat de Catalunya which allowed me to develop my Ph.D. thesis in a partnership of the research centre CIMNE and the company CITECHSA.

This work encompasses the result of the 4 years in CIMNE working in the field of Discrete Element Methods under the guidance of Prof. Eugenio Oñate.

First of all, I must thank my advisor Prof. Oñate for his support all over these years. His advice has been not only on the topic of research but also on the development of my career. He gave me great freedom in the decision of the research line and helped me having a rich international experience through my research stay abroad as well as the participation in several conferences. I have learned a lot from Prof. Oñate in many aspects and I consider myself fortunate to have had the opportunity to work with him.

I would like to mention Miguel Ángel Celigueta who has helped me a lot, specially during my first steps in CIMNE by allowing me take part in several ongoing project meetings. From him I learned most of what I know about coding in an efficient and organized manner. Later, J. Maria Carbonell became my second supervisor of the monograph helping me in the developments regarding the coupling with the solid mechanics code as well as in the elaboration of this document. For that and for his support along my monograph I would like to express my gratitude to him.

I would also like to thank Prof. Wriggers for finding me a seat in the Institute of Continuum Mechanics in Hanover where I learnt about contact mechanics and developed part of the monograph under his advisory. My stay in Germany has been very fruitful for my monograph but also for learning German. That is specially Tobias Steiner to thank, my Bürokollege and now my good friend. Danke!

From the CITECHSA side I want to acknowledge Natalia Alonso and María Angeles Viciana and thank the rest of the team as well.

Thanks to the DEM Team members: Salva, Ferran and Guillermo for the work done together and the intense discussions on the DEM. Thanks to Joaquín Irázabal who has become my closer collaborator; with him I share an article, a lot of developments and of course good moments in congresses during these years. I don't want to miss mentioning Jordi and Pablo who are exemplar engineers to me and Charlie who solved memory errors in my code uncountable times. Also Ignasi, Roberto, Kike, Pooyan, Riccardo, Antonia, Miguel, Abel, Anna, Adria, Javier Mora, Sonia S., Feng Chun, María Jesús and the rest of the CIMNE family. A little part of them is somehow in this monograph. Thank you!

Last but not least to my friends from Navas and specially to my family. Gracies Montse, Marina, Mery, Josep and Celia. 谢谢 Xiaojing for your necessary support and company during the last steps. This work was carried out with financial support from the programme Doctorats Industrials de la Generalitat de Catalunya, Weatherford Ltd. and the BALAMED project (BIA2012-39172) of MINECO, Spain.

## Abstract

This works encompasses a broad review of the basic aspects of the Discrete Element Method for its application to general granular material handling problems with special emphasis on the topics of particle-structure interaction and the modelling of cohesive materials. On the one hand, a special contact detection algorithm has been developed for the case of spherical particles representing the granular media in contact with the finite elements that discretize the surface of rigid structures. The method, named Double Hierarchy Method, improves the existing state of the art in the field by solving the problems that non-smooth contact regions and multi contact situations present. This topic is later extended to the contact with deformable structures by means of a coupled DE-FE method. To do so, a special procedure is described aiming to consistently transfer the contact forces, which are first calculated on the particles, to the nodes of the FE representing the solids or structures. On the other hand, a model developed by Oñate et al. for the modelling of cohesive materials with the DEM is numerically analysed to draw some conclusions about its capabilities and limitations.

In parallel to the theoretical developments, one of the objectives of the monograph is to provide the industrial partner of the doctoral programme, CITECHSA, a computer software called DEMPack (www.cimne.com/dem/) that can apply the coupled DE-FE procedure to real engineering projects. One of the remarkable applications of the developments in the framework of the monograph has been a project with the company Weatherford Ltd. involving the simulation of concrete-like material testing.

The monograph is framed within the first graduation (2012-13) of the Industrial Doctorate programme of the Generalitat de Catalunya. The monograph proposal comes out from the agreement between the company CITECHSA and the research centre CIMNE from the Polytechnical University of Catalonia (UPC).

## Resum

Aquest treball compen una amplia revisió dels aspectes basics del Metode dels Elements Discrets (DEM) per a la seva aplicació generica en problemes que involucren la manipulació i transport de material granular posant emfasi en els temes de la interacció partícula-estructura i la simulació de materials cohesius. Per una banda, s'ha desenvolupat un algoritme especialitzat en la detecció de contactes entre partícules esferiques que representen el medi granular i els elements finits que conformen una malla de superfície en el modelatge d'estructures rígides. El metode, anomenat Double Hierarchy Method, suposa una millora en l'estat de l'art existent en solucionar els problemes que deriven del contacte en regions de transició no suau i en casos amb múltiples contactes. Aquest tema és posteriorment estes al contacte amb estructures deformables per mitja de l'acoblament entre el DEM i el Metode dels Elements Finits (FEM) el qual governa la solució de mecanica de solids en l'estructura. Per a fer-ho, es descriu un procediment pel qual les forces de contacte, que es calculen en les partícules, es transfereixen de forma consistent als nodes que formen part de l'estructura o solid en qüestió. Per altra banda, un model desenvolupat per Oñate et al. per a modelar materials cohesius mitjancant el DEM és analitzat numericament per tal d'extreure conclusions sobre les seves capacitats i limitacions.

En paral${\displaystyle \cdot }$lel als desenvolupaments teorics, un dels objectius de la monografia és proveir al partner industrial del programa doctoral, CITECHSA, d'un software anomenat DEMpack (http://www.cimne.com/dem/) que permeti aplicar l'acoblament DEM-FEM en projectes d'enginyeria reals. Una de les aplicacions remarcables dels desenvolupaments en el marc de la monografia ha estat un projecte per l'empresa Weatherford Ltd. que involucra la simulació de tests en provetes de materials cimentosos tipus formigó.

Aquesta monografia doctoral s'emmarca en la primera promoció (2012-13) del programa de Doctorats Industrials de la Generalitat de Catalunya. La proposta de la monografia prové de l'acord entre l'empresa CITECHSA i el centre de recerca CIMNE de la Universitat Politecnica de Catalunya (UPC).

## 1. Introduction

Truesdell and Noll in the introduction of The Non-Linear Field Theories of Mechanics [1] state:

Whether the continuum approach is justified, in any particular case, is a matter, not for the philosophy or methodology of science, but for the experimental test…

The ones that agree on that statement may also agree that the same applies for the discontinuum approach in which the Discrete Element Method is framed on.

Before the introduction of the Discrete Element Method in the 70's, lot of effort has been placed in developing constitutive models for the macroscopic description of particle flows. However, the continuum based methods fail to predict the special rheology of granular materials which can rapidly change from a solid-like behaviour in zones where the deformation is small and rather homogeneous to a fluid-like behaviour with huge variation and distortion that can be concentrated in narrow areas such as shear bands. Within the DEM this behaviour, which is driven by the collisional and frictional mechanisms of the material, can be simulated at the grain level where each discrete element corresponds to a physical particle. The quality of the results depends then on the accuracy in the representation of the shape of the particles and their interaction.

The DEM is nothing else than Molecular Dynamics with rotational degrees of freedom and contact mechanisms. In its first conception, the method was designed for simulations of dynamic systems of particles where each element is considered to be an independent and non deformable entity which interacts with other particles by the laws of the contact mechanics and moves following the Newton-Euler equations.

The simplicity of the method is in contrast however, with the high computation cost which, in general, has associated to it due to the large number of particles needed in a real simulation and the time scales that have to be resolved. Imagine a hooper discharge problem which may require the computation of millions of particles simulated during tens of minutes when, at the same time, the phenomena that rules the problem lies in reproducing the behaviour of individual particles the interaction of which happens in distances several orders of magnitude smaller than their particle diameter. This implies that the necessary time steps to be used in the simulation have to be smaller than the characteristic contact duration.

In this sense, the implementation of the method using massive parallelization is something of crucial importance. Also the use of simple geometries such as spheres presents a great difference to other more complex geometries such as polyhedra, NURBS, etc. in the detection and characterization of the contacts. That is why still today the use of basic spheres is intensively used.

In many real applications involving granular materials, the interaction with structures and fluids are present. The employment of the FEM to simulate the structures involved in those industrial applications can provide better understanding of the problem and, therefore, could play an important role in the process of design optimization. To that end an efficient coupling of the method with a FEM-based solver for solids is of special interest.

Another field of interest of the application of the DEM is the simulation of material fracturing. The DEM as a discontinuum-based method has attractive features in contrast to continuum-based methods in problems where large deformations and fracture are involved. Many attemps have been done aiming to unify both the modelling of the mechanical behaviour of solid and particulate materials, including the transition from solid phase to particulate phase. Nowadays however, the DEM still presents many drawbacks and lack of reliability in the modelling of solids. Differently from other particle-based methods such as MPM, PFEM or SPH, the DEM shall not be regarded as a discretization method for the solution of PDE.

The interest in the Discrete Element Method has exponentially increased since the publication in 1979 of the first article by Cundall and Strack [2] and is still a hot topic nowadays. This can be seen in fig. 1 where the number of publications related to discrete element procedures from 1979 to 2016 are displayed. They were obtained from Google Scholar with the following keywords in the title of the article: 'Discrete Element Method/Model', or 'Distinct Element Method/Model', or 'Using a DEM' or 'A DEM' or 'With the DEM' or 'DEM Simulation'. This does not include all the publications related to DEM and may introduce other non related articles, however it gives a good image of the tendency of research in the field.

 Figure 1: Number of publications from 1979 to 2016 obtained from Google Scholar with the following keywords in the title of the article: 'Discrete Element Method/Model', or 'Distinct Element Method/Model', or 'Using a DEM' or 'A DEM' or 'With the DEM' or 'DEM Simulation'.

There is a great interest in the application of this method to a wide range of industrial problems.

### 1.1 DE-FE couplings

The term coupled DE-FE or combined DE-FE Method for soil and solid mechanics applications appears in the literature with different meanings and can be quite confusing. The most common ones are grouped here in 5 categories along with an example figure (Fig. 2). Other categories for DE-FE couplings are for instance coupling with fluids, thermal problems, etc.

1. Particle-structure interaction: The two domains are calculated separately and their communication is through contact models. This is the category in which the monograph is mainly focused on. It is developed in Section 4.
2. Two-scale models: These methods solve the problem at two different scales. The micro-macro transition is accomplished employing an overlapping zone to provide a smooth transition between a DE model (micro) and a FE material description (macro). The coupling is achieved by the imposition of kinematic constrains between the two domains. The original idea was presented by Xiao and Belytschko [3] for Molecular Dynamics, Wellmann [4] applied it to granular material while Rojek and Oñate [5] developed it for cohesive materials.
3.  (a) Particle-Structure (b) Two-scale. Taken from: Labra 101 (c) Projection onto a FE mesh (d) Embedded particles. Taken from: Zárate and Oñate 8 (e) Discretized DE. Taken from: Gethin et al. 9 Figure 2: Examples of different techniques that combine FE and DE methods
4. Projection techniques: Coarse-graining, averaging and other projection techniques are used to derive continuum fields out of discrete quantities. To do so, often a reference mesh is required either for the calculation or simply for the representation of the continuum results [6].
5. Embedded DE on FE: This technique consists on embedding (typically spherical) particles in the boundaries of FE models of solids and structures in order to detect and enforce the contact [7]. Recently, this technique has been applied to multi-fracturing in cohesive materials [8]. A FE-based method with failure or crack propagation models is combined with embedded particles that assist the detection and characterization of the contact forces.
6. FE discretization of discrete entities: This category involves methods that use a FE discretization to calculate the deformation of the particles and solve their interaction using a DEM-like technique [9]. A particular case is the so called DEM-Block method [10] which consists on a FE-based Method which elements are connected through breakable spring-like bounds imitating the cohesive DE models.

### 1.2 Objectives

This monograph has been developed in the framework of the first graduation of Doctorats Industrials de la Generalitat de Catalunya (Industrial Doctorates of Catalonia). The objectives defined for this work comprise an agreement between the research line determined by the research centre CIMNE in the Polytechnical University of Catalonia (UPC) and the business objectives of the society CITECHSA which is interested in the exploitation of a DEM-based software in its application to industrial engineering problems. In this regard, the objectives involve research, development of a code and educational and dissemination actions.

On the one hand, the research has to be focused in a deep revision of the state of the art of the Discrete Element Method in order to analyse and select the existing techniques that have to be adapted and implemented for the solution of the problems of interest which are basically three:

• General application of DEM to granular material handling problems
• Particle-structure interaction
• DE models for the simulation of cohesive materials

The research has to be conducted from a general point of view determining the advantages and drawbacks of the existing methods and proposing new developments that can improve the state of the art. The theoretical contributions will be communicated by dissemination actions.

The theoretical research in the above-mentioned topics have to accompanied by its implementation into the open-source code DEMpack (www.cimne.com/dem). The code will be developed with concerns on efficiency and parallelism as it is devised to be employed in real application projects. To that end, several GUIs for specific applications will be developed. This will be done forming part of a larger group of researchers that contribute to the development of the code.

Finally, the developments will be applied in ongoing projects of the research centre.

### 1.3 Organization of this work

The document is structured as follows:

After the introduction and the objectives, chapter two reviews the basic aspects of the Discrete Element Method that will set the basis upon which the developments in the monograph are established. It includes a revision of the most common contact models and integration methods. An assessment on performance, accuracy and stability is given to help choosing the most appropriate integration scheme. The treatment of clusters of spheres for the representation of non-spherical particles and the contact detection are also discussed in detail.

Chapter three is dedicated to the contact detection between spherical Discrete Elements and triangular or planar quadrilateral Finite Elements. The chapter starts with a complete review of the state of the art and follows with a thorough description of the strategy adopted for the global and local detection of contacts. The idea of using an intermediate fast intersection test is introduced and later proved to be efficient within an application example. Regarding the local resolution, the novel Double Hierarchy Method for contact with rigid boundaries is presented. The description of the methods is equipped with algorithm details, validation examples and limitations analysis.

The fourth chapter introduces the DE-FE coupling for the particle-structure interaction problem. After an introduction to the solid mechanics formulation employed, the coupled scheme is presented. The key point lies in the communication of the contact forces, which are calculated by the DE particles, to the nodes of the FEs. The described procedure proposes the distribution of the forces to all the FEs involved based on their area of intersection with the particles. Several examples show that this strategy improves the commonly used direct interpolation approach for the case of contacts with deformable solids or structures. The good functioning of the coupling is assessed by some tests with special attention placed on energy conservation.

The topic of DE modelling of cohesive materials such as concrete or rock is presented in chapter five. It begins with an overview of the state of the art of the methods available for this purpose together with a study of their limitations and capabilities. After, the model developed by Oñate, Santasusana et al. is described along with application examples where the numerical simulations and the laboratory tests are compared.

Chapter six is dedicated to the implementation of the code in the platform Kratos constituting the DEMpack software together with remarks on the efficiency and parallelitzation of the code.

Finally, the last chapter comprises the conclusions and the outlook of the work.

### 1.4 Related publications and dissemination

#### 1.4.1 Papers in scientific journals

• E. Oñate, F. Zárate, J. Miquel, M. Santasusana et al. Computational Particle Mechanics - Springer: Local constitutive model for the Discrete Element Method. Application to geomaterials and concrete.
• M. Santasusana, J. Irazábal, E. Oñate, J.M. Carbonell. Computational Particle Mechanics - Springer: The Double Hierarchy Method. A parallel 3D contact method for the interaction of spherical particles with rigid FE boundaries using the DEM.

#### 1.4.2 Communications in congresses

• M. Santasusana, E. Oñate, M.A. Celigueta, F. Arrufat, K. Valiullin, R. Gandikota. 11th. World Congress on Computational Mechanics (WCCM XI): A parallelized discrete element method for analysis of drill-bit mechanics problems in hard and soft soils.
• C.A. Roig, P. Dadvand, M. Santasusana, E. Oñate. 11th. World Congress on Computational Mechanics (WCCM XI): Minimal surface partitioning for particle-based models.
• M. Santasusana, E. Oñate, J.M. Carbonell, J. Irazábal, P. Wriggers. 4th. International Conference on Computational Contact Mechanics (ICCCM 2015): Combined DE/FE method for the simulation of particle-solid contact using a Cluster-DEM approach.
• E. Oñate, F. Arrufat, M. Santasusana, J. Miquel, M.A. Celigueta. 4th. International Conference on Particle-Based Methods (Particles 2015): A local constitutive model for multifracture analysis of concrete and geomaterials with DEM.
• M. Santasusana, E. Oñate, J.M. Carbonell, J. Irazábal, P. Wriggers. 4th. International Conference on Particle-Based Methods (Particles 2015): A Coupled FEM-DEM procedure for nonlinear analysis of structural interaction with particles.
• M. Santasusana, J. Irazábal, E. Oñate, J.M. Carbonell. 7th. European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS Congress 2016): Contact Methods for the Interaction of Particles with Rigid and Deformable Structures using a coupled DEM-FEM procedure.

### Participation in research projects

• Implementation of the DEM Technology, DEM material models and validation of the DEM Technology. Weatherford International LTD. July 2013 to March 2014.
• Nuevo proceso de voladura mediante técnicas predictives y adaptatives, eficaz y eficiente en la utilización de los recursos y materias primes, minimizando las emisiones. MINECO. January 2014 to June 2016.

## 2. The Discrete Element Method

The Discrete Element Method (DEM) was firstly introduced by Cundall in 1971 [11] for the analysis of the fracture mechanics problems. Afterwards, in 1979, Cundall and Strack [2] applied it to granular dynamics. The DEM in its original conception simulates the mechanical behaviour of a system formed by a set of particles arbitrarily disposed. The method considers the particles to be discrete elements forming part of a higher more complex system. Each discrete element has an independent movement; the overall behaviour of the system is determined by the application of contact laws in the interaction between the particles.

There exist two main approaches, namely the soft and the hard particle approach. The soft particle approach is a time-driven method where particles are allowed to inter-penetrate simulating small deformations due to contact. The elastic, plastic and frictional forces are calculated out of these deformations. The method allows accounting for multiple simultaneous particle contacts. Once the forces are calculated, the motion of the particles is earned from the application of the classical Newton's law of motion which is usually integrated by means of an explicit scheme. The hard-particle approach, on the other hand, is an event-driven method which treats the contacts as instantaneous and binary (no-multi contact). It uses momentum conservation laws and restitution coefficients (inelastic or frictional contacts) to determine the states of particles after a collision. These assumptions are only valid when the interaction time between particles is small compared to the time of free motion. A good review and comparison of the methods can be found in [12]. This thesis is developed using the soft-particle approach.

The DEM, as a particle method, has been used in a very wide range of applications. An important decision to take is to select which is the relation between a discrete element in the simulation and the physical particles or media in the reality. On the one hand, the one-to-one approach tries to assign a discrete element to every particle in the domain. The method describes all the contact and other interaction forces between particles with a model that only depends on the local relations and does not require fitting. On the other hand, a very common approach is to simulate granular matter or other media using discrete elements that represent a higher amount of material than just one particle. This technique, known as coarse-graining or up-scaling [13], represents a completely phenomenological approach which does require the fitting of parameters out of bulk experiments. Both techniques are used to simulate particulate matter that ranges from powder particles (${\textstyle \mu m}$) to the simulation of rock blocks (${\textstyle m}$).

Common applications of the Discrete Element Method are the simulation of granular mater in soil mechanics. A soil can deform as a solid or flow as a fluid depending on its properties and the situation. The use of DEM comes naturally as it can handle both behaviours of the soil and also account for discontinuous and very large deformations [14,15,4]. The DEM adapts also perfectly to the simulation of granular material handling in industrial processes. Some examples of applications are silo flows [16,12], screw-conveyors [17,18], vibrated beds [19,20], ball mill processes [21,22], etc.

Another application which is of special attention in this thesis is the particle-structure interaction problem. This category encompasses, among others, particle-tyre simulations [14,23], shot peening processes [24,25], impacts with flexible barriers [26], soil-structure interaction [4,27], etc. Some examples of applications are presented in section 6.3.

In particle-fluid flow modelling, the difficulty relies on the particulate phase rather than fluid phase. Therefore, a coupled CFD-DEM approach [28] is attractive because of its capability to capture the particle physics. This comprehends a large family of applications which includes gas fluidization, pneumatic conveying flows, particle coating, blast furnace, etc. [16]. Applications in civil and marine engineering are rock avalanches into water reservoirs [29], sediment and bed-load transportation in rivers and sea [30,31,32], etc. A comprehensive literature review on the applications of DEM to the simulation of particulate systems processes can be found in the work published by Zhu et al. [16].

In recent years the DEM has also been object of intense research to the study of multi-fracture and failure of solids involving geomaterials (soils, rocks, concrete), masonry and ceramic materials, among others. Some key developments can be found in [33,34,35,36]. In the cohesive models the contact law can be seen as the formulation of the material model on the microscopic level. Cohesive bonds can be broken, which allows to simulate the fractures in the material and its propagation. The analysis of solid materials within the DEM poses however, a number of difficulties for adequately reproducing the correct constitutive behaviour of the material under linear (elastic) and non-linear conditions (section 5).

### 2.1 Basic steps for DEM

From a computational point of view a basic DEM algorithm consists of three basic steps:

 Figure 3: Basic computational scheme for the DEM

After an initialization step, the time loop starts. First, the neighbouring particles for every discrete element needs to be detected (section 2.2) as well as the contact with rigid boundaries included in the simulation domain (chapter 3). Afterwards, for every contacting pair a the contact model is applied (section 2.5) to determine the forces and torques that have to be added to the rest of actions to be considered on a particle. Finally, given all the forces and the torques, the equations of motion are integrated and the particle's new position is usually calculated by means of an explicit time marching scheme (section 2.6). At this point, new contacts have to be detected and thus, the loop starts again. This sequence repeats over time until the simulation comes to an end.

### 2.2 Contact detection

Due to the method formulation, the definition of appropriate contact laws is fundamental and a fast contact detection is something of significant importance in DEM. Contact status between individual objects, which can be two DE particles or a DE particle and a boundary element (chapter 3), can be calculated from their relative position at the previous time step and it is used for updating the contact forces at the current step. The relative cost of the contact detection over the total computational cost is generally high in DEM simulations. Therefore, the problem of how to recognize all contacts precisely and efficiently has received considerable attention in the literature [37,38].

Traditionally, the contact detection is split into two stages: Global Neighbour Search and Local Contact Resolution. By the application of this split the computational cost can be reduced from ${\textstyle {\mathcal {O}}(N^{2})}$, in an all-to-all check, to ${\textstyle {\mathcal {O}}(N\cdot ln(N))}$.

#### Global Contact Search

It consists on locating the list of potential contact objects for each given target body. There are two main basic schemes: the Grid/Cell based algorithms and the Tree based ones, each of them with numerous available versions in the literature.

 (a) Grid/Cell-based structure (b) Tree-based structure

In the Grid based algorithms [39,40,41] a general rectangular grid is defined dividing in cells the entire domain (figure cell). A simple bounding box (rectangular or spheric) is adopted to circumscribe the discrete elements (of any shape) and is used to check in a approximate way which are the cells that have intersection with it. Those intersecting cells, store in their local lists the particles contained in the bounding boxes. The potential neighbours for every target particle are determined by selecting all the elements stored in the different cells where the bounding box of that target particle has been assigned to.

In the Tree based algorithms [42,43,44,45] each element is represented by a point ${\textstyle p}$ at coordinates ${\textstyle {\boldsymbol {X}}^{p}}$. Starting from a centred one, it splits the domain into two sub-domains. Points that have larger coordinate (${\textstyle {\boldsymbol {X}}_{i}\geq {{\boldsymbol {X}}_{i}}^{p}}$) are placed in one sub-domain while points with smaller coordinates (${\textstyle {\boldsymbol {X}}_{i}<{{\boldsymbol {X}}_{i}}^{p}}$) in the other sub-domain. The method proceeds for next points alternating every time the splitting dimension ${\textstyle i}$ and obtaining a tree structure like the one shown in figure tree. Once the tree is constructed, for every particle the nearest neighbours is determined following the tree in upwind direction.

Han et al [46] compared the most common Global Neighbour Search algorithms (cell-based and tree-based) in simulations with spherical particles. Numerical tests showed better performance for the cell based algorithms (D-Cell [40] and NBS [39]) over the tree-based ones (ADT [42] and SDT [44]), specially for large-scale problems. It should be noted also that the efficiency depends on the cell dimension and, in general, the size distribution can affect the performance. Han et al [46] suggest a cell size of three times the average discrete object size for 2D and five times for 3D problems. It is worth noting that, using these or other efficient algorithms, the cost of the Global Neighbour Search represents typically less than 5 percent of the total computation while the total cost of the search can reach values over 75 percent [14], specially when the search involves non-spherical geometries since it requires, in general, the resolution of a non-linear system of equations (see the case of superquadrics [4,47] or polyhedra [48,49,50]). In this sense, the focus should be placed on the Local Contact Resolution check rather than on optimizing the Global Neighbour Search algorithms.

#### Local Resolution Check

The local contact detection basically consists in determining which of the potential neighbours found during the global search algorithm constitute an actual contact with the target particle and to determine their contact characteristics (point of contact, normal direction, etc.). The case of spheres is trivial (fig. 4), contact exists if the following condition is met:

 ${\displaystyle \lVert {\boldsymbol {C}}_{i}-{\boldsymbol {C}}_{j}\rVert
(1)

and the normal and point of contact can be easily determined as it will be detailed in section 2.4.

 Figure 4: Spherical particles in contact

The problem of contact determination becomes complex and time consuming when other geometries such as superquadrics, polyhedra or NURBS are used to represent the particles or boundaries. A way to improve the efficiency is to take advantage of the temporal coherence. Normally the duration of a contact is encompassed by several calculation time steps and therefore the particle positions will only change a little bit. In this regard, it seems wise to perform the contact detection after several time steps instead of at every time step aiming to reduce the computational cost that it involves. However, if the contacts are not determined when the particles start to collide, the indentations will achieve high values which will lead to inaccurate results and numerical instabilities (section 2.6.4).

A possible solution for this issue is the use of a technique known as Verlet neighbouring lists [51,4]. It consists on using enlarged bounding boxes in the global search so that more remote particles are stored as well. This local Verlet list need no update during several time steps since the particles move only small distances every step. This way it can be assured that no contacts are missed along the simulation and the frequency of the search is reduced. This method is efficient for cases with high dispersion of particles.

In the framework of this dissertation, a basic cell-based algorithm [40] is chosen which has been parallelized using OMP. The geometries used for the particles are only spheres or clusters of spheres and thus the local detection is efficient. The treatment of the contact with FE representing rigid or deformable boundaries is extensively discussed in chapter 3.

### 2.3 Equations of motion

In the basic soft particle DEM approach the translational and rotational motion of particles are defined by the standard equations for the dynamics of rigid bodies. For the special case of spherical particles, these equations can be written as:

 ${\displaystyle {\mbox{m}}\,{\ddot {u}}={F}\,}$
(2)

 ${\displaystyle I\,{\dot {\boldsymbol {\omega }}}={\mathbf {T} }\,}$
(3)

where ${\textstyle {u}}$, ${\textstyle {\dot {u}}}$,${\textstyle {\ddot {u}}}$ are respectively the particle centroid displacement, its first and second derivative in a fixed coordinate system ${\textstyle X}$, m is the particle mass, ${\textstyle I}$ the inertia tensor, ${\textstyle \omega }$ is the angular velocity and ${\textstyle {\dot {\omega }}}$ the angular acceleration.

The forces {F} and the torques {T} to be considered at the equations of motion (eq.2 and eq. 3) are computed as the sum of:

1. all forces {F}^{ext} and torques {T}^{ext} applied to the particle due to external loads.
2. all the contact interactions with neighbouring spheres and boundary finite elements {F}^{ij}, j=1,\cdots , n^c, where i is the index of the element in consideration and j the neighbour index of the entities (particles or finite elements) being in contact with it.
3. all forces {F}^{damp} and torques {T}^{damp} resulting from external damping.

This can be expressed for every particle \end{array}[/itex]i as:

 ${\displaystyle {\mathbf {T} }_{i}={\mathbf {T} }_{i}^{\textrm {ext}}+\sum \limits _{j=1}^{n^{c}}{\mathbf {r} }_{c}^{ij}\times {\mathbf {F} }^{ij}+{\mathbf {T} }_{i}^{\textrm {damp}}}$
(5)

where ${\displaystyle {r}_{c}^{ij}}$ is the vector connecting the centre of mass of the ${\displaystyle i}$-th particle with the contact point ${\displaystyle {\boldsymbol {Pc}}^{ij}}$ with the ${\displaystyle j}$-th particle (eq. 8). ${\displaystyle {F}^{ij}}$ and ${\displaystyle {F}^{ji}}$ satisfy ${\displaystyle ({F}^{ij}=-{F}^{ji})}$. Fig. 5 shows contact forces between two spherical particles.

The rotational movement equation (3) is a simplified version of the Euler equations coming from the fact that a sphere has constant coefficients for its three principal inertia axes which are independent of the frame. The complete equations can be found in section 2.7 where the case of generic particle shapes is discussed.

### 2.4 Contact kinematics

The forces and torques that develop from a contact event are derived from the contact kinematics at the point of contact ${\displaystyle {Pc}^{ij}}$. The local reference frame in the contact point is defined by a normal ${\displaystyle {n}^{ij}}$ and a tangential ${\displaystyle {t}^{ij}}$ unit vectors as shown in figure 5.
 (a) Contact between two particles (b) Contact force decomposition Figure 5: Kinematics of the contact between two particles

The normal is defined along the line connecting the centres of the two particles and directed outwards from particle ${\displaystyle i}$.

 ${\displaystyle {n}^{ij}={\frac {{\boldsymbol {C}}_{j}-{\boldsymbol {C}}_{i}}{\lVert {\boldsymbol {C}}_{j}-{\boldsymbol {C}}_{i}\rVert }}}$
(6)

The indentation or inter-penetration is calculated as:

 ${\displaystyle \delta _{n}=R_{i}+R_{j}-({\boldsymbol {C}}_{j}-{\boldsymbol {C}}_{i})\cdot {n}^{ij}}$
(7)

where ${\displaystyle {\boldsymbol {C}}_{j}}$, ${\displaystyle {\boldsymbol {C}}_{i}}$ are the centre of the particles and ${\displaystyle R_{i}}$, ${\displaystyle R_{j}}$ their respective radius.

The vectors from the centre of particles to the contact point ${\displaystyle {\boldsymbol {r}}_{c}^{ij}}$ and ${\displaystyle {\boldsymbol {r}}_{c}^{ji}}$ are in general dependent on the contact model. they should take into account the contribution of each particle to the equivalent stiffness of the system. Eq. 8 describes the simple case of two linear springs with different Young's modulus set in serial:

 ${\displaystyle {\begin{array}{l}{\boldsymbol {r}}_{c}^{ij}=\left(R_{i}+{\frac {E_{j}}{E_{i}+E_{j}}}\,\delta _{n}\right){n}^{ij}\\\end{array}}}$
(8)

The position of the contact point can then be determined from any of the particles:

 ${\displaystyle {\begin{array}{l}{\boldsymbol {Pc}}^{ij}={\boldsymbol {C}}_{i}+{\boldsymbol {r}}_{c}^{ij}={\boldsymbol {C}}_{j}+{\boldsymbol {r}}_{c}^{ji}\\\end{array}}}$
(9)

The velocity ${\displaystyle {v}^{ij}}$ at the contact point is determined by eq. 10 taking into account the angular and translational velocities of the contacting particles, as shown in fig. 5.

 ${\displaystyle {v}^{ij}=\left({\boldsymbol {\omega }}_{j}\times {\boldsymbol {r}}_{c}^{ji}+{v}_{j}\right)-\left({\boldsymbol {\omega }}_{i}\times {\boldsymbol {r}}_{c}^{ij}+{v}_{i}\right)}$
(10)

In case of contact with a boundary ${\displaystyle b}$, the velocity of the rigid (or deformable) structure at the contact point has to be determined. If finite elements are used to discretize the boundaries, typically the velocities can be interpolated from the nodal velocities by means of the shape functions ${\displaystyle N_{k}}$ (see chapter 4). Equation 10 is then modified to:

 ${\displaystyle {v}^{ib}=\sum _{k=0}^{n_{b}}N_{k}({\boldsymbol {Pc}}^{ib})\cdot {v}^{k}-\left({\boldsymbol {\omega }}_{i}\times {\boldsymbol {r}}_{c}^{ij}+{v}_{i}\right)}$
(11)

The velocity at the contact point can be decomposed in the local reference frame defined at the contact point as:

 ${\displaystyle {v}_{n}^{ij}=\left({v}^{ij}\cdot {n}^{ij}\right)\cdot {n}^{ij}}$ (12.a) ${\displaystyle {v}_{t}^{ij}={v}^{ij}-{v}_{n}^{ij}}$ (12.b)

And thus, the definition of the tangential unit vector becomes:

 ${\displaystyle {t}^{ij}={\frac {{v}_{t}^{ij}}{\lVert {v}_{t}^{ij}\rVert }}}$
(13)

Now the contact force }{F}^{ij}\mbox{ between the two interacting spheres ${\displaystyle i}$ and ${\displaystyle j}$ can be decomposed into its normal ${\displaystyle {F}_{n}^{ij}}$ and tangential ${\displaystyle {F}_{t}^{ij}}$ components (Fig. 5):

 ${\displaystyle {F}^{ij}={F}_{n}^{ij}+{F}_{t}^{ij}=F_{n}{n}^{ij}+F_{t}{t}^{ij}}$
(14)

The forces ${\displaystyle F_{n}}$, ${\displaystyle F_{t}}$ are obtained using a contact constitutive model. Standard models in the DEM are characterized by the normal ${\displaystyle k_{n}}$ and tangential ${\displaystyle k_{t}}$ stiffness, normal ${\displaystyle d_{n}}$ and tangential ${\displaystyle d_{t}}$ local damping coefficients at the contact interface and Coulomb friction coefficient ${\displaystyle \mu }$ represented schematically in Fig. 6 for the case of two discrete spherical particles.

 Figure 6: DEM standard contact rheology

Some of the most common models are detailed in the next section 2.5. The models used in a combined DE-FE strategy are described in Chapter 4.

### 2.5 Contact models

The contact between two particles poses in general a complex problem which is highly non-linear and dependent on the shape, material properties, relative movement of the particles, etc. Theoretically, it is possible to calculate these forces directly from the deformation that the particles experience during the contact [52]. In the framework of the DEM however, simplified models are used which depend on a few contact parameters such as the particles relative velocity, indentation, radius and material properties such as the Young's modulus and Poisson's ratio toghether with some parameters that summarize the local loss of energy during the contacts.

The most common model is the so-called linear spring-dashpot model (LS+D) proposed by Cundall and Strack [2] which has an elastic stiffness device and a dashpot which introduces viscous (velocity-dependent) dissipation. This model, while being the simplest one, happens to yield nice results as described in the work from Di Renzo and Di Maio [53] for the case of elastic collisions and in the work of Thornton [54] for the case of inelastic collisions. This model is described in section 2.5.1.

In a second level of complexity, we find models that derive from the theory of Hertz-Mindlin and Deresiewicz. Hertz [55] proposes that the relationship between the normal force and normal displacement is non-linear. Mindlin and Deresiewicz [56] proposed a general tangential force model where the force-displacement relationship depends on the whole loading history as well as on the instantaneous rate of change of the normal and tangential force or displacement. This model was adapted to the DEM by Vu-Quoc and Zhang [57] and later by Di Renzo and Di Maio [53]. This model is quite complicated and requires high computational effort. Other simplified models exist [58,53,54] which consider only the non-slip regime of the Mindlin theory [59]. The model presented in section 2.5.2 is the simplified model by Thornton et al. [54], labeled HM+D.

Other models exist in literature which introduce plastic energy dissipation in a non-viscous manner. This includes the semi-latched spring force-displacement models of Walton and Braun [60] which uses, for the normal direction, different spring stiffnesses for loading and unloading. Similarly, Thornton [61] introduced a model in which the evolution of the contact pressure can be approximated by an elastic stage up to some limit followed by a plastic stage.

Unless the contrary is specified, the HM+D contact law will be used in examples of the thesis. In general, the criterion suggested here is to employ this model with the real material parameters whenever the physics of the contact have influence in the simulation results. In other cases, where the details of the contacts are not relevant, both linear and Hertzian contact laws can be used as a mere penalty technique being the stiffness value a trade-off between simulation time and admissible interpenetration. The model presented for the cohesive materials in chapter 5 is an extension of the linear law (LS+D).

#### 2.5.1 Linear contact law (LS+D)

The model presented here corresponds to a modification of the original model from Cundall and Strack [2] in which the damping force is included in the way the contact rheology has been presented (figure 6).

#### Normal force

In the basic linear contact law the normal contact force ${\displaystyle F_{n}}$ is decomposed into the elastic part ${\displaystyle F_{ne}}$ and the damping contact force ${\displaystyle F_{nd}}$:

 ${\displaystyle F_{n}=F_{ne}+F_{nd}}$
(15)

The damping part is a viscous force which models the loss of energy during a contact. It also serves as a numerical artifact that helps to decrease oscillations of the contact forces which is useful when using an explicit time scheme.

#### Normal elastic force

The elastic part of the normal compressive contact force ${\displaystyle F_{ne}}$ is, in the basic model, proportional to the normal stiffness ${\displaystyle k_{n}}$ and to the indentation (or interpenetration) ${\displaystyle \delta _{n}}$ (eq. 7) of the two particle surfaces, i.e.:

 ${\displaystyle F_{ne}=k_{n}\delta _{n}}$
(16)

Since no cohesive forces are accounted in the basic model. eq. 16 holds only if ${\displaystyle \delta _{n}>0,}$, otherwise ${\displaystyle F_{ne}=0}$. The cohesive contact will be considered in Chapter 5.

#### Normal contact damping

The contact damping force is assumed to be of viscous type and given by:

 ${\displaystyle F_{nd}=c_{n}\cdot {\dot {\delta }}_{n}}$
(17)

where ${\displaystyle {\dot {\delta }}_{n}}$ is the normal relative velocity of the centres of the two particles in contact, defined by:

 ${\displaystyle {\dot {\delta }}_{n}=-({\dot {C}}_{j}-{\dot {C}}_{i})\cdot {n}^{ij}}$
(18)

The damping coefficient ${\displaystyle c_{n}}$ is taken as a fraction ${\displaystyle \xi }$ of the critical damping ${\displaystyle c_{c}}$ for the system of two rigid bodies with masses ${\displaystyle m_{i}}$ and ${\displaystyle m_{j}}$ connected by a spring of stiffness ${\displaystyle k_{n}}$ with:

 ${\displaystyle c_{n}=\xi c_{c}=2\xi {\sqrt {m_{eq}k_{n}}}}$
(19)

with ${\displaystyle 0<{\xi }\leq 1}$ and ${\displaystyle m_{eq}}$ is the equivalent mass of the contact,

 ${\displaystyle m_{eq}={\frac {m_{i}m_{j}}{m_{i}+m_{j}}}}$
(20)

The fraction ${\displaystyle \xi }$ is related with the coefficient of restitution ${\displaystyle e_{n}=-{\dot {\delta }}_{n}^{\textrm {after}}/{\dot {\delta }}_{n}^{\textrm {before}}}$, which is a fractional value representing the ratio of speeds after and before an impact, through the following expression (see [62]):

 ${\displaystyle \xi ={\frac {-\ln {e_{n}}}{\sqrt {\pi ^{2}+{\ln ^{2}{e_{n}}}}}}}$
(21)

#### Contact duration

The equation of motion describing the collision of particles with the LS+D model in the normal direction is achieved solving the differential equation resultant from the application of equation 2 in a frame centred at the point of contact:

 ${\displaystyle m_{eq}\,{\ddot {\delta }}_{n}=-(k_{n}\delta _{n}+c_{n}{\dot {\delta }}_{n})}$
(22)

Eq. 22 can be rewritten as [62]:

 ${\displaystyle {\ddot {\delta }}_{n}+2\Psi ({\dot {\delta }}_{n})+\Omega _{0}^{2}\,\delta _{n}=0}$
(23)

Where ${\displaystyle \Omega _{0}={\sqrt {k_{n}/m_{eq}}}}$ is the frequency of the undamped harmonic oscillator and ${\displaystyle \Psi =\xi \Omega _{0}=c_{n}/(2m_{eq})}$ is the part accounting for the energy dissipation.

The solution of the eq. 22 for the initial conditions ${\displaystyle \delta _{n}=0}$ and ${\displaystyle {\dot {\delta }}_{n}=v_{0}}$ and for the sub-critical damped case1 (${\displaystyle \Omega _{0}^{2}-\Psi ^{2}>0}$ or ${\displaystyle \xi {<1}}$) reads:

(1) The cases of critical and super-critical damping yield to other solutions which can be found in [63]

 ${\displaystyle \delta _{n}(t)=\left(v_{0}/\Omega \right)\,e^{-\Psi t}\,\sin {(\Omega t)}\qquad {\textrm {with}}\qquad \Omega ={\sqrt {\Omega _{0}^{2}-\Psi ^{2}}}}$
(24)

And the relative normal velocity of the spheres:

 ${\displaystyle {\dot {\delta }}_{n}(t)=\left(v_{0}/\Omega \right)\,e^{-\Psi t}\,\left(-\Psi \sin {(\Omega t)}+\Omega \cos {(\Omega t)}\right)}$
(25)

Now the contact duration can be determined from the condition ${\displaystyle \delta _{n}(t_{c})=0}$, which combined with eq. 24 gives:

 ${\displaystyle t_{c}=\pi /\Omega }$
(26)

Note that the contact duration does not depend on the initial approaching velocity ${\displaystyle {\dot {\delta }}_{n}(t)}$ which is obviously wrong as the formulation is not derived from the theory of elasticity [52] (see section 4.4 for more details).

The coefficient of restitution can be rewritten as:

 ${\displaystyle e_{n}={\frac {-{\dot {\delta }}_{n}(t_{c})}{{\dot {\delta }}_{n}(0)}}=e^{-\pi \Psi /\Omega }}$
(27)

The inverse relationship allows the determination of the parameter ${\displaystyle c_{n}}$ of the model from the restitution coefficient ${\displaystyle e_{n}}$, with the intermediate calculation of ${\displaystyle \Psi }$:

 ${\displaystyle \Psi ={\frac {-\ln {e_{n}}}{\sqrt {\pi ^{2}+{\ln ^{2}{e_{n}}}}}}\,\Omega _{0}}$
(28)

Finally, the maximum indentation can be obtained from the condition ${\displaystyle {\dot {\delta }}_{n}(t)=0}$:

 ${\displaystyle \delta _{max}=(v_{0}/\Omega _{0})e^{-{\frac {\Psi }{\Omega }}\arctan {(\Omega /\Psi )}}}$
(29)

#### Note on tensional forces

It has been appointed by different authors [62,64,54] that this simple model presents unrealistic tension force when the particles are separating if the damping force is large enough (Fig. 7). Normally in the implementation of the codes the normal force is constrained to be exclusively positive, i.e ${\displaystyle F_{n}\geq 0}$ always, as no tractions occur in frictional cohesion-less contacts. In this situation the definition of the contact duration should be modified as it has been derived by Schwager and Pöschel [64].
 ] Figure 7: The different stages of a normal collision of spheres with a viscous damped model. Taken from: Fig. 1 in Schwager and Pöschel [64]

The determination of the damping coefficients and the maximum indentation vary accordingly. It is not possible to derive an explicit expression for the damping coefficient ${\displaystyle c_{n}}$ in function of the restitution coefficient ${\displaystyle e_{n}}$. Fitting curves are proposed in [54].

#### Tangential frictional contact

In the original model from Cundall and Strack [2] the relationship between the elastic shear force ${\displaystyle F_{t}}$ and the relative tangential displacement ${\displaystyle \Delta s}$ is defined through a regularized Coulomb model. The update of the tangential force at time step ${\displaystyle n+1}$ reads:

 ${\displaystyle F_{t}^{n+1}=\min \left(\mu F_{n},\;F_{t}^{n}+k_{t}\Delta s^{n+1}\right)}$
(30)

Several authors (including the original paper) calculate the increment of tangential displacement at a given time step ${\displaystyle n}$, ${\displaystyle \Delta s^{n}}$, as ${\displaystyle \lVert {v}_{t}^{ij,n}\rVert \cdot \Delta t}$. In our in-house code implementation it is calculated from the integration of the relative displacement and rotation in the local frame:

 ${\displaystyle {\Delta s}^{n}=\lVert {u}^{ij}\cdot {t}^{ij}\rVert }$ (31.a) ${\displaystyle {u}^{ij}=\left({\boldsymbol {\Theta }}_{j}\times {\boldsymbol {r}}_{c}^{ji}+{u}_{j}\right)-\left({\boldsymbol {\Theta }}_{i}\times {\boldsymbol {r}}_{c}^{ij}+{u}_{i}\right)}$ (31.b)

In the original paper [2] the damping is included only during the non-sliding phase (${\displaystyle F_{t}\leq \mu F_{n}}$) and it is applied afterwards as an extra force which opposes the relative velocity. The magnitude of the damping force is evaluated as ${\displaystyle c_{t}\cdot \lVert {v}_{t}^{ij}\rVert }$ where ${\displaystyle c_{t}}$ the tangential damping coefficient. In other authors' works and also in our code implementation it is chosen to include the dissipation in the check for sliding. In case of sliding (${\displaystyle F_{t}=\mu Fn}$), extra decision on how to distribute the resultant tangential force in elastic and dissipative part have to be taken. This will not be discussed here. Eq. 30 modifies as:

 ${\displaystyle {F}_{t}^{\textrm {trial}}={F}_{t}^{n}+k_{t}\Delta s^{n+1}\,{t}^{ij}+c_{t}{v}_{t}^{ij,n+1}}$ (32.a) ${\displaystyle {F}_{t}^{n+1}=\min \left(\mu F_{n},\;\lVert {F}_{t}^{\textrm {trial}}\rVert \right){\frac {{F}_{t}^{\textrm {trial}}}{\lVert {F}_{t}^{\textrm {trial}}\rVert }}}$ (32.b)

The previous time step forces are transferred from its previous local coordinate frame to the new local contact frame with a rotation of the force vector (section 2.7.1).

#### Selection of the stiffness and damping parameters

The selection of the normal stiffness ${\displaystyle k_{n}}$ is, in the LS+D model, a design parameter. The general rule of thumb is that the value of ${\displaystyle k_{n}}$ should be large enough to avoid excessive particle inter-penetration but at the same time should be small enough to permit reasonable simulation time steps (section 2.6.4) [65].

Cundall and Strack [2] investigated several values for the relation ${\displaystyle \kappa =kt/kn}$ in the range ${\displaystyle [2/3,1]}$, obtained from the following expression:

 ${\displaystyle \kappa ={\frac {2(1-\nu )}{2-\nu }}}$
(33)

The values for the damping in the original paper [2] are selected as a proportion ${\displaystyle \beta }$ of the respective stiffnesses:

 ${\displaystyle c_{n}=\beta k_{n}}$ (34.a) ${\displaystyle c_{t}=\beta k_{t}}$ (34.b)

Normally, the selection of ${\displaystyle \beta }$ will be based on the desired restitution coefficient through eq. 19 and 21. Alternatively, Schäfer [66] suggests a value of ${\displaystyle k_{t}}$ equal to two-sevenths of the normal stiffness coefficient and a damping ${\displaystyle c_{t}}$ as half of the normal damping coefficient. Thornton [67], in his turn, suggests a value of ${\displaystyle k_{n}}$ that yields the same contact duration as the one predicted by the Hertzian theory (section 2.5.2).

(1) The cases of critical and super-critical damping yield to other solutions which can be found in [63]

#### 2.5.2 Hertzian contact law (HM+D)

As introduced in section 2.5, there exist in literature several contact laws under the framework the Hertzian contact theory [55]. The model chosen for this dissertation is an adaptation of the one referred as HM+D model in the work by Thornton [54] due to its balance between simplicity and accuracy in both elastic [53,67] and inelastic collisions [54]. This is a model based on the original one by Tsuji [58] in which the tangential spring is provided by the no slip theory of Mindlin [59].

The magnitude of the normal force can be calculated as:

 ${\displaystyle F_{n}={\frac {2}{3}}k_{n}\delta _{n}+c_{n}{\dot {\delta }}_{n}}$
(35)

The tangential update has two branches whether the normal force is increasing (loading phase) or decreasing (unloading case). For the loading phase the tangential force is increased as usual due to the tangential displacement (Eq. 36.a). In the unloading phase, however (Eq. 36.b), the tangential force must be reduced (even with no tangential displacement) due to the reduction in the contact area. The interpretation of this is that the previous tangential force can not longer be supported [54].

 ${\displaystyle F_{te}^{n+1}=F_{te}^{n}+k_{t}^{n+1}\Delta s^{n+1}{\textrm {for}}\qquad \Delta F_{n}\geq 0}$ (36.a) ${\displaystyle F_{te}^{n+1}=F_{te}^{n}\left({\frac {k_{t}^{n+1}}{k_{t}^{n}}}+k_{t}^{n+1}\right)\Delta s^{n+1}{\textrm {for}}\qquad \Delta F_{n}<0}$ (36.b)

Finally, the check for sliding is performed restricting the maximum tangential force to the Coulomb's friction limit:

 ${\displaystyle F_{t}^{\textrm {trial}}=F_{te}^{n+1}+c_{t}{v}_{t}^{ij}}$ (37.a) ${\displaystyle F_{t}^{n+1}=F_{t}^{\textrm {trial}}\qquad {\textrm {if}}\qquad F_{t}^{n+1}<\mu F_{n}}$ (37.b) ${\displaystyle F_{t}^{n+1}=\mu F_{n}\qquad {\textrm {if}}\qquad F_{t}^{n+1}\geq \mu F_{n}}$ (37.c)

The stiffness parameters were described by Tsuji [58] following from the Hertz theory [55] and the works of Mindlin and Deresiewicz [56]:

 ${\displaystyle k_{n}=2E^{*}{\sqrt {R_{eq}\delta _{n}}}}$ (38.a) ${\displaystyle k_{t}=8G^{*}{\sqrt {R_{eq}\delta _{n}}}}$ (38.b)

The same for the damping parameters:

 ${\displaystyle c_{n}=2\xi {\sqrt {m_{eq}k_{n}}}}$ (39.a) ${\displaystyle c_{t}=2\xi {\sqrt {m_{eq}k_{t}}}}$ (39.b)

The expressions presented here (eq. 9.2 and 9.2) are a generalization to the case of two spheres ${\displaystyle i}$ and ${\displaystyle j}$ colliding with different values of ${\displaystyle R}$, ${\displaystyle E}$, ${\displaystyle \nu }$ and ${\displaystyle m}$. This generalization includes the case of a sphere ${\displaystyle i}$ colliding with a fixed wall ${\displaystyle j}$ which will be discussed in section 2.5.3.

 ${\displaystyle R_{eq}=R_{i}R_{j}/(R_{i}+R_{j})}$ (40.a) ${\displaystyle m_{eq}=m_{i}m_{j}/(m_{i}+m_{j})}$ (40.b) ${\displaystyle E_{eq}^{*}=\left((1-\nu _{i}^{2})/E_{i}+(1-\nu _{j}^{2})/E_{j}\right)^{-1}}$ (40.c) ${\displaystyle G_{eq}^{*}=\left((2-\nu _{i})/G_{i}+(2-\nu _{j})/G_{j}\right)^{-1}}$ (40.d)

Although the selection of the stiffness has here a physical meaning, it is common practice however, to diminish its value to increase the calculation speed in simulations where the correct contact duration and rebound angles are not of capital importance. The derivation of the force-displacement relationship and the collision time by the Hertzian theory are described in the Appendix A.

#### 2.5.3 Contact with rigid boundaries

Rigid boundaries are commonly introduced in a DE simulation to model the interaction of particles with mechanical components that can be either fixed or have an imposed rigid body motion. Although they are normally discretized with a FE mesh for contact detection purposes (section 3.1), they are not calculated by means of a FE procedure. The rheology of a particle ${\displaystyle i}$ contacting a FE ${\displaystyle j}$ is presented in figure 8.
 Figure 8: DE-FE standard contact rheology

Same as for DE/DE contact, Hertzian contact law is preferred to model the contacts or impacts in a physical basis. Alternatively the linear contact law can still be used as basic penalty method. The adaptation of the presented Hertzian contact law to the case of rigid boundaries is straightforward, it simply requires the particularization of the equivalent contact parameters summarized in 9.2 setting: ${\displaystyle R_{j}\rightarrow \infty }$ and ${\displaystyle m_{j}\rightarrow \infty }$. The normal stiffness of the wall is left as an input parameter so that a certain elasticity of the wall can be modelled. Since the tangential displacement of the wall will be in most cases much smaller than the particle's one, it is recommended to be set ${\displaystyle G_{j}\rightarrow \infty }$ [58]. The equivalent values become:

 ${\displaystyle R_{eq}=R_{i}}$ (41.a) ${\displaystyle m_{eq}=m_{i}}$ (41.b) ${\displaystyle E_{eq}^{*}=\left((1-\nu _{i}^{2})/E_{i}+(1-\nu _{j}^{2})/E_{j}\right)^{-1}}$ (41.c) ${\displaystyle G_{eq}^{*}=G_{i}/(2-\nu _{i})}$ (41.d)

The stiffness and damping parameters are modified accordingly inserting these equivalent values in eq. 9.2 and eq. 9.2. The friction value to be employed in this case is a new parameter to be introduced, which is characteristic of the contact between the two materials involved and might be different from the particle-particle friction.

Additionally, special contact laws can be applied which model other effects such as a specific dynamic response, wear, plasticity, thermal coupling, etc. [68,69,70].

#### 2.5.4 Rolling friction

It should be noted that the use of spherical particles to represent real materials may lead to excessive rotation. To avoid this effect the rolling resistance approach has been used. This approach consists in imposing a virtual resistive torque which is proportional to the normal contact force and opposites the rolling direction. The rolling resistance torque ${\displaystyle {T}^{r}}$ is defined as;

 ${\displaystyle {T}^{r}=-\eta _{r}R_{r}|{F}^{n}|{\frac {{\boldsymbol {\omega }}^{rel}}{|{\boldsymbol {\omega }}^{rel}|}}}$
(42)

where ${\displaystyle \eta _{r}}$ is the rolling resistance coefficient that depends on the material, ${\displaystyle R_{r}}$ is the smallest radius of the DEs in contact and ${\displaystyle {\boldsymbol {\omega }}^{rel}}$ the relative angular velocity of both DEs. Note that ${\displaystyle R_{r}=R_{i}}$ for the case where particle ${\displaystyle i}$ is in contact with a wall (${\displaystyle R_{j}\rightarrow \infty }$).

An improvement to the classical Rolling Resistance Model A presented by Wensrich and Katterfeld [71] has been developed by Irazábal [72] in order to avoid the instabilities that appear when ${\displaystyle {\boldsymbol {\omega }}^{rel}}$ is close to ${\displaystyle 0}$.

### 2.6 Time integration

The equations of motion introduced in section 2.3 can be numerically integrated to obtain a solution of the problem. Traditionally there are two strategies to achieve this: a) An explicit scheme where the information at the current (or previous steps) suffices to predict the solution at the next step. b) An implicit scheme, which requires the solution of a non-linear system of equations to compute the state at the new time step. The disadvantage of the explicit schemes is that they require the time step to be below a certain limit in order to be stable. Implicit schemes instead, are unconditionally stable and thus, allow for larger time steps.

Some analysis on both implicit and explicit methods for discrete element simulations showed that the second one is generally preferable [73,74]. Implicit algorithms turn to be not efficient for DEM simulations because of the nature of the dynamics of particles where relatively large motions are simulated combined with very small characteristic relative displacements between particles during contact events. In order to correctly capture the dynamics of the contact, the time resolution should be several times smaller than the duration of these contacts [74]. Under this condition, the explicit integration yields sufficient accuracy and the time step is generally below its stability limits (see section 2.6.4). Following the same reasoning, low order explicit schemes are usually preferred rather than higher order ones. Another important outcome of the use of an explicit integration is the easier parallelization of the code and the avoidance of linearization and employment of system solvers.

In other situations where the same contacts are kept for large simulation times, such as cohesive models for DE (chapter 5), the use implicit schemes can be advantageous. Otherwise, the stiffness matrices have to be rebuild, in general, at each iteration and time step due to the formation and destruction of contacts. Amongst the most popular implicit approaches in DEM is the Discontinuous Deformation Analysis [75].

#### 2.6.1 Explicit integration schemes

In the present dissertation an explicit integration is used. Next, four different one-step integration algorithms with similar computational cost are described and compared in this section. The derivation of these methods comes from the application of the Taylor series approximation to the second order differential equations of motion (2) that describes the problem.

 ${\displaystyle f(t+\Delta t)=f(t)+{\frac {f'(t)}{1!}}\Delta t+{\frac {f''(t)}{2!}}{\Delta t}^{2}+{\frac {f'''(t)}{3!}}{\Delta t}^{3}+...}$
(43)

#### Forward Euler

The forward difference approximation of the first derivative of a function reads as:

 ${\displaystyle f'(t)={\frac {1}{\Delta t}}(f(t+\Delta t)-f(t))}$
(44)

The terms can be rearranged to obtain an integration formula:

 ${\displaystyle f(t+\Delta t)=f(t)+\Delta tf'(t)}$
(45)

which is applied to the integrate the acceleration and the velocity respectively:

 ${\displaystyle {\dot {u}}^{n+1}={\dot {u}}^{n}+\Delta t\,{\ddot {u}}^{n}}$ (46) ${\displaystyle {u}^{n+1}={u}^{n}+\Delta t\,{\dot {u}}^{n}}$ (47)

The truncation error of the Taylor expansion approximations are of ${\displaystyle {\mathcal {O}}(N^{2})}$. Hence, the method is referred to as a first order approximation of the displacement and velocities.

#### Symplectic Euler

The Symplectic Euler is a modification of the previous method which uses a backward difference approximation for the derivative of the position:

 ${\displaystyle f'(t)={\frac {1}{\Delta t}}(f(t)-f(t-\Delta t))}$
(48)

The algorithm is as follows:

 ${\displaystyle {\dot {u}}^{n+1}={\dot {u}}^{n}+\Delta t\,{\ddot {u}}^{n}}$ (49) ${\displaystyle {u}^{n+1}={u}^{n}+\Delta t\,{\dot {u}}^{n+1}}$ (50)

This way a higher accuracy and order of convergence can be achieved as it is shown in the numerical convergence analysis performed in the following section 2.6.3.

#### Taylor Scheme

The Taylor schemes are a family of integration methods which make use of the Taylor expansion (43) to approximate the next values of the variable of interest. If the series are truncated at the first derivative for the velocity and at the second derivative for the position, the following integration rule is obtained:

 ${\displaystyle {\dot {u}}^{n+1}={\dot {u}}^{n}+\Delta t\,{\ddot {u}}^{n}}$ (51) ${\displaystyle {u}^{n+1}={u}^{n}+\Delta t\,{\dot {u}}^{n}+{\frac {1}{2}}{\Delta t}^{2}\,{\ddot {u}}^{n}}$ (52)

Which is a first order integrator for the velocity and a second order integrator for the position.

#### Velocity Verlet

This algorithm is sometimes simply called Central Differences [76,38] and some other times it is interpreted as the velocity form of the Verlet algorithm [74,77]. It also coincides with the special case of the Newmark-beta method [78] with ${\displaystyle \beta =0}$ and ${\displaystyle \gamma =1/2}$. The derivation presented here is the same as it is described by Belytschko in [76]. The central difference formula is written as:

 ${\displaystyle f'(t)={\frac {1}{\Delta t}}(f(t+1/2\Delta t)-f(t-1/2\Delta t))}$
(53)

Applying it to the velocity at an intermediate position ${\displaystyle n+1/2}$:

 ${\displaystyle {\dot {u}}^{n+1/2}={\frac {1}{\Delta t}}({u}^{n+1}-{u}^{n})}$
(54)

and to the acceleration at the time step ${\displaystyle n}$:

 ${\displaystyle {\ddot {u}}^{n}={\frac {1}{\Delta t}}({\dot {u}}^{n+1/2}-{\dot {u}}^{n-1/2})}$
(55)

Inserting equation 54 and its counterpart for the previous time step (${\displaystyle {v}^{n-1/2}}$) into equation 55, the central difference formula for the second derivative of the displacement is obtained:

 ${\displaystyle {\ddot {u}}^{n}={\frac {1}{{\Delta t}^{2}}}({u}^{n+1}-2{u}^{n}+{u}^{n-1})}$
(56)

The algorithm follows from the rearrangement of equations 54 and 55

 ${\displaystyle {\dot {u}}^{n+1/2}={\dot {u}}^{n-1/2}+\Delta t\,{\ddot {u}}^{n}}$ (57) ${\displaystyle {u}^{n+1}={u}^{n}+\Delta t\,{\dot {u}}^{n+1/2}}$ (58)

Since it may be necessary to have both velocity and position evaluated at every time step of the discretization, a split in the calculation of ${\displaystyle {\dot {u}}^{n+1/2}}$ can be performed.

 ${\displaystyle {\dot {u}}^{n}={\dot {u}}^{n-1/2}+1/2\,\Delta t\,{\ddot {u}}^{n}}$ (59) ${\displaystyle {\dot {u}}^{n+1/2}={\dot {u}}^{n}+1/2\,\Delta t\,{\ddot {u}}^{n}}$ (60)

The implementation of the method is summarized in the following table:

 Initialization of the scheme. ${\displaystyle n=0}$, ${\displaystyle {\ddot {u}}^{0}={F}^{0}/m}$ while ${\displaystyle t Update step: ${\displaystyle n=n+1}$, ${\displaystyle t=t+\Delta t}$ First velocity update: ${\displaystyle {\dot {u}}^{n+1/2}={\dot {u}}^{n}+1/2\;\Delta t\;{\ddot {u}}^{n}}$ Position update: ${\displaystyle {u}^{n+1}={u}^{n}+\Delta t\;{\dot {u}}^{n+1/2}}$ Calculate forces ${\displaystyle {F}^{n+1}={F}\left({u}^{n+1},{\dot {u}}^{n+1/2}\right)}$ Calculate acceleration: ${\displaystyle {\ddot {u}}^{n+1}={F}^{n+1}/m}$ Second velocity update: ${\displaystyle {\dot {u}}^{n+1}={\dot {u}}^{n+1/2}+1/2\;\Delta t\;{\ddot {u}}^{n+1}}$

This is the selected scheme for the examples in this dissertation.

#### 2.6.2 Integration of the rotation

The particular case of spherical particles simplifies the equations for the rotation of rigid bodies yielding to equation 3. Some authors [79,80,60] adapt a simple central difference scheme to integrate the equations:

 ${\displaystyle {\dot {\omega }}_{i}^{n}={\frac {{T}_{i}^{n}}{I_{i}}}\,,}$ (61) ${\displaystyle {\omega }_{i}^{n+1/2}={\omega }_{i}^{n-1/2}+{\dot {\omega }}_{i}^{n}{\Delta t}}$ (62)

The vector of incremental rotation ${\displaystyle {\Delta \theta }^{n+1}}$ is then calculated as:

 ${\displaystyle {\Delta \theta }_{i}^{n+1}={\omega }_{i}^{n+1/2}{\Delta t}}$
(63)

Knowledge of the incremental rotation suffices to update the tangential contact forces. If necessary, it is also possible to track the rotational position of particles, as detailed in section 2.7.1.

#### 2.6.3 Accuracy analysis

In this section the error of the different integration methods previously introduced is addressed by means of accuracy and convergence analysis. Three cases representative of translational motion occurring in a DEM simulation are analysed here: free parabolic motion, normal contact between two spheres using a linear contact law and normal contact between two spheres using a Hertzian contact law. The description of the test examples is in figure 9. A similar analysis has been performed by Samiei [74] for the comparison of some explicit schemes against implicit integration.

The case of rotational motion is analysed in section 2.7.3 where a higher order scheme is implemented for the case of a generic rigid body which can be also applied to the spheres. It is shown that the integration of the rotation equation requires higher order schemes for similar levels of accuracy as the one-step methods.

 (a) Set-up parabolic motion (b) Set-up normal contact Figure 9: Examples for the accuracy and convergence analysis on time integration schemes

#### Parabolic motion analysis

An initial upwards velocity of ${\displaystyle 1.0\,m/s}$ is set to a particle situated at the origin of coordinates which moves freely only under the effect of gravity which is set to ${\displaystyle -10\,m/s^{2}}$ during ${\displaystyle 0.2}$ seconds.

A numerical integration of the problem is performed with the presented methods and compared against the analytical solution. The time step is chosen to be a tenth of the total time so that the error of the methods can be easily observed.

 Figure 10: Vertical displacement of a sphere under gravity using ${\displaystyle 10}$ time steps
 Figure 11: Velocity of a sphere under gravity using ${\displaystyle 10}$ time steps
As expected, the velocity is perfectly integrated for any of the analysed schemes since the acceleration is constant over time (figure 11). The position (figure 10) is integrated perfectly by the Taylor Scheme and Velocity Verlet which are second order schemes in displacement.
 Figure 12: Convergence in velocity and displacement for different integration schemes

Figure 12 shows that the Forward Euler and Symplectic Euler schemes have a linear convergence when integrating the position. The convergence is omitted for other schemes and for the velocity since the algorithmic error is null.

#### Normal contact analysis with the LS+D model

Two spheres are set in space with tangential contact (no indentation) and without the effect of the gravity. One of the spheres approaches the other one with an initial velocity in the direction of the vector joining the spheres' centres as depicted in figure 9c. The linear contact law introduced in section 2.5.1 is applied.

The expression for the maximum indentation (eq. 29) for the non-damped case ${\displaystyle (\Psi =0)}$ turns into:

 ${\displaystyle \delta _{max}=v_{0}{\sqrt {\frac {m_{eq}}{k_{n}}}}}$
(64)

And the contact duration (eq. 26):

 ${\displaystyle t_{c}=\pi {\sqrt {\frac {m_{eq}}{k_{n}}}}}$
(65)

The simulation is carried out for the different schemes with a time step corresponding to a contact resolution1 (${\displaystyle CR}$) of 10, i.e. the time step corresponds to a tenth of the contact duration. The parameters of the simulation are summarized in the following Table 2:

 Contact law Linear Contact Law (section 2.5.1) Radius ${\displaystyle 0.01\;m}$ Density ${\displaystyle 100\;kg/m^{3}}$ ${\displaystyle k_{n}}$ ${\displaystyle 520.83\;kN/m}$ Restitution coeff. ${\displaystyle 1.0}$ ${\displaystyle V_{0}}$ ${\displaystyle 0.5\;m/s}$ Contact time ${\displaystyle 4.17\cdot {10}^{-3}\,s}$ ${\displaystyle CR}$ ${\displaystyle 10}$
 Figure 13: Indentation during the collision of two spheres using LS+D with ${\displaystyle CR=10}$
 Figure 14: Velocity during the collision of two spheres using LS+D with ${\displaystyle CR=10}$
Both Symplectic Euler and Velocity Verlet accurately approximate the indentation (Fig. 13). Regarding the velocity, the Verlet scheme is the one with superior accuracy over the other schemes (Fig. 14).
 Figure 15: Convergence in velocity and displacement for the FE and SE schemes

The numerical results for the maximum indentation as well as the exit velocity of the contact have been taken as the measures to evaluate the error for different time steps. Both F.E. and Taylor schemes showed linear convergence in displacement and velocity (Fig. 15). On the other hand, S.E. and V.V. showed quadratic convergence for the displacement and velocities.

#### Normal contact analysis with the HM+D model

Finally, the same test is carried out using a Hertzian contact law (section 2.5.2). The derivation of the contact time duration and other properties of the Hertzian contact are detailed in Appendix A. The simulation parameters are summarized in Table 3. The different schemes are tested with a ${\displaystyle CR=10}$ and the results for the indentation evolution and its time derivative are plotted in Fig. 16 and Fig. 17 respectively.

 Contact law Hertzian Contact Law (section 2.5.2) Radius ${\displaystyle 0.01\;m}$ Density ${\displaystyle 100\;kg/m^{3}}$ Young's modulus ${\displaystyle 1\cdot {10}^{5}\;kN/m^{2}}$ Poisson's ratio ${\displaystyle 0.2}$ Restitution coeff. ${\displaystyle 1.0}$ ${\displaystyle V_{0}}$ ${\displaystyle 0.5\;m/s}$ Contact time ${\displaystyle 1.99\cdot {10}^{-3}\,s}$ ${\displaystyle CR}$ ${\displaystyle 10}$
 Figure 16: Indentation during the collision of two spheres using HM+D with ${\displaystyle CR=10}$
The same conclusions of the linear case can be drawn for the Hertzian contact: the Symplectic Euler and Velocity Verlet accurately approximate the indentation (Fig. 16) while the other schemes present some error. Regarding the velocity, the better scheme is clearly the Verlet scheme (Fig. 17).
 Figure 17: Velocity during the collision of two spheres using LS+D with ${\displaystyle CR=10}$

In terms of convergence, the velocity presented even a higher order than quadratic for the Verlet scheme. It shall be noticed however, that the error of this variable for the selected time steps is too small to draw conclusions on the scheme convergence.

 Figure 18: Convergence in velocity and displacement for different integration schemes

(1) The concept of contact resolution defined as ${\displaystyle CR=t_{c}/\Delta t}$ is discussed in section 2.6.4.

#### 2.6.4 Stability analysis

There are many factors that can cause instabilities in a Discrete Element simulation. The first basic requisite for the time step, in a DEM simulation, is to be stable in terms of the integration scheme. Another significant source of instabilities is the lack of accuracy in the determination of the formation of contacts. In this sense, quantities such as the velocity of the particles and the search frequency play a great role in the overall stability and are not sufficiently studied in the literature. While most of the authors merely perform a scheme stability analysis [81] for the determination of the time step, a large safety factor is applied which reduces the estimated value. This reinforces the idea of using a time step based on the concept of contact resolution [12,82] defined as the number of steps used to resolve a contact event, ${\displaystyle CR=t_{c}/\Delta t}$.

#### Stability of the integration scheme

Explicit integration schemes present a limitation in the time step in order to be numerically stable ${\displaystyle \Delta t\leq \Delta t_{cr}}$. Belytschko [76] shows that the critical time step ${\displaystyle {\Delta t}_{cr}}$ for a central difference method is determined by the highest natural frequency of the system }\omega _{max}\mbox{ as:

 ${\displaystyle {\Delta t}_{cr}={\frac {2}{\omega _{max}}}}$
(66)

Exact calculation of the highest frequency }\omega _{max}\mbox{ requires the solution of the eigenvalue problem defined for the whole system of connected rigid particles. In an approximate solution procedure, an eigenvalue problem can be defined separately for every rigid particle using the linearized equations of motion. The maximum frequency is estimated as the largest of the natural frequencies of the mass-spring systems defined for all the particles with one translational and one rotational degree of freedom:

 ${\displaystyle \omega _{max}=\operatorname {max} _{i}\omega _{i}}$
(67)

And the natural frequency for each mass-spring system (contact) is defined as:

 ${\displaystyle \omega _{i}={\sqrt {\frac {k}{m_{i}}}}}$
(68)

being ${\displaystyle k}$ the spring stiffness and ${\displaystyle m_{i}}$ the mass of particle ${\displaystyle i}$. Now, for the case with no damping, it is possible to rewrite the critical time step as:

 ${\displaystyle \Delta t_{cr}=\operatorname {min} _{i}\,2{\sqrt {\frac {m_{i}}{k}}}}$
(69)

The effective time step is considered as a fraction of the critical time step:

 ${\displaystyle \Delta t=\beta \Delta t{cr}}$
(70)

The fraction ${\displaystyle \beta \in [0,1]}$ has been studied by different authors. O'Sullivan and Bray in [81] recommend values close to ${\displaystyle \beta =0.17}$ for 3D simulation, and ${\displaystyle \beta =0.3}$ for the 2D case.

If damping exists, the critical time increment is modified with the fraction of the critical damping ${\displaystyle \xi }$ corresponding to the highest frequency }\omega _{max}\mbox{ in the following way [76]:

 ${\displaystyle {\Delta t}_{cr}={\frac {2}{\omega _{max}}}\left({\sqrt {1+\xi ^{2}}}-\xi \right)}$
(71)

Further details are given in section 4.4.1 where the critical time step for a explicit finite element procedure is discussed.

#### Example of the scheme stability

An example is presented here to show the performance of the different schemes for time steps near the critical one and smaller. A sphere of radius ${\displaystyle R=4\,mm}$ and density ${\displaystyle 2.000\,kg/m}$ oscillates between two parallel plates which are separated ${\displaystyle 7\,mm}$ using a linear contact law with stiffness ${\displaystyle k_{n}=1\,N/m}$. The sphere presents an initial indentation with the top plate of ${\displaystyle 1\,mm}$ (fig. 19a). The example tries to simulate the instability effects that can occur locally in a system with dense particle packings.

The linear mass-spring system has a theoretical frequency1 of ${\displaystyle \omega ={\sqrt {2k_{n}/m}}=61.08rad/s}$ which yields to a critical time step ${\displaystyle {\Delta t}_{cr}=0.03275\,s}$. The results for the four schemes are presented (fig. 19) using time steps: ${\displaystyle \Delta t=0.03275\,s}$, ${\displaystyle \Delta t=0.00300\,s}$ and ${\displaystyle \Delta t=0.00010\,s}$.

The results show how the Velocity Verlet is the only scheme which has an acceptable performance in the limit of the critical time step (fig. 19b) as it is a second order scheme. It was found that for a slightly larger time step the Velocity Verlet scheme becomes also unstable as predicted by the criterion in eq. 69. Symplectic Euler, which showed properties similar to a second order scheme in terms of accuracy, does not unstabilize but presents a wrong prediction of the amplitude. As it can be seen in figure 19c the first order schemes are still unstable even for a time step which is ten times smaller than the critical one, being Forward Euler the most unstable one. Finally, in figure 19d it is shown that all methods converge to the analytical solution as the time step diminishes.
 (a) Setup of the example (b) Position evolution for ${\textstyle \Delta t=0.03275}$ (c) Position evolution for ${\textstyle \Delta t=0.00300}$ (d) Position evolution for ${\textstyle \Delta t=0.00010}$ Figure 19: Setup and results for the position of the sphere between the plates

#### Stability due to lack of accuracy

The lack of accuracy can produce instabilities in a DEM simulation. The easiest way to explain it is to imagine a particle travelling with a very large velocity towards another particle or a wall; while the critical time step was shown to be independent of the velocity (eq. 71), a large velocity will imply inaccuracy in the detection of the contact and this translates into an indentation that can be unboundedly large and thus yielding to an unrealistic increase in the energy. This can also be interpreted as an insufficient resolution of the contact.

An example of this effect is found in the work by Ketterhagen et al. [12] where an analysis of how the time step affects the mean stress tensor measurements in two-dimensional granular shear flow simulations is performed. For a time step small enough the simulation results for the stress tensor (or any other variable) should be independent of the time step size. The studies performed using a linear contact model and several stiffness values showed that for a ${\displaystyle CR=15}$ the error in the stress measurement was below ${\displaystyle 2.5\%}$ while higher time steps yielded a sudden increase in the error up to values above ${\displaystyle 10\%}$. These inaccuracies may introduce instabilities as it was shown by their results which are referenced here in Fig. 20.

 ] Figure 20: Stress measurement error in shear flow simulations. Taken from: Fig. 4 in Ketterhagen et al. [12]

A commonly accepted approach as an alternative to the critical time step criterion (section 2.6.4) is to select the time step of the simulation in function of the characteristic duration of the contacts, i.e, by means of the contact resolution. No agreement is found when addressing a correct value for the ${\displaystyle CR}$, several authors recommend values around ${\displaystyle CR=50}$ which is quite conservative (See [83,84,60,65]), Keterhagen et al. recommends a contact resolution of ${\displaystyle CR=33}$ while others like Dury [85] use larger time steps: ${\displaystyle CR=15}$. O'Sullivan [81] determines values of ${\displaystyle CR}$ in the range ${\displaystyle [6-10]}$ using a central differences scheme with regular mono-disperse (same radius) meshes.

Summarizing, there is not a unique solution for the problem of selecting a suitable time step. It depends on many factors such as the mesh, the integration scheme, the type of simulation, the material parameters, the contact law, etc. Our suggestion is to estimate a characteristic contact time of the problem and then select a time step based on the ${\displaystyle CR}$ criterion in the range ${\displaystyle [10-50]}$ depending on the conditions of the problem and the accuracy desired. This will be in general much lower than the critical time step.

(1) The ${\displaystyle 2}$ multiplying the stiffness comes from the fact that this is not a single mass-spring system, instead two plates are contributing to the stiffness of the system.

#### 2.6.5 Computational cost

From the accuracy and stability analysis it is clear that Velocity Verlet and Symplectic Euler are much superior than the Forward Euler and Taylor Scheme, being Velocity Verlet the best one among these four one-step schemes. The final aspect to take into consideration is the computational cost of the method. Simulations in real applications involve millions of particles and can also comprehend millions of time evaluations. The example described in section 6.2.2 is used here to calculate ${\displaystyle 1.000}$ time steps. It includes approximately 30.000 spheres contacting among them and also with around ${\displaystyle 2.500}$ rigid finite elements. The test has been run in a personal computer with an Intel Core i7 processor of 4 Gb RAM and 2.93 GHz.

 Scheme F. Euler Taylor S. Euler V. Verlet Time (s) 169.61 170.04 169.64 174.28

The results showed similar computational times for the four schemes. Velocity Verlet performed ${\displaystyle 2.7\%}$ slowlier which is insignificant considering the advantages found in terms of accuracy and stability in the integration of velocities. Obviously, it will vary in every computer but in general lines we determine that it is worth to employ a Velocity Verlet scheme.

### 2.7 Particle shape

In granular matter, the effects of the particle geometry are crucial in the behaviour of the particles as a bulk or as individuals [86]. Often, the phenomenological approach is considered and the granular media are modelled with spheres as it is the cheapest and most efficient option for simulating a large amount of particles [79]. Alternatively, if we want a method which is purely based on contact and other interaction forces, the real geometry of the particles have to be well represented.

Among the most common methods there is the use superquadrics, which permits a wide range of symmetric convex shapes [4], the Granular Element Method [87], which uses NURBS to represent the particles and, finally, the use of clusters or agglomeration of spheres [88]. The last one is chosen in this work since it provides great balance between shape representation accuracy1 and efficiency in terms of computational cost. Furthermore, it is the most versatile method in terms of particle shape and can naturally include angularities. The contact forces and torques are evaluated as usual on every sphere through eq. 4 and eq. 5. The contribution from every sphere is then gathered and translated to the centre of gravity of the rigid body altogether with the additional torque yielding from the application of the this force from the centre of every particle ${\displaystyle i}$ to the centre of the cluster ${\displaystyle {\mathbf {x} }_{cm}}$ through the distance vector ${\displaystyle {\mathbf {r} }_{i}^{p}={\boldsymbol {C}}_{i}-{\mathbf {x} }_{cm}}$.

 ${\displaystyle {\mathbf {F} }=\sum \limits _{i=1}^{np}{\mathbf {F} }_{i}}$ (72.a) ${\displaystyle {\mathbf {T} }=\sum \limits _{i=1}^{np}{\mathbf {T} }_{i}+\sum \limits _{i=1}^{np}{\boldsymbol {r}}_{i}^{p}\times {\mathbf {F} }_{i}}$ (72.b)
 Figure 21: Discretization of a rigid body using a cluster approach with spheres on the surface or overlapping in the interior

Once the total force ${\displaystyle {\mathbf {F} }}$ and torque ${\displaystyle {\mathbf {T} }}$ of the rigid body is obtained, the classical Newton's second law for the translation and the Euler rotation equations have to be solved in order to obtain the full motion of the rigid body (section 2.7.2). These equations can be integrated in an explicit way, preferably with a second or higher order scheme (section 2.7.3).

(1) The use of sphere cluster can introduce artificial friction due to the irregularities in the surface meshed by spheres. This problem is discussed in [89].

#### 2.7.1 Representation of the rotation

There are three ways which are very popular to represent rotations in the DEM: the use of Euler Angles, the use of rotation matrices and the use of quaternions. A review of the advantages and drawbacks of the methods can be found in [90].

The use of quaternions represents a clear advantage. It avoids the singularity problems that Euler angles present, it is more compact and it has less memory requirements than storing rotation matrices. Furthermore, the rotation operations are done in a more efficient way than using rotation matrices.

A rotation matrix ${\displaystyle R}$ is a ${\displaystyle 3\times {3}}$ orthogonal matrix which transforms a vector or a tensor from one coordinate system to another one as follows:

 ${\displaystyle {\mathbf {v} }'=R{\mathbf {v} }}$ (73.a) ${\displaystyle {\mathbf {A} }'=R{\mathbf {A} }R^{T}}$ (73.b)

Given a rotation of ${\displaystyle \theta }$ degrees over a unitary vector ${\displaystyle {\boldsymbol {u}}}$, the rotation matrix is constructed as follows:

 ${\displaystyle R={\begin{bmatrix}\cos \theta +u_{x}^{2}\left(1-\cos \theta \right)&u_{x}u_{y}\left(1-\cos \theta \right)-u_{z}\sin \theta &u_{x}u_{z}\left(1-\cos \theta \right)+u_{y}\sin \theta \\u_{y}u_{x}\left(1-\cos \theta \right)+u_{z}\sin \theta &\cos \theta +u_{y}^{2}\left(1-\cos \theta \right)&u_{y}u_{z}\left(1-\cos \theta \right)-u_{x}\sin \theta \\u_{z}u_{x}\left(1-\cos \theta \right)-u_{y}\sin \theta &u_{z}u_{y}\left(1-\cos \theta \right)+u_{x}\sin \theta &\cos \theta +u_{z}^{2}\left(1-\cos \theta \right)\end{bmatrix}}}$
(74)

A quaternion can summarize the same information just using 4 scalars. It is defined in the complex number system as:

 ${\displaystyle q=q_{0}+q_{1}{\boldsymbol {i}}+q_{2}{\boldsymbol {j}}+q_{3}{\boldsymbol {k}}}$
(75)

or in a compact form:

 ${\displaystyle q=[q_{0},\;{\boldsymbol {q}}]}$
(76)

Defining its conjugate as ${\displaystyle q^{*}=[q_{0},-{\boldsymbol {q}}]}$, the norm of a quaternion can be expressed:

 ${\displaystyle \lVert q\rVert ={\sqrt {qq^{*}}}}$
(77)

and its inverse:

 ${\displaystyle q^{-1}={\frac {q^{*}}{\lVert q\rVert }}}$
(78)

Now, given a rotation of ${\displaystyle \theta }$ degrees over a unitary vector ${\displaystyle {\boldsymbol {u}}}$, the resulting unit quaternion reads:

 ${\displaystyle q=\cos(\theta /2)+\sin(\theta /2){\boldsymbol {u}}}$
(79)

And the conversion from quaternions to a rotation matrix is the following:

 ${\displaystyle R={\begin{bmatrix}1-2(q_{2}^{2}+q_{3}^{2})&2q_{1}q_{2}-2q_{0}q_{3}&2q_{0}q_{2}+2q_{1}q_{3}\\2q_{1}q_{2}+2q_{0}q_{3}&1-2(q_{1}^{2}+q_{3}^{2})&2q_{2}q_{3}-2q_{0}q_{1}\\2q_{1}q_{3}-2q_{0}q_{2}&2q_{0}q_{1}+2q_{2}q_{3}&1-2(q_{1}^{2}+q_{2}^{2})\end{bmatrix}}}$
(80)

By using unit quaternions the intermediate transformation to a rotation matrix can be skipped and the rotation can be directly applied to vectors and tensors. The specification of unit quaternions is important in order to preserve lengths during rotational transformations. The rotations are applied in the following way:

 ${\displaystyle {\mathbf {v} }'=q{\mathbf {v} }q^{-1}}$ (81.a) ${\displaystyle {\mathbf {A} }'=\left(q\left(q{\mathbf {A} }q^{-1}\right)^{T}q^{-1}\right)^{T}}$ (81.b)

To do so, the multiplication operation needs to be employed. Given two quaternions ${\displaystyle p}$ and ${\displaystyle q}$ the multiplication yields a new quaternion ${\displaystyle t}$:

 ${\displaystyle t=pq=[p_{0}q_{0}-{\boldsymbol {p}}{\boldsymbol {q}},\;p_{0}{\boldsymbol {q}}+q_{0}{\boldsymbol {p}}+{\boldsymbol {p}}\times {\boldsymbol {q}}]}$
(82)

The vector involved in a quaternion multiplication (eq. 81.a) is treated as a quaternion ${\displaystyle {\boldsymbol {v}}=[0,\;{\boldsymbol {v}}]}$ with a null scalar part. The tensor multiplication (eq. 81.b) can be simply done treating the tensor as an assembly of vectors that are being multiplied subsequently. Note that the multiplication of quaternions is not commutative since it involves a cross product. A extended review on quaternion algebra can be found in [91].

#### 2.7.2 Rigid body dynamics

In a rigid body the distance between two material points is constant over time. Any spatial movement undergone by a rigid body can be described with the displacement of the centre of mass plus a rotation over some axis passing through the centre of gravity.
 Figure 22: A generic rigid body

For sake of convenience the spatial description of the body will be used identifying the position of every material point ${\displaystyle P}$ in time ${\displaystyle t}$ with its spatial position ${\displaystyle {\boldsymbol {x}}(t)}$ referred to global inertial reference system ${\displaystyle {\mathbf {X} },{\mathbf {Y} },{\mathbf {Z} }}$. In its turn, the superscript ${\displaystyle '}$ as in ${\displaystyle {\boldsymbol {x}}'(t)}$ denotes a quantity expressed with respect to the body fixed frame ${\displaystyle {\boldsymbol {x}}',{\boldsymbol {y}}',{\boldsymbol {z}}'}$. The temporal dependence will be dropped in the following developments for clarity.

The definition of the centre of mass of a body enclosed by the domain ${\displaystyle \Omega }$ supposing constant ${\displaystyle \rho }$ density is:

 ${\displaystyle {\mathbf {x} }_{cm}:={\frac {1}{m}}\int _{\Omega }{\rho \,{\mathbf {x} }\,d\Omega }}$
(83)

Defining ${\displaystyle {\mathbf {r} }:={\mathbf {x} }-{\mathbf {x} }_{cm}}$. The velocity and acceleration can be obtained:

 ${\displaystyle {\mathbf {{v}({\boldsymbol {x}})} }={\dot {\mathbf {x} }}_{cm}+{\boldsymbol {\omega }}\,\times \,{\boldsymbol {r}}}$ (84.a) ${\displaystyle {\mathbf {{a}({\boldsymbol {x}})} }={\ddot {\mathbf {x} }}_{cm}+{\dot {\boldsymbol {\omega }}}\,\times {\boldsymbol {r}}+{\boldsymbol {\omega }}\,\times \,({\boldsymbol {\omega }}\,\times \,{\boldsymbol {r}})}$ (84.b)

The linear and angular momentum are defined as:

 ${\displaystyle {\boldsymbol {L}}(t):=\int _{\Omega }{\rho \,{\mathbf {v} }\,d\Omega }}$ (85.a) ${\displaystyle {\boldsymbol {H}}(t):=\int _{\Omega }{{\boldsymbol {r}}\,\times \,\rho \,{\mathbf {v} }\,d\Omega }}$ (85.b)

and the balance expressions for linear and angular momentum read:

 ${\displaystyle {\dot {\boldsymbol {L}}}(t)={\mathbf {F} }(t)}$ (86.a) ${\displaystyle {\dot {\boldsymbol {H}}}(t)={\mathbf {T} }(t)}$ (86.b)

Now, the expression for the translational motion is obtained combining equation 85.a and 84.b onto the equation of balance of linear momentum 86.a yielding the classical Newton's second law of motion:

 ${\displaystyle {\begin{array}{l}{\mathbf {F} }={\boldsymbol {\dot {L}}}=m\,{\ddot {\mathbf {x} }}_{cm},\\\end{array}}}$
(87)

Likewise, the expression for the rotational motion is achieved plugging equation 85.a into equation 86.a and evaluating the temporal derivative. The expression of the Euler equations is found with the use of eq. 84.b and eq. 84.a onto the balance of linear momentum (eq. 86.a).

 ${\displaystyle {\mathbf {T} }={\boldsymbol {\dot {H}}}={\mathbf {I} }\cdot {\boldsymbol {\dot {\omega }}}+{\boldsymbol {\omega }}\times {\boldsymbol {I}}\cdot {\boldsymbol {\omega }}}$
(88)

Where ${\displaystyle {\mathbf {I} }}$ is the inertia tensor which is defined as:

 ${\displaystyle {\mathbf {I} }=\int _{\Omega }\rho \,\left(\left({\boldsymbol {r}}\cdot {\boldsymbol {r}}\right){\boldsymbol {1}}-{\boldsymbol {r}}\otimes {\boldsymbol {r}}\right)\,d\Omega }$
(89)

Note that the inertia tensor depends on the reference axis. Only in a body fixed frame the tensor has constant components. If we set this frame 88 in the so-called principal axis of inertia the tensor diagonalizes and the Euler equations can be expressed component-wise as:

 ${\displaystyle {T'}_{x}{=I'}_{x}\,{{\dot {\omega }}'}_{x}+(I'_{z}-I'_{y})\omega '_{z}\omega '_{y}}$ ${\displaystyle {T'}_{y}{=I'}_{y}\,{{\dot {\omega }}'}_{y}+(I'_{x}-I'_{z})\omega '_{x}\omega '_{z}}$ ${\displaystyle {T'}_{z}{=I'}_{z}\,{{\dot {\omega }}'}_{z}+(I'_{y}-I'_{x})\omega '_{y}\omega '_{x}}$
(90)

#### 2.7.3 Time integration of rotational motion in rigid bodies

The integration of the rotation needs different schemes than the ones presented for the translational motion in section 2.6.1 due to the higher complexity of the equations. The strategy described here is an adaptation of the scheme presented by Munjiza et al. [92] and Wellman [4]. The modification consists in the use of quaternions instead of rotational matrices in the integration scheme which makes the calculations more efficient in terms of computational cost and memory storage. The point of departure is the balance of angular momentum in the following form:

 ${\displaystyle {\mathbf {T} }={\frac {d{\boldsymbol {H}}}{dt}}}$
(91)

Munjiza et al. [92] introduced the idea that the change in angular momentum can be approximated by increments due to the change in the external torques at every time step. Actually, this assumption adapts perfectly to the temporal discretization used in DEM where the forces and torques are evaluated in discrete time steps.

 ${\displaystyle {\begin{array}{l}{\mathbf {H} }^{n+1}={\mathbf {H} }^{n}+\Delta t\,{\mathbf {T} }^{n}\\\end{array}}}$
(92)

This yields a constant angular momentum throughout a time step. The angular velocities can be approximated from the definition of the angular momentum expressed in the following way:

 ${\displaystyle {\mathbf {H} }={\mathbf {I} }\cdot {\boldsymbol {\omega }}}$
(93)

The key here is not to derive a constant angular velocity from the relation ${\displaystyle {\boldsymbol {\omega }}={\mathbf {I} }^{-1}\cdot {\mathbf {H} }}$ but approximate it using a higher order scheme such as a fourth-order Runge-Kutta.

Normally, the torques ${\displaystyle {\boldsymbol {T}}}$ and angular momentum ${\displaystyle {\boldsymbol {H}}}$ are expressed in global coordinates while the inertia tensor ${\displaystyle {\boldsymbol {I}}'}$ is naturally stored in the local body-fixed frame where it is diagonal with constant coefficients. Applying the quaternion tensor rotation described in equation 81.b, the local inertia tensor ${\displaystyle {\boldsymbol {I}}'}$ can be expressed in global coordinates, ${\displaystyle {\boldsymbol {I}}}$. Now, the angular velocity can be obtained from equation 93 as:

 ${\displaystyle {\begin{array}{l}{\boldsymbol {\omega }}=\left(\left(q\left(q\,{\mathbf {I} }'\,q^{-1}\right)^{T}q^{-1}\right)^{T}\right)^{-1}\cdot {\mathbf {H} }\\\end{array}}}$
(94)

where ${\displaystyle q}$ is the quaternion defining the transformation between local and global coordinates1. Instead of calculating it directly, a four order Runge-Kutta scheme is applied for the determination of an average angular velocity ${\displaystyle {\bar {\boldsymbol {\omega }}}}$ during the time step:

 ${\displaystyle {\boldsymbol {\omega }}_{1}:={\boldsymbol {\omega }}^{n}}$ (95.a) ${\displaystyle {\boldsymbol {\omega }}_{k}:=\left(\left(q_{k}\left(q_{k}\,{\mathbf {I} }'\,q_{k}^{-1}\right)^{T}q_{k}^{-1}\right)^{T}\right)^{-1}\cdot {\mathbf {H} }^{n+1}\qquad k\in [2,4]}$ (95.b) ${\displaystyle {\bar {\boldsymbol {\omega }}}:=1/6\,({\boldsymbol {\omega }}_{1}+2{\boldsymbol {\omega }}_{2}+2{\boldsymbol {\omega }}_{3}+{\boldsymbol {\omega }}_{4})}$ (95.c)

where the values of the transformation quaternions ${\displaystyle q_{k}}$ are:

 ${\displaystyle q_{2}:=q({\boldsymbol {\omega }}_{1},\Delta t/2)q^{n}}$ (96.a) ${\displaystyle q_{3}:=q({\boldsymbol {\omega }}_{2},\Delta t/2)q^{n}}$ (96.b) ${\displaystyle q_{4}:=q({\boldsymbol {\omega }}_{3},\Delta t)q^{n}}$ (96.c)

Once the average angular velocity during the time step ${\displaystyle {\bar {\boldsymbol {\omega }}}}$ is obtained, the final update predicts the velocity at the new step as:

 ${\displaystyle q^{n+1}=q({\bar {\boldsymbol {\omega }}},\Delta t)q^{n}}$ (97.a) ${\displaystyle {\boldsymbol {\omega }}^{n+1}=\left(\left(q^{n+1}\left(q^{n+1}\,{\mathbf {I} }'\,{(q^{n+1})}^{-1}\right)^{T}{(q^{n+1})}^{-1}\right)^{T}\right)^{-1}\cdot {\mathbf {H} }^{n+1}}$ (97.b)

The quaternions expressed in the form ${\displaystyle q({\boldsymbol {a}},b)}$ (eq. 96.a, 96.b, 96.c and 97.a ) represent incremental rotations that are derived from the application of constant angular velocities ${\displaystyle {\boldsymbol {a}}}$ during a fraction of time ${\displaystyle b}$. First, the corresponding rotation angles are calculated:

 ${\displaystyle \Delta {\boldsymbol {\theta }}({\boldsymbol {a}},b)=b\cdot {\boldsymbol {a}}}$
(98)

The unitary vector defining the rotation is ${\displaystyle {\boldsymbol {u}}_{\theta }=\Delta {\boldsymbol {\theta }}\,(\lVert \Delta {\boldsymbol {\theta }}\rVert )^{-1}}$ and its magnitude ${\displaystyle \lVert \Delta {\boldsymbol {\theta }}\rVert }$. With these two quantities the mapping ${\displaystyle \Delta {\boldsymbol {\theta }}({\boldsymbol {a}},b)\rightarrow q({\boldsymbol {a}},b)}$ can be achieved applying equation 79.

#### Direct explicit integration

Some codes perform a direct forward explicit integration of the equations of motion which is presented here. Eq. 11.3.1 is expressed in a diagonalized local frame where the ${\displaystyle '}$ superscript has been dropped for clarity:

 ${\displaystyle \omega _{x}^{n+1}=\omega _{x}^{n}+{\frac {\Delta t}{I_{x}}}\left(T_{x}^{n}-(I_{z}-I_{y})\,\omega _{z}^{n}\,\omega _{y}^{n}\right)}$ (99.a) ${\displaystyle \omega _{y}^{n+1}=\omega _{y}^{n}+{\frac {\Delta t}{I_{y}}}\left(T_{y}^{n}-(I_{x}-I_{z})\,\omega _{x}^{n}\,\omega _{z}^{n}\right)}$ (99.b) ${\displaystyle \omega _{z}^{n+1}=\omega _{z}^{n}+{\frac {\Delta t}{I_{z}}}\left(T_{z}^{n}-(I_{y}-I_{x})\,\omega _{y}^{n}\,\omega _{x}^{n}\right)}$ (99.c)

#### Rotation integration benchmark

In the works of Munjiza et al. [92] and Lillie [93] an example which can be analytically solved is run with the presented scheme using rotation matrices instead of quaternions. They showed that the scheme rapidly yields accurate results. Here the same example is reproduced to check the good implementation of the ${\displaystyle RK-4}$ scheme using quaternions and also to show its superiority against a direct explicit integration.

A cylinder of ${\displaystyle 1.5\,m}$ height and ${\displaystyle 0.5\,m}$ radius with a density of ${\displaystyle 1\,kg/m3}$ is set to freely rotate in the space with an initial angular velocity of ${\displaystyle {\boldsymbol {\omega }}_{0}=[0,1,100]\;rad/s}$ during ${\displaystyle 0.5}$ seconds. Since the initial axis of rotation does not coincide with any of the principal directions ${\displaystyle {{\boldsymbol {e}}'}_{1}}$, ${\displaystyle {{\boldsymbol {e}}'}_{2}}$, ${\displaystyle {{\boldsymbol {e}}'}_{3}}$, the resulting rotational motion presents the so called torque free precession which is characterized by a varying rotational velocity ${\displaystyle {\boldsymbol {\omega }}}$ and inertia tensor ${\displaystyle {\boldsymbol {I}}}$ (in global coordinates).

As figures 24 and 25 show, the ${\displaystyle RK-4}$ scheme is much more accurate than the direct integration. Even using a time step ten times smaller for the direct integration than for the ${\displaystyle RK-4}$, the last scheme performs better. Both methods proved to have convergence to the analytical solution when smaller time steps were used.
 Figure 23: Cylinder set-up
 Figure 24: Integration results for local ${\displaystyle {\boldsymbol {\omega }}_{x}}$
 Figure 25: Integration results for local ${\displaystyle {\boldsymbol {\omega }}_{y}}$

We introduced the implementation of the ${\displaystyle RK-4}$ with quaternions in order to have a scheme that is much more efficient in computational cost compared to the original one using rotation matrices and it handles the storage of the rotations with less than half of the memory. Therefore, and taking into account the poor accuracy of the direct approximation, we highly recommend the use of the ${\displaystyle RK-4}$ method for the integration of the rotations both for spherical and non-spherical particles.

(1) It shall be noted that in the case of spherical particle we can skip this transformation since the inertia is diagonal and constant in every reference system.

### 2.8 Mesh generation

Several industrial processes in which the particle flow is simulated do not require an initial mesh but an inlet and possibly an outlet. However, in a general case, an initial configuration of particles is required and a thus a generation tool has to be employed. Normally, a heterogeneous mesh is desired with a specific granulometry or size distribution. To that end several techniques exist which are based on different principles.

A first family of methods, known as Lily-pound methods [94,95,96] insert particles in random locations checking if intersections occur, if so, a new location is determined. On the other hand, the advancing front techniques [97,98,99] collocate the particles layer by layer starting from the boundaries or the interior of the domain presenting a better control on the desired size distribution. Different modifications exist which attempt to improve the packing of these techniques like in [100].

In the framework of the thesis the GiD sphere mesher developed by Labra [101] has been used for the generation of the sphere meshes. Its principle is based on a first collocation of particles with a later rearrangement technique [102,103] which corrects the inclusions generated being able to achieve dense packings. The reduction of the porosity is solved with the minimization of a distance function with every particle and its neighbours.

Some other meshing techniques rely on a DEM pre-simulation to fill the domain with an inlet or pushing boundaries in a expanded domain where the particles are initially set. These are the techniques used for the generation of meshes based on clusters of spheres in the framework of this thesis.

### 2.9 Basic DEM flowchart

 Figure 26: Basic DEM flowchart

## 3. The Double Hierarchy (${\displaystyle H^{2}}$) Method for DE-FE contact detection3. The Double Hierarchy ( H 2 {\displaystyle H^{2}} ) Method for DE-FE contact detection

This chapter presents a detailed description of the contact detection between discrete elements and finite elements. First, the state of the art of the existing methods for modelling the contact with boundaries is reviewed as well as the specific DE-FE collision detection methods. Later, the Double Hierarchy Method [104], a novel method developed for the interaction with rigid structures, is thoroughly described including implementation details together with validation examples.

As it will be shown, the literature lacks of a flexible method that computes efficiently the contact between particles and FE, allowing for multi-contact problems and providing continuity of forces in non-smooth contact regions. The objective of this new method is to provide a robust, versatile and efficient procedure which can tackle the above-mentioned problems and be implemented in any DEM code allowing parallel computation.

The method presented here adapts perfectly to the case of spherical particles (including clusters of spheres) contacting triangles or quadrilaterals belonging to the rigid boundaries included in a DEM simulation. The discussion on how to upgrade this method to the case of deformable structures will be presented in chapter 4.

### 3.1 State of the art

Several solutions have been reported for the inclusion of boundaries to the discrete element method. Among the simplest ones is the glued-sphere approach [105], which approximates any complex geometry (i.e. a rigid body or boundary surface) by a collection of spherical particles so it retains the simplicity of particle-to-particle contact interaction. This approach, however, is geometrically inaccurate and computationally intensive due to the introduction of an excessive number of particles. A second simple approach (used in some numerical codes, e.g., ABAQUS) is to define the boundaries as analytical surfaces. This approach is computationally inexpensive, but it can only be applied in certain specific scenarios, where the use of infinite surfaces does not disturb the calculation. A more complex approach which combines accuracy and versatility is to resolve the contact of particles (spheres typically) with a finite element boundary mesh. These methods take into account the possibility of contact with the primitives of the FE mesh surface, i.e., facet, edge or vertex contact. The term FE will be used in this dissertation when referring to the geometry elements (triangles, quadrilaterals, etc.) which are used to discretize the boundaries even if they are not used for the calculation of a deformable solid.

Horner et al. [14] and Kremmer and Favier [80] developed the first hierarchical contact resolution algorithms for contact problems between spherical particles and triangular elements, while Zang et al. [106] proposed similar approaches accounting for quadrilateral facets. Dang and Meguid [27] upgraded the method introducing a numerical correction to improve smoothness and stability. Su et al. [107] developed a complex algorithm involving polygonal facets under the name of RIGID which includes an elimination procedure to resolve the contact in different non-smooth contact situations. This approach, however, does not consider contact with entities of different surfaces at the same time (multiple contacts) leading to an inaccurate contact interaction. The upgraded RIGID-II method presented later by Su et al. [108] and also the method proposed by Hu et al. [109] account for the multiple contact situations, but they have a complex elimination procedure with many different contact scenarios to distinguish, which is difficult to code in practice. Chen et al. [110] presented a simple and accurate algorithm which covers many situations. Their elimination procedure, however, requires a special database which strongly limits the parallel computation.

In the framework of this monograph, the Double Hierarchy Method (${\textstyle H^{2}}$) [104], has been developed. It consists in a simple contact algorithm based on the FE boundary approach. It is specially designed to resolve efficiently the intersection of spheres with triangles and planar quadrilaterals but it can also work fine with any other higher order planar convex polyhedra. A two layer hierarchy is applied upgrading the classical hierarchy method presented by Horner [14]; namely hierarchy on contact type followed by hierarchy on distance. The first one, classifies the type of contact (facet, edge or vertex) for every contacting neighbour in a hierarchical way, while the distance-based hierarchy determines which of the contacts found are valid or relevant and which ones have to be removed.

Industrial applications may involve a large number of particles and also a fine definition of the boundaries which, using boundary FE, would turn into large number of conditions to check. The selected algorithm works efficiently in parallel computations as will be shown in chapter 6. This is a clear advantage over the above-mentioned publications which algorithms are mostly serial. Exceptions are Nakashima [23] whose method is presumably parallelizable and Zang [106] and Su [108] which remark the importance of the future parallelization of their algorithms.

Summarizing, the contact search framework presented is designed to satisfy the following requirements:

1. Include poly-disperse elements for both: FEs and DEs.
2. Allow different FE geometries and primitives (triangle, quadrilateral, polygon).
3. Ensure contact continuity in non-smooth regions (edges and vertices).
4. Resolve multiple contacts and contact with different entities simultaneously.
5. Need low memory storage.
6. Be simple, fast and accurate.
7. Be fully parallelizable.

Table 5 summarizes the strengths and drawbacks of the reviewed contact detection methods. Methods which have a elimination procedure to remove the invalid contacts (RIGID-II [108], Hu et al. [109], Chen et al. [110] and ${\textstyle H^{2}}$) are the most accurate. They treat the cases with large indentations (relative to the FE size) and provide a solution to the contact continuity in non-smooth boundary regions. These methods have, however, some limitations due to the fact that the real deformed geometry of the sphere is not represented in the DEM. Due to this fact, error in the contact detection in concave transitions is common for all these methods (including the ${\textstyle H^{2}}$). This is analysed in section 3.5.

 Glued Anal. Hierarchy RIGID RIGID-II Hu Chen ${\displaystyle H^{2}}$ [105] [14,80,106,27] [107] [108] [109] [110] [104] Wide size rate DEs/FEs - - ${\displaystyle \times }$ ✓ ✓ ✓ ${\displaystyle \times }$ ✓ Contact elem. typologies ${\displaystyle \times }$ - ✓ ✓ ✓ ${\displaystyle \times }$ ${\displaystyle \times }$ ✓ Boundary shape variety ✓ ${\displaystyle \times }$ ✓ ✓ ✓ ✓ ✓ ✓ Multi-contact ✓ - ✓ ${\displaystyle \times }$ ✓ ✓ ✓ ✓ Simple ✓ ✓ ✓ ${\displaystyle \times }$ ${\displaystyle \times }$ ${\displaystyle \times }$ ✓ ✓ Efficient ${\displaystyle \times }$ ✓ ${\displaystyle \times }$ ✓ ${\displaystyle \times }$ ${\displaystyle \times }$ ✓ ✓ Accurate ${\displaystyle \times }$ ${\displaystyle \times }$ ✓ ${\displaystyle \times }$ ✓ ✓ ${\displaystyle \times }$ ✓ Low storage ✓ ✓ ${\displaystyle \times }$ ${\displaystyle \times }$ ${\displaystyle \times }$ ${\displaystyle \times }$ ✓ ✓ Upgradable to CSM ${\displaystyle \times }$ ${\displaystyle \times }$ ✓ ✓ ✓ ✓ ✓ ✓ Large indentation ${\displaystyle \times }$ ✓ ${\displaystyle \times }$ ${\displaystyle \times }$ ✓${\textstyle \ast }$ ✓${\textstyle \ast }$ ${\displaystyle \times }$ ✓${\textstyle \ast }$ Contact continuity ${\displaystyle \times }$ - ✓${\textstyle \ast }$ ${\displaystyle \times }$ ✓${\textstyle \ast }$ ✓${\textstyle \ast }$ ✓${\textstyle \ast }$ ✓${\textstyle \ast }$

Symbol (✓) implies that the method satisfies the property while (${\textstyle \times }$) indicates that the method does not satisfy the property. Symbol (-) denotes that the property does not apply to that method and (✓${\textstyle \ast }$) means that, the method satisfies the property upon some limitations.

### 3.2 DE-FE contact detection algorithm

The strategy of dividing the search into global and local stages also applies to the DE-FE collision detection. In the same way, the methods described in section 2.2 regarding the global search can be also used here. The cell-based algorithm presented in [40] has been selected for the global search due to its simplicity and the possibility to be parallelized.

As it has been appointed in section 2.2, the most expensive part of the collision detection lies on the local resolution which can reach values over 75 percent of the simulation when non-spherical elements are involved [14]. To that end, a specialized algorithm has been developed for the case of collision between spheres (particles) and triangles or quadrilaterals (boundary elements) which is particularly efficient. Moreover, a further split of the Local Contact Resolution is performed: a) A Fast Intersection Test, b) Full characterization of valid contacts. Figure 29 shows the different stages of the search.

#### 3.2.1 Global Search algorithm

The main purpose of the Global Search is to determine through a fast rough search which are the potential neighbours for every element in the domain. A cell-based algorithm [40] is chosen here which has been parallelized in OMP and adapted for the DE-FE search. The FE domain is selected to build the search bins taking advantage of the fact that usually the spatial distribution of the FEs is more regular and in some cases fixed. As an additional feature, the Search Bins is built dynamically considering only the elements belonging to the intersection of the bounding boxes of the DEs and FEs domains, ${\textstyle FE\in \Omega _{I}}$ and ${\textstyle DE\in \Omega _{I}}$. Fig. 27a shows how the intersection evolves as long as the simulation goes on. On the other hand, only the DEs inside the intersection domain (${\textstyle \Omega _{I}}$) will look for their neighbours. This reduces significantly the contact pairs to be checked afterwards and, therefore, the global search performance is increased.

In the global search, every FE and DE has an associated Bounding Box (${\textstyle FE_{BBX},DE_{BBX}}$) that is used to tag the position of the elements on the Search Bins and rapidly check for potential neighbours. This is done using a hash table structure as depicted in fig. 27d which relates each cell to the bounding box ${\textstyle FE_{BBX}}$ that fall into it. Rectangular hexahedral bounding boxes encompassing both types of elements are chosen here.

The steps needed to perform the neighbouring search at the Global Search level are:

1. Set the bounding box of the intersection of domains ${\textstyle \Omega _{I}}$ (fig. 27a).
2. Set the bounding box for every ${\textstyle FE\in \Omega _{I}}$ (fig. 27b).
3. Generate the Search Bins based on the size and position of the bounding boxes ${\textstyle FE_{BBX}}$ of the ${\textstyle FE\in \Omega _{I}}$ (fig. 27c).
4. Place every FE in the Search Bins (based on their associated bounding box ${\textstyle FE_{BBX}}$ coordinates) and build the hash table (fig. 27d).
5. Set the bounding box for every ${\textstyle DE\in \Omega _{I}}$ (fig. 27e).
6. For every DE particle ${\textstyle \in \Omega _{I}}$ obtain the FE potential neighbours in the Search Bins. Check the intersection of the ${\textstyle DE_{BBX}}$ with the ${\textstyle FE_{BBX}}$ of the FEs lying in the surrounding cells (fig. 27f).
7. Apply the Local Resolution Method to the pairs with intersecting bounding boxes (fig. 27g).
 (a) Evolution of ${\textstyle \Omega _{I}}$ (b) ${\displaystyle FE_{BBX}\in \Omega _{I}}$ (c) Bins over FEs ${\textstyle \in \Omega _{I}}$ (d) Hash table (e) ${\displaystyle DE_{BBX}\in \Omega _{I}}$ (f) Intersection cells (g) Local Contact Resolution Figure 27: Global search stages

#### 3.2.2 Local Contact Resolution

Normally a full characterization takes place after the global search and determines completely the contact status of each potential contact pair. In this monograph a split is suggested:

• Fast Intersection Test: First, the actual contacting pairs are determined. This has to be fast because there are many FE potential neighbours in the adjacent cells to be checked. Therefore, all detailed contact computations such as determining the type of contact, the contact point, normal direction, etc. are skipped. On the other hand, a good accuracy in the determination of the contacting neighbours is needed. It should be avoided to fill the contact pool with FE which do not have contact and have to be eliminated or treated subsequently. This procedure is described in detail in section 3.3.
• Full contact characterization: A more expensive check takes place which determines the type of contact of every neighbour, which are the relevant contacts and which ones have to be removed in order to avoid instabilities or redundant contact evaluations in non-smooth regions and contact transitions. All the detailed contact characteristics are fully determined at this stage for each one of the valid neighbouring entities.

The split gives the code higher modularity, i.e. any other contact characterization can be applied for the contacting entities. Moreover, in the in-house code Kratos, the split yields also higher efficiency (see table 24 in chapter 6). This is due to the fact that the full characterization is a much more expensive procedure than the simple Fast Intersection Test, and at the same time, the first group of FE potential neighbours is very large in comparison to the group of FE with contact.

In order to demonstrate this, an example of a horizontal mixer with approximately ${\textstyle 30\,k}$ DEs and ${\textstyle 2.5\,k}$ FEs has been run for ${\textstyle 0.5}$ second, i.e. ${\textstyle 1.5}$ turns of the helical blades (full description in section 6.2.2). The cumulative counts of the following quantities is computed:

1. FE Potential neighbours: The number of times the Fast Intersection Test (section 3.3) is called (number of FE potential neighbours to be checked) averaged over the number of particles.
2. FE with contact: The average number of FE per particle that yield a positive result (have intersection with sphere) in the Fast Intersection Test.
3. Entity with valid contact: The average number of relevant entities per particle determined by the ${\textstyle H^{2}}$ Method.
 Figure 28: Counts of FE checks in different stages

Fig. 28 presents the results which show that the number of FE Potential Neighbours to be treated is large compared to the FE with actual contact, a ratio of ${\textstyle 30:1}$. Additionally, as it will be shown in chapter 6, the improvement in performance showed in Table 24 it can be concluded that it is a good choice to perform the split which additionally brings modularity to the code.

Figure 29 summarizes the stages in which the neighbour finding is divided.
 Figure 29: Neighbour finding scheme

### 3.3 Fast Intersection Test

An efficient algorithm designed to determine the intersection of spheres contacting triangles or planar quadrilaterals is described here. Some of the procedures existing in the computer graphics bibliography [111,112] have been adapted to the case where the facet contact (inside of the FE) occurs in a substantial higher frequency compared to edge and vertex geometrical contact types. See [109] where the type of contact frequency (facet, edge, vertex) is determined for different number of particles and relative sizes.

The test works for any planar convex polygons of N sides. For every DE ${\textstyle \in \Omega _{I}}$ we loop over the FE potential neighbours provided by Global Neighbour Searching algorithm. Every FE with valid contact is stored in an array for every DE.

#### 3.3.1 Intersection test with the plane containing the FE

The first check is to determine whether the particle intersects the ${\textstyle \pi ^{m}}$ plane formed by the ${\textstyle m-th}$ planar finite element ${\textstyle {\mbox{(e)}}^{m}}$. This is represented in fig. 30.

 Figure 30: Intersection of a DE particle with a plane formed by a plane FE

The outward-pointing normal of the plane can be calculated with the cross product ${\textstyle {\boldsymbol {T}}}$ of any pair of edges taken counter-clockwise. This can be written in the following form, using the permutation tensor ${\textstyle \epsilon _{ijk}}$ on two edges formed, for example, by the three consecutive vertices ${\textstyle {\boldsymbol {\mathrm {v} }}^{1}}$, ${\textstyle {\boldsymbol {\mathrm {v} }}^{2}}$, ${\textstyle {\boldsymbol {\mathrm {v} }}^{3}}$:

 ${\displaystyle T_{i}=\epsilon _{ijk}(\mathrm {v} _{j}^{2}-\mathrm {v} _{j}^{1})\cdot (\mathrm {v} _{k}^{3}-\mathrm {v} _{k}^{2})}$
(100)

which has to be normalized to unit length to obtain the normal to the plane ${\textstyle {\boldsymbol {n}}}$:

 ${\displaystyle {\boldsymbol {n}}={\frac {\boldsymbol {T}}{\lVert {\boldsymbol {T}}\rVert }}}$
(101)

In the case of a zero-thickness wall which can have contact at both sides of the FE, the sense of the normal will be set such that points outwards to each particle centre. Once the normal is defined, the distance of the DE centre ${\textstyle {\boldsymbol {C}}}$ to the plane ${\textstyle \pi ^{m}}$ can be determined taking any known point of the plane, namely a vertex ${\textstyle \mathrm {v} ^{a}}$, as

 ${\displaystyle d_{\pi }=\sum \limits _{i=1}^{3}\left(n_{i}\cdot C_{i}-n_{i}\cdot \mathrm {v} _{i}^{a}\right)}$
(102)

The distance ${\textstyle d_{\pi }}$ should be compared to the radius ${\textstyle R}$. If and only if ${\textstyle \left|d_{\pi }\right|\leq R}$, the contact between the sphere and the FE is possible. In this case, we proceed with the next checks. Otherwise, the contact with the current FE is discarded and we will jump to check the next potential FE neighbour.

#### 3.3.2 Inside-Outside test

The purpose of this test is to determine whether the contact is inside the FE (facet contact) or outside (edge, vertex or no contact). It applies to the cases which ${\textstyle \left|d_{\pi }\right|\leq R}$. A modification of the Inside-Outside status check [113] is used. The projection ${\textstyle C_{\pi ^{m}}}$ of the centre ${\textstyle {\boldsymbol {C}}}$ of a DE onto the plane ${\textstyle \pi ^{m}}$ formed by an element ${\textstyle {\mbox{(e)}}^{m}}$ with normal ${\textstyle n}$ can be calculated as

 ${\displaystyle \mathbf {C} _{\pi ^{m}}=\mathbf {C} -d_{\pi }\cdot \mathbf {n} }$
(103)

The next step is to evaluate whether the projection ${\textstyle C_{\pi ^{m}}}$ lies inside or outside the FE ${\textstyle {\mbox{(e)}}^{m}}$ with respect to every edge ${\textstyle {\boldsymbol {e}}^{a}}$ formed with the vertices ${\textstyle {\boldsymbol {\mathrm {v} }}^{a}}$ and ${\textstyle {\boldsymbol {\mathrm {v} }}^{a+1}}$ (${\textstyle {\boldsymbol {\mathrm {v} }}^{N}={\boldsymbol {\mathrm {v} }}^{0}}$) (See fig. 31). For every edge ${\textstyle {\boldsymbol {e}}^{a}}$ we compute the cross product sign ${\textstyle s^{a}}$ as

 ${\displaystyle {\boldsymbol {e}}^{a}={\boldsymbol {\mathrm {v} }}^{a+1}-{\boldsymbol {\mathrm {v} }}^{a}}$
(104)

 ${\displaystyle s^{a}=\left({\boldsymbol {e}}^{a}\times ({\boldsymbol {C}}_{\pi ^{m}}-{\boldsymbol {\mathrm {v} }}^{a})\right)\cdot \mathbf {n} }$
(105)
If the product is positive, the projection point ${\textstyle C_{\pi ^{m}}}$ turns to be inside the triangle with respect to that edge. The loop proceeds with the next edges. If the same result is found for every edge, contact occurs with the facet of the FE (Inside) and so the contact is assured. Otherwise, if for any edge an Outside status is found, the loop aborts automatically and no contact with facet can be found. The current value of the edge index ${\textstyle a}$ is stored in an auxiliary variable ${\textstyle index_{e}}$ which will be used in the next step where contact with vertices or edges is checked.
 Figure 31: Inside-Outside check of the projection point edge by edge

#### 3.3.3 Intersection test with an edge

This test is needed for the cases where ${\textstyle \left|d_{\pi }\right|\leq R}$ but the Inside-Outside test failed. Here we use the idea that the edge contact can not happen to be on the edges where the Inside-Outside check yield a Inside status. Therefore, it is recommendable to test the edges ${\textstyle {\boldsymbol {e}}^{a}}$ with ${\textstyle a\in \left[{index}_{e},N\right]}$ starting from the vertex which failed in the previous test and skipping the previous ones (Note that the edge check is the most expensive one). This approach has also been used by Chen et al. [110].

First, the shortest distance ${\textstyle d_{e}}$ between the edge ${\textstyle {\boldsymbol {e}}^{a}}$ and the particle centre ${\textstyle {\boldsymbol {C}}}$ should be calculated and compared to the radius ${\textstyle R}$. The distance is calculated finding out the contact point ${\textstyle {\boldsymbol {Pc}}}$, as

 ${\displaystyle d_{e}=\lVert {\boldsymbol {Pc}}-{\boldsymbol {C}}\rVert }$
(106)

 ${\displaystyle {\boldsymbol {Pc}}={\boldsymbol {\mathrm {v} }}^{a}+p{\frac {{\boldsymbol {e}}^{a}}{\lVert {\boldsymbol {e}}^{a}\rVert }}}$
(107)

 ${\displaystyle {\boldsymbol {e}}^{a}={\boldsymbol {\mathrm {v} }}^{a+1}-{\boldsymbol {\mathrm {v} }}^{a}}$
(108)

where ${\textstyle p}$ is the distance resulting from the projection of the vector connecting the centre ${\textstyle {\boldsymbol {C}}}$ and the vertex ${\textstyle {\boldsymbol {\mathrm {v} }}^{a}}$ onto the edge ${\textstyle {\boldsymbol {e}}^{a}}$:

 ${\displaystyle p=({\boldsymbol {C}}-{\boldsymbol {\mathrm {v} }}^{a})\cdot {\boldsymbol {e}}^{a}}$
(109)
 Figure 32: Intersection of a DE particle with an edge

If ${\textstyle d_{e}>R}$ the contact with this edge is not possible and the check starts again with the next edge ${\textstyle {\boldsymbol {e}}^{a+1}}$. Otherwise, if ${\textstyle d_{e}\leq R}$ we determine where the ${\textstyle {\boldsymbol {Pc}}}$ lies, along the edge, with the help of ${\textstyle \eta }$, defined as:

 ${\displaystyle \eta ={\frac {p}{\lVert {\boldsymbol {e}}^{a}\rVert }}}$
(110)

The case of ${\textstyle 0\leq \eta \leq 1}$ implies edge contact. Therefore contact is found and the Fast Intersection Test finishes yielding a positive result. The FE neighbour is saved to the current DE and the algorithm proceeds to check the next FE potential neighbour.

Otherwise, if this test failed for the current edge ${\textstyle {\boldsymbol {e}}^{a}}$, the connecting vertices (${\textstyle {\boldsymbol {\mathrm {v} }}^{a}}$ and ${\textstyle {\boldsymbol {\mathrm {v} }}^{a+1}}$) have to be evaluated. A value of ${\textstyle \eta <0}$ indicates that the check has to be done with ${\textstyle {\boldsymbol {\mathrm {v} }}^{a}}$; on the other hand, for ${\textstyle \eta >1}$ the vertex to be tested is ${\textstyle {\boldsymbol {\mathrm {v} }}^{a+1}}$.

#### 3.3.4 Intersection test with a vertex

For the vertex ${\textstyle {\boldsymbol {\mathrm {v} }}^{a}}$ under consideration the squared distance to the DE centre ${\textstyle {\boldsymbol {C}}}$ is calculated:

 ${\displaystyle {d_{{\boldsymbol {\mathrm {v} }}^{a}}}^{2}=\sum \limits _{i=0}^{i<3}\left({\boldsymbol {C_{i}}}-{\boldsymbol {\mathrm {v} }}_{i}^{a}\right)^{2}}$
(111)

If ${\textstyle {d_{\mathrm {v} ^{a}}}^{2}\leq R^{2}}$, then the Fast Intersection Test yields a positive result and the test finishes. Otherwise, the test moves on with the next edge ${\textstyle {\boldsymbol {e}}^{a+1}}$ and its subsequent vertices.

We recall that the purpose of this Fast Intersection Test is merely to determine whether there is intersection or not between the DE sphere and the FE planar convex polygon. An intersection found with a vertex or edge does not ensure that this is the actual contact point. In this case, however, we omit at this stage further checks with subsequent edges or vertices where the contact point can happen to be closer.

#### 3.3.5 Fast Intersection Test algorithm

 Parallel loop over all DE, check FE potential neighbours.
 ${\textstyle (1)}$ Intersection with plane containing the FE ${\displaystyle {\mbox{(e)}}^{m}}$ Calculate normal outwards ${\textstyle {\boldsymbol {n}}={\frac {\boldsymbol {T}}{\lVert {\boldsymbol {T}}\rVert }}}$, ${\textstyle \,T_{i}=\epsilon _{ijk}(\mathrm {v} _{j}^{2}-\mathrm {v} _{j}^{1})\cdot (\mathrm {v} _{k}^{3}-\mathrm {v} _{k}^{2})}$. Calculate distance to plane ${\textstyle d_{\pi }=\sum \limits _{i=1}^{3}\left(n_{i}\cdot C_{i}-n_{i}\cdot \mathrm {v} _{i}^{a}\right)}$. if( ${\textstyle \left|d_{\pi }\right|>R}$ ): ${\textstyle \Rightarrow }$ Go to ${\textstyle (4)}$ (${\textstyle False}$). else: ${\textstyle \Rightarrow }$ Calculate ${\textstyle \mathbf {C} _{\pi ^{m}}=\mathbf {C} -d\cdot \mathbf {n} }$ and Go to ${\textstyle (2)}$.
 ${\textstyle (2)}$ Inside-Outside test Initialize ${\textstyle {index}_{e}=0}$ and Inside-Outside flag = In. loop over every edge ${\textstyle {\boldsymbol {e}}^{a}={\boldsymbol {\mathrm {v} }}^{a+1}-{\boldsymbol {\mathrm {v} }}^{a}}$ with ${\textstyle a\in \left[0,N\right]}$. calculate ${\textstyle s^{a}=\left({\boldsymbol {e}}^{a}\times ({\boldsymbol {C}}_{\pi ^{m}}-{\boldsymbol {\mathrm {v} }}^{a})\right)\cdot \mathbf {n} }$. if(${\textstyle s^{a}<0}$): ${\textstyle \Rightarrow }$ Inside-Outside = Out. Break loop. Save ${\textstyle {index}_{e}=a}$. Go to ${\textstyle (3)}$. else(${\textstyle s^{a}\geq 0}$): ${\textstyle \Rightarrow }$ Continue with next edge. if(Inside-Outside flag == In): ${\textstyle \Rightarrow }$ Go to ${\textstyle (4)}$ (${\textstyle True}$). else: ${\textstyle \Rightarrow }$ Go to ${\textstyle (3)}$.
 ${\textstyle (3)}$ Intersection with Edge and Vertex loop over every edge ${\textstyle {\boldsymbol {e}}^{a}}$ with ${\textstyle a\in \left[{index}_{e},N\right]}$. Calculate projection: ${\textstyle p=({\boldsymbol {C}}-{\boldsymbol {\mathrm {v} }}^{a})\cdot {\boldsymbol {e}}^{a}}$. Calculate the contact point: ${\textstyle {\boldsymbol {Pc}}={\boldsymbol {\mathrm {v} }}^{a}+p{\frac {{\boldsymbol {e}}^{a}}{\lVert {\boldsymbol {e}}^{a}\rVert }}}$. Calculate distance to edge ${\textstyle d_{\boldsymbol {e}}=\lVert {\boldsymbol {Pc}}-{\boldsymbol {C}}\rVert }$. if(${\textstyle d_{e}>R}$): ${\textstyle \Rightarrow }$ Continue with next edge. else: Calculate ${\textstyle \eta ={\frac {p}{\lVert {\boldsymbol {e}}^{a}\rVert }}}$. if(${\textstyle 0\leq \eta \leq 1}$ ): ${\textstyle \Rightarrow }$ Go to ${\textstyle (4)}$ (${\textstyle True}$). if(${\textstyle \eta <0}$): ${\textstyle \Rightarrow }$ ${\textstyle d_{{\mathrm {v} }^{a}}^{2}=\sum \limits _{i=0}^{i<3}\left({\boldsymbol {C_{i}}}-{\boldsymbol {\mathrm {v} }}_{i}^{a}\right)^{2}}$. if(${\textstyle d_{{\mathrm {v} }^{a}}^{2}\leq R^{2}}$): ${\textstyle \Rightarrow }$ Go to ${\textstyle (4)}$ (${\textstyle True}$). else: ${\textstyle \Rightarrow }$ check next edge. if(${\textstyle \eta >1}$): ${\textstyle \Rightarrow }$ ${\textstyle d_{{\mathrm {v} }^{a+1}}^{2}=\sum \limits _{i=0}^{i<3}\left({\boldsymbol {C_{i}}}-{\boldsymbol {\mathrm {v} }}_{i}^{a+1}\right)^{2}}$. if(${\textstyle d_{{\mathrm {v} }^{a+1}}^{2}\leq R^{2}}$): ${\textstyle \Rightarrow }$ Go to ${\textstyle (4)}$ (${\textstyle True}$). else: ${\textstyle \Rightarrow }$ check next edge. Go to ${\textstyle (4)}$ (${\textstyle False}$).
 ${\textstyle (4)}$ Contact Found (${\displaystyle True/False}$) ${\textstyle True}$: ${\textstyle \Rightarrow }$ Store ${\textstyle {\mbox{(e)}}^{m}}$ as FE with contact and Continue. ${\textstyle False}$: ${\textstyle \Rightarrow }$ Stop! No contact.

The presented algorithm applies to any planar convex polygons of N sides.

### 3.4 The Double Hierarchy Method

The application of constitutive contact laws such as the Hertz-Mindlin (section 2.5.2) requires that the contact surfaces are smooth and present a unique normal at each point. In the DE-FE contact, usually, the original geometry presents regions where this requirement is not fulfilled. Moreover, even the smooth surfaces loose this feature when they are discretized by means of FEs. In these situations a special treatment of the non-smooth regions should be applied under the requirement of some conditions to ensure reasonable results. The following conditions were also analysed in the work by Wellmann [4]:

• The contact constitutive model will be applied normally when the contact is on the facet and will vanish when there is no interpenetration between the elements.
• There should be no discontinuities in the contact force when a contact point evolves from facet to edge and the other way round in order to avoid non-physical results and numerical instabilities.
• The energy should be conserved in an elastic frictionless impact.

The use of the present contact determination algorithm helps the selected contact model ensuring these objectives as it will be shown through the validation examples in section 3.6.

This procedure is applied to the list of FE with contact that the Fast Intersection Test has generated for every particle. In the case of no previous fast check this operation could be directly applied as a Local Contact Resolution with the disadvantage that many potential FE have to be tested. It is developed in two different stages:

• Contact Type Hierarchy (section 3.4.1): where for every FE with contact the entity with higher priority is determined.
• Distance Hierarchy (section 3.4.2): the elimination procedure takes place determining which contact points have distance priority over others which are redundant or false and have to be eliminated.

#### 3.4.1 Contact Type Hierarchy

The basis of this procedure is that each primitive has hierarchy over its sub-entities, i.e., a facet of a ${\textstyle N}$-sides polygon has hierarchy over the ${\textstyle N}$ edges that compose it. In turn each of the edges ${\textstyle {\boldsymbol {e}}^{a}}$ has hierarchy over its two vertices ${\textstyle {\boldsymbol {\mathrm {v} }}^{a},{\boldsymbol {\mathrm {v} }}^{a+1}}$. Figure 33a outlines the Contact Hierarchy for a triangle. The algorithm is organized as a sequence of three entity-checking levels. If a particle is in contact with the facet of a FE the contact search over its edges and vertices, which are in a lower hierarchy level, is discarded (see fig. 33). Otherwise, if contact with the FE facet does not exist, the contact check should continue over the sub-entities. Similarly, at the edges level, any contact with an edge cancels out further contact checks for those two vertices belonging to that edge. It does not cancel out, however, the contact check with the other edges because they are at the same hierarchy level. Table 7 in section 3.4.1 displays the pseudocode of the contact Type detection.

 (a) Contact Type Hierarchy for a triangle (b) Contact with facet. Edges and vertices are discarded from contact check Figure 33

Every time a new contact entity is determined by the Contact Type Hierarchy, the Distance Hierarchy (section 3.4.2) takes place immediately after. The Distance Hierarchy will determine if the new contacting entity found is redundant or non-valid, if it cancels out the previously found ones or if it is a new valid contacting entity to be considered for the DE.

For any valid contact entity the geometrical contact characteristics that will be stored are:

• The contact Point ${\textstyle {\boldsymbol {Pc}}}$.
• The FE nodal weights.
• The contact type: Facet, Edge or Vertex.

Note that some of the geometrical characteristics such as the distance, the normal vector or the contact local axis can be recalculated later when the contact constitutive law is applied and, thus, it is optional to store them here at this stage.

#### Facet level

The check proceeds in the same way as explained in section 3.3, checking for the intersection of the DE with the plane formed by the FE (section 3.3.1). If the Fast Intersection Test has been performed previously ${\textstyle \left|d_{\pi }\right|\leq R}$ is necessarily true since contact has been found for this FE. Otherwise, if no previous Fast Intersection Test has been carried out, this condition applies now to discard FE without contact.

Next, the Inside-Outside test (section 3.3.2) has to be performed. This test will tell us whether the projection ${\textstyle C_{\pi ^{m}}}$ (equation 103) lies on the facet (inside the FE) or it is outside, contacting with the edges or vertices. Fig. 34 shows two examples where the projection ${\textstyle C_{\pi ^{m}}}$ is inside and outside the FE facet.

 (a) ${\displaystyle C_{\pi ^{m}}}$ inside the facet (b) ${\displaystyle C_{\pi ^{m}}}$ outside the facet Figure 34: Example of projection ${\displaystyle C_{\pi ^{m}}}$ inside and outside the FE facet

The values of the cross product sign ${\textstyle s^{a}}$ obtained from equation 105 for every edge ${\textstyle {\boldsymbol {e}}^{a}}$ are used to obtain the weights of the shape function at the contact point. The areas needed for the calculation are simply one half of the cross product sign: ${\textstyle \Delta _{a}=s^{a}/2}$.

The weights of the nodal shape functions on the contact point are then calculated. For a triangle:

 ${\displaystyle N_{1}={\frac {\Delta _{2}}{{\hat {\Delta }}_{T}}},\quad N_{2}={\frac {\Delta _{3}}{{\hat {\Delta }}_{T}}},\quad N_{3}={\frac {\Delta _{1}}{{\hat {\Delta }}_{T}}}\qquad {\textrm {where}}\qquad {\hat {\Delta }}_{T}=\Delta _{1}+\Delta _{2}+\Delta _{3}}$
(112)

For 4-nodded convex quadrilaterals (fig. 35 the following expression can be applied as introduced in Zhong [114]):

 Figure 35: Triangular areas for the calculation of shape function values in a planar convex quadrilateral

 ${\displaystyle N_{1}={\frac {\Delta _{2}\Delta _{3}}{{\hat {\Delta }}_{Q}}},\quad N_{2}={\frac {\Delta _{3}\Delta _{4}}{{\hat {\Delta }}_{Q}}},\quad N_{3}={\frac {\Delta _{4}\Delta _{1}}{{\hat {\Delta }}_{Q}}},\quad N_{4}={\frac {\Delta _{1}\Delta _{2}}{{\hat {\Delta }}_{Q}}}}$ ${\displaystyle {\textrm {where}}\qquad {\hat {\Delta }}_{Q}=(\Delta _{1}+\Delta _{3})(\Delta _{2}+\Delta _{4})}$
(113)

Note that if any of the cross product signs ${\textstyle s^{a}}$ evaluated with respect to the edge ${\textstyle {\boldsymbol {e}}^{a}}$ yields a negative value the check stops since the projection of the centre ${\textstyle C_{\pi ^{m}}}$ lies outside. The current edge index ${\textstyle {index}_{e}}$ is stored and it will be the first to be checked as it has been appointed in section 3.3.3.

If the projection ${\textstyle C_{\pi ^{m}}}$ (equation 103) lies inside the facet, it becomes the contact point ${\textstyle {\boldsymbol {Pc}}}$. Due to the highest hierarchy level of the facet, the Contact Type Hierarchy finishes here for this FE. The Distance Hierarchy is now called and all the necessary contact characteristics are saved.

#### Edge level

Here the edge check (section 3.3.3) has to be applied for every edge ${\textstyle {\boldsymbol {e}}^{a}}$ with ${\textstyle a\in \left[{index}_{e},N\right]}$ in a ${\textstyle N}$-sided FE starting with the first edge that yielded an outside status at the Facet level.

When contact with the edge ${\textstyle {\boldsymbol {e}}^{a}}$ is found the check at the lower level for the vertices associated to it, ${\textstyle {\boldsymbol {\mathrm {v} }}^{a}}$ and ${\textstyle {\boldsymbol {\mathrm {v} }}^{a+1}}$, is discarded (fig. 17.1.2). The contact check with the following edges can not be discarded, however, since they are at the same hierarchy level in terms of Contact Type. The Distance Hierarchy will determine the validity of the new contact and eliminate or substitute previous ones. This is a key difference with the Fast Intersection Test where the check automatically stops once a contact entity is found.

 (a) Contact with edge. Vertices belonging to that edge are discarded (b) Weights for an edge contact in a triangle

The nodal weights can be obtained from the ${\textstyle \eta }$ parameter (equation 110) at the edge ${\textstyle {\boldsymbol {e}}^{a}}$. Fig. 17.1.2 shows graphically how ${\textstyle \eta }$ is determined,

 ${\displaystyle N_{a}=1-\eta ,\quad N_{a+1}=\eta \quad (N_{N}=N_{0})}$
(114)

Equation 114 gives the values at the nodes connected to the edge ${\textstyle {\boldsymbol {e}}^{a}}$. The rest of nodes have a null value for its shape functions. If the edge contact check failed but the distance ${\textstyle d_{e}}$ (equation 106) is lower than the radius (${\textstyle d_{e}\leq R}$) the closest vertex (based on the calculation of ${\textstyle \eta }$) will be checked. The check will proceed in any case (found edge, found vertex or none) with the next edges.

#### Vertex level

The vertex check is described in section 3.3.4. Fig. 36 illustrates why the edge ${\textstyle {\boldsymbol {e}}^{a}}$ has hierarchy over its two vertices ${\textstyle {\boldsymbol {\mathrm {v} }}^{a},{\boldsymbol {\mathrm {v} }}^{a+1}}$ but not over the non-contiguous one ${\textstyle {\boldsymbol {\mathrm {v} }}^{a+2}}$. The shape function weights are ${\textstyle 1}$ for the found vertex and ${\textstyle 0}$ for the rest.

 Figure 36: Contact with edge and vertex. When contact exists with edge ${\displaystyle {\boldsymbol {e}}^{1}}$ it can also exist with vertex ${\displaystyle {\boldsymbol {\mathrm {v} }}^{3}}$

As usual the Distance Hierarchy is called after the contact is detected and, if the contact is valid, its characteristics are stored.

#### Contact Type Hierarchy scheme

The scheme of Table 7 assumes that the Fast Intersection Test has taken place already. For every DE the first loop is over the found neighbours. The check can be performed in parallel for every particle in the model.

 loop over every FE with contact neighbour ${\textstyle {\mbox{(e)}}^{m}}$.
 ${\textstyle (1)}$