## Conceptual Model

Discrete element method analyzes the interactions between a finite number of distinct blocks. Each block, in the implementation used in this work, is elas-tically deformable, while non-linear deformation is concentrated in the contacts.

All the equilibrium and constitutive equations have been obtained under the hypothesis of finite deformations. The deformation of a single block has been considered infinitesimal but the one of the whole system, due to possible roto-translation ofits parts, can be large compared to the usual finite elements problems.

Each block or body is composed by one or more 8-node hexahedrons depending on the complexity of the geometry. In the contact routine each body is identified and interact with the others with its external surface.

The external surface of a block is given by the external faces of the finite elements forming it (Fig. 2).

At the beginning of the analysis all the possible pairs of contacting faces are determined and stored in an array. The uncertainty is introduced by the roughness of this first check: there can be surfaces included in a contact pair which are not actually touching each other but are simply close enough to trick the contact criterion. During the calculation, at each time step, each contact pair is considered separately and is marked as active or not-active. A contact pair is considered active only if one surface is penetrating into the other and some other criteria established to avoid singular situations are satisfied. This procedure has to be repeated during the calculation every n time-steps to follow the physical evolution of the system.

Each active contact pair represents a contact between two blocks. The procedure to calculate the

element Mock combination boundary

Figure 2. 2D representation of blocks: (a) single finite element; (b) elements forming a block; (c) block composition; (d) block boundary.

forces belonging to each contact consists in the following:

- The larger surface is identified as the target-surface and the other as the hitting-surface.

- If the distance between the 4 nodes of the hitting-surface is higher than a tolerance, a regular grid of virtual nodes is created on it.

- The position of each node, real or virtual, of the hitting surface, relative to the target one, is calculated. If the node is not penetrating, the pair hitting-node/target-surface is considered non-active.

- If the node is penetrating the target surface, normal and shear forces are calculated and added to the forces acting on the two surfaces.

The normal force at contact is calculated using the penalty method. In this approach a small penetration is accepted during the impact of two blocks. After the penetration, a reaction force, proportional to the penetration and the penalty stiffness, orthogonal to the target surface, starts to counteract the superposition of the two bodies. Choosing an adequate penalty stiffness the penetration can be limited to an acceptable threshold.

The shear force has been modelled with elastoplas-tic behaviour. Once a hitting point penetrates for the first time in the target surface, the natural coordinates Xho of the normal projection of the point onto the surface are stored together with the point-surface contact pair. At the successive time-steps has to be checked if the point is still penetrating the surface. In this case the new coordinates of the hitting-point Xhi have to be calculated. The shear force has been modelled with elasto-plastic behaviour. At first the distance d between Xh0 and Xhi is calculated and multiplied

Figure 3. 2D representation of non-active and active contact pair: (a) lock and bounding box (it includes the block and permit a first rough contact test); (b) non-active contact pair (superposition of bounding-boxes, no block interpenetration); (c) active contact pair (superposition of bounding-boxes and block interpenetration).

Figure 3. 2D representation of non-active and active contact pair: (a) lock and bounding box (it includes the block and permit a first rough contact test); (b) non-active contact pair (superposition of bounding-boxes, no block interpenetration); (c) active contact pair (superposition of bounding-boxes and block interpenetration).

by the shear stiffness Kg. The value obtained is confronted with the absolute value of the normal penalty force Fn multiplied by the static friction coefficient

If dKs > • Fn, the shear force is still in the elastic range, its absolute value is given by d Kg, the direction and verse is given by the vector connecting Xhi and XH0.

If d• Kg < • Fn, the shear force, calculated as an elastic reaction to the displacement, exceeds the static friction force. In this case the direction and verse of the shear force are the ones calculated in the elastic range while the absolute value is scaled to • Fn.

The contact forces, together with the external ones, are used to calculate the displacements and the deformations of the blocks at the successive time-step. Central difference method has been implemented to solve the dynamic differential system given by