# Dissipative contacts and realistic block shapes for modeling rock avalanches

Vincent Richefeu, Guilhem Mollon, Dominique Daudon, Pascal Villard

► To cite this version:

Vincent Richefeu, Guilhem Mollon, Dominique Daudon, Pascal Villard. Dissipative contacts and realistic block shapes for modeling rock avalanches. Engineering Geology, 2012, 149-150, pp.78-92. 10.1016/j.enggeo.2012.07.021. hal-02002821

HAL Id: hal-02002821

https://hal.univ-grenoble-alpes.fr/hal-02002821

Submitted on 2 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dissipative contacts and realistic block shapes for modeling rock avalanches

Vincent Richefeu *, Guilhem Mollon, Dominique Daudon, Pascal Villard

UJF-Grenoble 1, Grenoble-INP, CNRS UMR 5521, 3SR Lab, Grenoble F-38041, France

A specific contact model was used in a code (discrete element method) to explain the mechanisms of energy dissipation by collision and friction during the propagation of a granular mass on a slope. The numerical model focused on both realistic block shapes and the relevance of the single-collision law that make physical parameters easily assessable. Identification of the contact parameters was carried out by means of the digital image analysis of two-body collisions. To this purpose, the dropping of single blocks on a flat surface was filmed from two angles with high speed cameras. The digital images acquired during the rebound were then analyzed to extract accurately the block kinematics (3D trajectory and velocities). The contact parameters were optimized by minimizing an error function obtained by comparison between the numerical predictions and the experimental results. Once the parameters were set, a simulation of the collective behavior of the release of piled and randomly poured bricks under the same conditions as those released experimentally as described in the literature, was carried out. The satisfactory match between the experiments and the numerical predictions showed that (i) the proposed collision laws are sufficient to describe with accuracy the energy dissipation that occurs during binary collision or during mass propagation, (ii) the optimization procedure enables correct identification of the parameters, and (iii) the initial layout of the blocks is of primary importance in this process.

## 1. Introduction

Understanding and predicting rock avalanches are key elements in risk management when developing mountainous areas. Due to the complexity of the mechanisms involved, the morphology of the deposit or the propagation distance of a rock mass is difficult to estimate accurately at the moment. On the one hand, the characteristics of the geological base (the initial fracture network which determines the shape and size of the blocks) or the topography of the avalanche area can be determined by on site observations. On the other hand, the phenomena of energy dissipation by friction, fragmentation, lamination, attrition or seismic wave generation are much more complex to estimate.

To introduce these dissipation mechanisms, the physical means by which energy is dissipated in normal and tangential directions is usually dealt with using a coefficient of restitution (CoR) and a coefficient of friction, respectively. The overall damping rate is first defined as a ratio between the velocity magnitude of a particle just after and just before impact. Typical values can range from 0.2 to 0.3 for soft contacts (e.g. impact on scree or loose soils Habib, 1976, 1977; Descoeudres, 1997; Pfeiffer and Bowen, 1989; Evans and Hungr, 1993), and from 0.7 to 0.8 for hard contacts. The CoR values vary considerably from one

site to another, even for similar slopes. This is due to the huge amount of parameters that are involved. In fact, besides the role played by the characteristics of the blocks (weight, size, shape) and those of the slope of the ground, an important part of the CoR is to be attributed to the kinematics of the impact (direction, speed and slope angle) (Wu, 1985; Chau et al., 1999; Wu et al., 2003; Heidenreich, 2004).

Another definition of the CoR makes use of a distinction between the normal, the tangential, and the rotational components of velocity (Wu, 1985; Descoeudres, 1997; Okura et al., 2000a). The values of the tangential CoR, typically ranging from 0.7 to 0.9, are greater than those of the normal CoR, which vary between 0.1 and 0.6 depending on the nature of the soil (Fornaro et al., 1990; Giani et al., 2004). The threshold for the rolling resistance in relation to the CoR ranges from 0.4 to 0.85 (Azzoni and De Freitas, 1995; Giani et al., 2004). In the field of trajectography, a CoR greater than 1 can be used according to some authors (see e.g. Evans and Hungr, 1993; Paronuzzi, 2009). Although this seems to be unphysical, a value larger than 1 is actually possible due to a different definition of the CoR.

Moreover, experimental studies of block impacts on slopes showed that increasing the impact velocity leads to a greater loss of energy and therefore a reduction in the CoR (Bozzolo et al., 1988; Urciuoli, 1988). This can be explained in some cases by the fact that for high velocities, a block penetrates further into the substrate which then hardens (Heidenreich, 2004). There exist many other possible explanations since the field of contact mechanics is very extensive. Keys for the

<footnote>* Corresponding author.
E-mail address: vincent.richefeu@hmg.inpg.fr (V. Richefeu).</footnote>
analysis of the phenomena involved in the impacts are provided in the literature (see e.g. Walton and Braun, 1986; Johnson, 1987; Tsuji et al., 1992b; Thornton et al., 2011).

Another parameter of analysis widely used is the CoR expressed in terms of energy rate, i.e., a value that expresses the energy after impact compared with that before impact (Bozzolo and Pamini, 1986; Azzoni and De Freitas, 1995; Chau et al., 1999). This definition using kinetic energy rates led to a ratio of squared velocities and thus equals the square of the CoR based on velocity. In the case of a vertical release, this coefficient can be seen as the potential energy of the block before and after impact and thus as the maximum height of the block (after impact) divided by its initial height.

The numerical models commonly used to simulate the propagation of granular materials, are continuum models based on fluid mechanics (Savage and Hutter, 1989; Hungr, 1995; Mangeney-Castelnau et al., 2003; Pirulli, 2009) or discrete element models (Cundall and Strack, 1979; Okura et al., 2000b; Taboada and Estrada, 2009). In an engineering context, rock flows and in particular their run out distances must be assessed by means of continuum models. What is meant here is more specifically shallow-layer models for which the granular mass is treated as a homogeneous pseudo-fluid (sometimes incompressible) using Saint-Venant equations. The rheological models used characterize specific dissipation phenomena located at the base of the flow or within the granular mass: frictional and Voellmy models (Voellmy, 1955; Hungr, 1995; Pirulli and Mangeney, 2008). The most sophisticated ones can simulate 3D flows (Denlinger and Iverson, 2004; Hungr and Evans, 2004; McDougall and Hungr, 2004) by considering the material carried along the path (Hungr and Evans, 2004; McDougall and Hungr, 2005; Pirulli, 2009; Sautier et al., 2010). The shallow-layer models are in general accurate if the dissipation parameters are retro analyzed. Unfortunately, the parameters obtained are not always satisfying from a physical point of view (if the suitable rheology is not explicitly taken into account).

Discrete models have the advantage in modeling the deterministic movements of a set of interacting particles. Many authors have used this model to study the propagation of granular materials (Campbell et al., 1995; Calvetti et al., 2000; Cleary and Prakash, 2004; Linares-Guerrero et al., 2007; Staron, 2008; Valentino et al., 2008; Favier et al., 2009; Taboada and Estrada, 2009). Taking into account realistic shapes is possible using polyhedra, spheres or rigid sets of spheres (clumps). The motion of each particle is governed by the fundamental principle of dynamics. Several types of dissipation models are proposed in the literature (e.g. viscous contact models Campbell et al., 1995; Cleary and Prakash, 2004) but few remain relevant to the mechanisms observed in collision. Some models (Cundall, 1987) assume that a rate of the driving force can be removed by means of a damping coefficient. This artifact of calculation, which has no real meaning, dissipates energy in an arbitrary manner. It affects both the kinematics of free bodies and bodies that interact with each other. While it hides much of the physical reality, this ratio has sometimes been used to study the propagation of granular materials (Calvetti et al., 2000; Tommasi et al., 2008; Valentino et al., 2008).

Discrete element approaches are not based on strong assumptions concerning the flow and "collective mechanisms" of dissipation. Unlike in the case of continuum models, the velocity profiles, the dissipation due to impacts and the friction in the bulk or bottom of the flow, the specific dissipations due to the terrain, etc. are not presupposed. According to (Salciarini et al., 2009), discrete models enable accurate modeling of the propagation phenomenon without reaching a level of discretization such that all existing mechanisms of interaction are translated. Discrete models are thus valuable tools that help in the derivation of the continuum model of dissipation.

The model proposed here is based (i) on realistic block shapes and (ii) on the definition of simple interaction laws that make physical parameters easily assessable. These laws incorporate the mechanisms of energy dissipation globally.

To validate this approach, which is extremely difficult to carry out in the case of rock avalanches, simulations were performed that replicate laboratory experiments conducted under idealized test conditions (Manzella and Labiouse, 2009). Model parameters, optimized by tests conducted on single brick release, were used to simulate the collective behavior of a set of bricks on a slope.

Before applying the model to realistic cases, this study aimed precisely to assess if the error induced by a very simple representation of each single impact has an influence on the overall behavior of the rock flow or if this error disappears owing to the high number of impacts. The numerical results presented in this paper tend to confirm that this assumption is acceptable.

One advantage of the numerical model is that it gives access to quantities which are difficult to assess experimentally at any point of the granular mass: the velocity and rotation of bodies, the energy dissipated by friction or collisions with the slope or within the mass movement, the nature of flow and the geometry of the final deposit.

## 2. Numerical model

A classical discrete element method to simulate rock avalanches was chosen. Although many references exist on the topic of smooth-DEM, most of them focus on quasistatic loadings and/or spherical bodies. In this section a numerical model dedicated to brick flow on rigid planes implemented within the C++ toolkit DEMbox (see Computational Granular Physics Gateway<fnref n="1" />) is presented. This section focuses on the implementation of particle shapes and contact laws, the block motion being detailed in Appendix A.

### 2.1. Particle shapes

In the rock avalanche problematic, the shape of the blocks is of primary importance and this must be taken into account explicitly in the model. Different strategies are possible (e.g. polyhedra, clumps) but spheropolyhedra were chosen. This choice has several advantages including highly simplified contact detection (Alonso-Marroquin, 2008). Fig. 1 shows an example of a spheropolyhedron. The shape is defined by a set of vertices interconnected by edges (lines) and faces (plane polygons corresponding to triangles in Fig. 1). The vertex positions of a body *i* are given by reference to the mass center $O_i$ in the body framework $R_i$ defined by the main directions of the inertial tensor. The body shape is then defined by sweeping a sphere of radius $r_i$ along each point of its edges and faces. From a mathematical viewpoint, these block shapes can be seen as the Minkowsky sum of a polyhedron and a sphere (van den Bergen, 2003). In practice, the contact position, the overlap and the local frame are determined by considering a few basic geometric computations based on the distances between points, lines and planes. This geometric trick allows the contact area between spheropolyhedra to be defined by a set of contact points resulting from elementary intersection tests involving the swept sphere radii: (1) vertex-vertex, (2) vertex-edge, (3) edge-edge, and (4) vertex-face. One can better appreciate the benefit of this method when considering for example the face-face intersection test: The latter is simply replaced by a set of edge-edge and vertex-face tests. The spheropolyhedra method has many other benefits such as the ability to define concave and/or hollow shapes. Also the contact normals are well-defined.

Since a spheropolyhedron is a rigid body, only the time evolution of the mass center position and overall rotation is computed (see Appendix A). The movement of the entities that compose the spheropolyhedron (namely the slave bodies) are governed by the relations of rigid motion.

<footnote n="1">¹ www.cgp-gateway.org</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Definition of a spheropolyhedron illustrated in the case of a tetrahedron.</figcaption></figure>

## 2.2. Contact force laws

The model proposed for a rock avalanche is based on body shapes which are more realistic than spheres as seen previously. For the contact force-laws, a simple formulation which nevertheless incorporated the energy dissipation due to block impacts was needed.

Considering the huge amount of uncertainties related to a natural event, it seemed totally impossible to predict the exact behavior of each particle in the flow. Thus, it was decided to focus on the energy loss linked to each impact, rather than to reproduce the exact physical phenomena related to this impact.

The energy loss may result from very complex physical mechanisms (heat production, wave propagation...) that are beyond understanding in the case of collective behavior. Moreover, the local mechanisms do not need to be precisely identified, especially since it will be necessary to identify the parameters involved. Minimalist laws were opted for, where only the rate of energy loss together with friction are required to dissipate the kinematic energy of the blocks.

In other words, we have therefore chosen to consider a coarser scale to take into account force transmission and dissipation mechanisms in a granular assembly that flows and then stops. It is important to stress here that this coarser scale is not unphysical but ignores some physical mechanisms that are involved at smaller scales.

It will be seen in the sequel that the simple laws proposed here are sufficient to satisfactorily describe the main rebound patterns. Obviously, the model here is not able to reproduce more complex behaviors such as rocking blocks (Bourgeot et al., 2006) or acoustic wave propagation in confined granular systems (Somfai et al., 2005).

Rock avalanches involve dynamic block movements. For this reason, damping models which affects the block movements with an artificial parachute cannot be used since it would lead to unphysical behavior. Another solution is to account for a local viscous damping at contact level. This solution was also rejected because, although it introduces a viscosity parameter which can be connected to a well defined dissipation rate in the case of single contact (Tsuji et al., 1992a), it is ill-defined in the particular case of the multiple contacts involved in the interactions of spheropolyhedra. More precisely, the effective mass $m_{eff}$ involved in the critical viscosity $2\sqrt{m_{eff}k_n}$ is not well defined for complex shapes and should depend on the positions of the contact points and their number.

The simplest formulation for the normal force $f_n$ makes use of a linear elastic law with two different stiffnesses in the case of loading or unloading (respectively $k_n^+$ and $k_n^-$) (Banton et al., 2009). When the overlap $h_n$ increases (i.e. $\Delta h_n \ge 0$), the normal force increment reads:

$$ \Delta f_n = k_n^+ \Delta h_n. \qquad (1) $$

Otherwise, if $\Delta h_n < 0$ and $h_n > 0$, the force $f_n$ is given by:

$$ f_n = k_n^- h_n. \qquad (2) $$

Fig. 2(a) is a plot of this force-law which illustrates the role of incremental loading.

It should be noted that the literature provides other contact laws that also introduce energy dissipation by means of a difference of normal stiffness for loading and unloading (see e.g. Luding et al., 2003; Oger et al., 1998). These laws have all the desired properties for rock flow simulations, and provide accurate results for the impacts. Unfortunately, this leads to numerical issues at the end of the flow where persistent contacts tend to "oscillate", and an additional damping (or smaller time step) is required.

On the contrary, the contact law shown in Fig. 2(a) has the ability to dissipate naturally the energy during these oscillations. It is thus much more convenient for the numerical modeling of rock flows under consideration here. It has the added advantage of using a single parameter for normal energy dissipation without the need to postulate on all phenomena actually involved (e.g. viscosity, plasticity). Since the focus is not on the actual value of the repulsive force, the "force jump" used in the proposed contact law is not an issue. This sounds unphysical from a contact mechanics point of view, but the apparent lack of realism in the definition of contact force laws is irrelevant for the problem treated here.

Now consider a normal impact, that is – in this case of smooth contact formulation – a contact loading up to a given overlap $h_n^m$, immediately followed by an unloading until $f_n = 0$. After the impact, one part $k_n^-(h_n^m)^2/2$ has been restored, and since the maximum energy that can be restored is $k_n^+(h_n^m)^2/2$, the dissipation rate in the normal direction is:

$$ e_n^2 = k_n^- / k_n^+. \qquad (3) $$

In this definition, the square indicates that the dissipation rate is expressed in terms of energy rate.

The friction forces $f_t \vec{t}$ in a local frame – the tangent direction $\vec{t}$ being opposite to the sliding direction – is incrementally updated as a function of the increment of relative displacement $\Delta h_t$ in the sliding direction (during a time step $\Delta t$). The relative rotation of the contact normal is accounted for by the strategy developed in (Hart et al., 1988). However for the sake of simplicity, the tangential law is

<figure><img src="image_3.png" /><figcaption>Fig. 2. Force laws: (a) Normal contact force as a function of the overlap, (b) tangent friction force as a function of the sliding displacement.</figcaption></figure>
presented here as in the 2D case (i.e. without taking into account the relative rotations due to twisting):

$$f_t = \min \left\{ \sum_{t_0}^{t} k_t \Delta h_t ; \mu f_n \right\} \quad (4)$$

where $k_t$ is the tangential elastic stiffness, $\mu$ is the coefficient of friction, and $t_0$ is the time at contact onset. This common force-law is plotted in Fig. 2(b).

### 3. Estimation of the model parameters

#### 3.1. Single-brick experiments

The geometrical characteristics of the small scale experiment (Manzella and Labiouse, 2009) to be reproduced are well known (see Section 1), but the physical parameters of the interaction model cannot be determined directly from the experimental results reported in this paper. Additional experiments are therefore performed using the same materials (small bricks and plastic supports called “forex”) as the ones used in (Manzella and Labiouse, 2009). Each of these experiments consists in observing the impact of a single brick on a support after a fall, either from a specifically designed launcher (allowing reproducible launches) or by launching it manually in various positions. The fall, the impact, and the rebound are filmed at 1000 images per second by two high-speed cameras positioned along orthogonal directions. Both brick-support and brick-brick impacts are tested. In the first case the impact surface is identical to the plastic support used by (Manzella and Labiouse, 2009), and in the second case it is composed of several bricks aligned and stuck onto a support. The experimental device is described in Fig. 3.

The fall and the rebound of the launched bricks are filmed at an average time interval of 200 ms (roughly 100 ms before impact and 100 ms after), i.e. 200 shots per camera. Fig. 4 provides an example of some of these shots for each camera with a 20 ms time step. The point of impact O of the brick is used as the origin of the spatial system of axis (x, y, z), and the instant of impact is used as the origin of the time axis. y is the vertical axis (pointing upwards), while x is the axis of camera 1 (pointing from the camera to point O) and z is the axis of camera 2 (pointing from O to the camera). It is therefore assumed that the two cameras have perfectly orthogonal horizontal axes, meeting in O. As shown in Fig. 4, four points (located at visible corners of the bricks and called A₁, A₂, A₃ and A₄) are followed during the image analysis. Points A₁ and A₂ are visible from both cameras, while A₃ and A₄ only appear on cameras 2 and 1 respectively. For each of the 200 shots of each camera, the positions of the points on the images (expressed in pixels) are estimated, either manually or using a digital image correlation technique. This operation provides

the trajectories of A₁, A₂, and A₄ (respectively A₁, A₂, and A₃) projected on camera 1 (respectively on camera 2), with a time step of 1 ms.

#### 3.2. Synchronization and scaling

The synchronization of the images consists in determining the exact time coordinate of each one. The cameras are not synchronous, but separated by a constant unknown time interval, and a common time of origin is chosen (moment of the impact). The position of the center of mass of the brick is roughly estimated as the middle of segment A₂A₄ (respectively segment A₂A₃) on camera 1 (respectively on camera 2). Provided the free-flight movement of the center of mass is parabolic (i.e. the velocities involved are low enough to disregard air friction), this movement is estimated for each camera from the points measured. The intersection of the parabolas before and after impact provides a satisfying estimate of the moment of impact, and therefore of the time of origin (Fig. 5).

Estimating the positions of points A₁ to A₄ in the physical system of axis (O, x, y, z) requires knowledge of the scale (expressed in meters/pixel) on each camera. This scale is different from one launch to another, depending on the position of the brick. It is assessed by measuring (in pixels) the apparent length of an edge of the brick the real length of which (in meters) is known. If this edge belongs to a normal plan to the camera axis (i.e. if one may consider that both its extremities are the same distance from the camera), the scale is obtained directly by dividing the actual length by the apparent length on the image. For a given impact, this scaling is performed for each camera, and the positions of points A₁ to A₄ with respect to O are evaluated. These coordinates are linked to the time coordinates of each image (estimated by the synchronization method presented earlier). This step provides measurements of x₁(t), y₁(t), z₁(t), x₂(t), y₂(t), z₂(t), x₃(t), y₃(t), y₄(t), and z₄(t). The redundancy on the coordinates y₁(t) and y₂(t) (provided by both cameras) makes manual corrections of the image scales possible in the case of poor correspondence.

#### 3.3. Determination of the parameters of the contact law

The experimental kinematics of the bricks before and after impact can be assessed from the trajectories of points A₁ to A₄ (see Appendix B). The brick (31 × 15 × 8 mm) is represented using a spheropolyhedron (the Minkowski sum of a parallelepiped and a small sphere with a 1 mm radius) provided in Fig. 6, and composed of 8 spheres (corners), 12 cylinders (edges), and 6 planes (faces). The determination of the four parameters ($e_n^2$, $\mu$, $k_n$, and $k_t/k_n$) of the contact law is performed by retro analysis on a simulation of the impact using a discrete element model which implements the contact law described earlier. The experimental kinematics of the brick being defined, the parameters of the fall (initial position,

<figure><img src="image_4.png" /><figcaption>Fig. 3. Experimental device used for identification of the parameters.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Example of images acquired from the orthogonal cameras: (a) camera 1, (b) camera 2.</figcaption></figure>

rotation, velocity, and angular velocity of the brick) are used as initial conditions in this simulation. For a given set of the four parameters of the contact law, the model provides a trajectory of the brick

after the impact, which may be compared to the experimental trajectory. The parameters of the contact law will be correct if the correspondence between the experimental rebound and the one

<figure><img src="image_6.png" /><figcaption>Fig. 5. Example of shot synchronization (the impact time is chosen as the time of origin).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Shape of a brick modeled on a spheropolyhedron.</figcaption></figure>
simulated is optimal. This correspondence is obtained by minimizing the error function $E_2$:

$$E_2 = \sum_{i=1}^{8} \sum_{t=0}^{\tau} d_i(t). \qquad (5)$$

In this expression, $i$ is an index defining one of the 8 corners of the brick, and $d_i(t)$ represents at a given time $t$ the distance between the positions of this point obtained from both the identified motion and from the simulated one. The error function is therefore a sum over time of the distances between the positions of the 8 corners obtained experimentally and numerically.

The choice of 8 points induces some redundancy but does not complicate the optimization process since it only modifies the computation of the error function. When dealing with retro analysis, redundancies are always welcome because they lead to more discriminating error functions and therefore to better estimates of the unknown parameters. The sum in Eq. (5) is computed from the moment of impact and during a given time $\tau$ which should be neither too short (the experimental and numerical trajectories would not have time to diverge and the error would have no meaning) nor too long (the error would only represent the correspondence between the trajectories of the center of mass, the influence of rotation becoming imperceptible). Taking this into consideration a reference time $\tau_{ref}$ can be defined as the ratio between the length of the brick and its characteristic velocity. For the tests here, $\tau_{ref}$ is about 0.01 s, and the normalized time is thus $\tau/\tau_{ref} = 0.1/0.01 = 10$. In other words, the time considered in the error function $E_2$ corresponds to a travel length 10 times longer than the length of the brick.

It should be noted that the minimization process only relies on kinematical measurements before and after impact, and therefore totally ignores any measurements of force. This is consistent with the fact that the chosen contact law is only meant to reproduce the kinematics of the impacting objects, considering the contact as a “physical black-box” (even though this contact has a duration in the numerical process). To that extent, it may be compared with the resolution of contacts performed in the contact dynamics approach (see e.g. Moreau, 1994; Jean, 1999; Radjai and Richefeu, 2009).

For a given impact, the minimization of $E_2$ leads to an optimum on the parameters of the contact law, which therefore provides the best numerical representation of the experimental impact. However, it was observed that, from one impact to another, the sets of optimal parameters were not strictly similar. Several reasons for this may be considered:

* The modeling of a contact involving numerous complex physical phenomena with only four parameters makes it impossible to account for certain influences, such as those of impact velocity or angle of attack.

* It is not possible to simulate accurately the exact geometry of the brick, especially the abrasion of the corners which leads to a contact surface that differs from the one introduced in the model.

* The methodology of identification of the experimental trajectory of the brick suffers from uncertainties, such as those stemming from the assumption of a constant scale for a given impact despite the relative motion between the brick and the camera.

Despite these observations, it is interesting to limit the number of contact parameters, and the proposed model still appears relevant. Moreover, the objective of this model is to simulate the collective motion of a rock avalanche, not to simulate accurately the motion of a single particle for which perfect repeatability cannot be attained. Indeed, a very small change in the initial conditions (position, rotation, velocities, geometry, sharpness of the corners) has a dramatic effect on the brick trajectory after the first impact. Actually, such repeatability is not mandatory when dealing with natural phenomena such as rock flows where the conditions of impacts (contact

geometries, material properties, etc.) are extremely diverse even for a single event. It is very difficult to model accurately all these conditions with only a few parameters. A better approach would be to find the parameters that lead to the best representation of the overall behavior of the flow, ignoring the approximations related to the representation of single impacts. It was therefore decided to diversify the impact conditions during calibration instead of repeating the same conditions several times. Given this choice, the optimal parameters will not be perfect for each single impact, but will be as good as possible for the wide variety of impact conditions. This choice cancels out the need to evaluate the respective parts of error related to experimental measurements or changes in the conditions of impact. However, in the perspective of applying this method to natural events, the present authors believe that the uncertainty linked to measurement is significantly lower than that linked to the impact conditions.

It may therefore be estimated that global parameters obtained from a number of single impacts should provide a representative behavior of the particles, despite a small inaccuracy. The global error function for a given type of contact is the sum of the error functions computed for the different experiments performed for this contact type (four launches for a brick-support contact, two launches for a brick-brick contact), and is given by:

$$E_{2G} = \sum E_2. \qquad (6)$$

The minimization of this global error leads to a global optimal set of parameters for the representation of a given type of impact. A representation of this error is provided in Fig. 7 for the brick-support contact, with respect to the four parameters $e_n^2, \mu, k_n$, and $k_t/k_n$. The error is plotted using curves of equal values of $E_{2G}$ in two planes: plane ($e_n^2, \mu$) with $k_n = 10^5$ and $k_t/k_n = 0.42$, and plane ($k_n, k_t/k_n$) with $e_n^2 = 0.53$ and $\mu = 0.46$.

Concerning the two parameters related to energy dissipation ($e_n^2, \mu$), Fig. 7a shows that a clear optimum appears, under the shape of a “well” of error around the point of coordinates $e_n^2 = 0.53$ and $\mu = 0.46$. These two values are therefore probably able to represent properly the average behavior of a brick-support contact. For the two stiffness parameters, the level curves of Fig. 7b show a less accentuated optimum. Note that the gray-level gradient of Fig. 7b is highly dilated when compared with Fig. 7b to ease the location of the minimum value; the surface is actually almost flat. The parameter of normal stiffness $k_n$ does not seem excessively critical between $10^4$ and $10^6$, whereas the parameter of tangential stiffness $k_t/k_n$ has a clearer optimum around 0.42. Similar response surfaces are plotted for the brick-brick contact. The optimum parameters obtained for the two contact types are provided in Table 1. It clearly appears that a contact between two bricks is much more dissipative than a contact between a brick and the plastic support.

It should be noted that the performed calibration is probably shape-dependent and also depends on the locus of impact. For our brick shape, the obtained parameters for each launch are quite similar and the shape-dependency is thus limited. This is surely not the case for shapes with a high aspect ratio.

Fig. 8 presents a summary of the results obtained for an example of brick-support impact (corresponding to Fig. 4). For each camera and for a time step of 20 ms, it provides successively the positions of the brick captured by the camera (left), the trajectories estimated by minimization of the error function $E_1$ (middle), and the rebound provided by the numerical simulation (right) for the set of parameters provided by minimization of $E_{2G}$ and given in Table 1. It clearly appears that the globally optimal parameters (determined after several impacts) are able to simulate, with satisfactory accuracy, a single impact.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Contour lines of the global error function $E_{2G}$ (brick/support contact): (a) {$k_n$; $k_t/k_n$} parameter-space, (b) {$e_n^2$; $\mu$} parameter-space.</figcaption></figure>

# 4. Validation of the model

## 4.1. Description of the experiment chosen for validation

In order to assess the ability of the proposed model to predict the behavior of a rock avalanche, a validation is undertaken by attempting to reproduce some experimental results from (Manzella and Labiouse, 2009). A series of launches of a large number of small bricks were carried out on a device composed of two rectangular boards (3 × 4 m) of forest (a sort of plastic), linked by a hinge. The first board was fixed

<table><caption>**Table 1**
Optimal parameters obtained by minimizing $E_{2G}$.</caption><thead><tr><th></th><th>e<sub>n</sub><sup>2</sup></th><th>μ</th><th>k<sub>n</sub></th><th>k<sub>t</sub>/k<sub>n</sub></th></tr></thead><tbody><tr><td>Brick/support</td><td>0.53</td><td>0.46</td><td>10<sup>5</sup></td><td>0.42</td></tr><tr><td>Brick/brick</td><td>0.13</td><td>0.86</td><td>10<sup>5</sup></td><td>0.27</td></tr></tbody></table>

and horizontal, while the second one was inclined at a user-defined angle. A rectangular box (height 0.2 m, width 0.4 m, depth 0.6 m) was filled with a given amount of material, and positioned at a determined height on the inclined plane. A trap was open to release the tested material. This material traveled along the slope until it landed on the horizontal plane. A high-speed camera was used during the experience in order to evaluate the position and velocity of the mass-front. The morphology and the dimensions of the final deposit of material were determined by a fringe projection technique. Several parameters were studied: the nature of the material released (gravel, small bricks), the slope angle, the fall height, the volume of material, and the type of support. The validation proposed in the present article focuses on two experiments carried out on small bricks, of average dimensions 31 × 15 × 8 mm, with a density evaluated at 1700 kg.m⁻³. For both experiments, the slope angle was fixed at 45°, and the launch height was equal to 1 m (Fig. 9). During the first test, a volume of 40 l of bricks was poured randomly into the launching box. This represents
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Positions of the brick each 20 ms: (left) camera shots; (middle) identification by minimization of $E_1$; and (right) numerical model with an optimal set of parameters obtained by minimization of $E_{2G}$. (a) Camera 1, (b) camera 2.</figcaption></figure>

an apparent density of 1000 kg.m⁻³, while during the second test they were piled orderly to reach also a volume of 40 l (apparent density of 1600 kg.m⁻³). It should be noted that the materials used in Section 4

for the determination of the contact parameters (small bricks and
forex support) were identical to the ones used during the experiments
by (Manzella and Labiouse, 2009) (Fig. 10).

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Proportions of the reference experimental device and definition of measured parameters.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Materials used by (Manzella and Labiouse, 2009): (a) bricks randomly poured and (b) bricks piled orderly.</figcaption></figure>

## 4.2. Simulation of randomly poured bricks

In a first attempt to validate the proposed model, the experiment involving randomly poured bricks is reproduced numerically. The initial state of the experiment is reproduced by launching simultaneously 6300 virtual bricks (Fig. 6) in a fictive box, each brick being oriented in a random direction. In order to reach a volume of material of 40 l (corresponding to an apparent density of 1000 kg.m⁻³), several trials are needed with various values of the contact parameters. When the target density is achieved, the contact parameters are restored to their optimal values (Table 1), and the fictive box is installed in conformity with the geometry of the experimental device. The fall of the bricks starts when the lower face of the box is deleted.

Fig. 11a and b shows several stages of the avalanche, respectively a vertical view and a perspective view. In Fig. 12, a comparison between the experimental and numerical results is provided in terms of the contour of the deposit (in plane and in elevation), and of the mass-front velocity on the horizontal plane. A very satisfying correspondence appears between the numerical and experimental deposits projected in a horizontal plane, but this correspondence is not as accurate for the projection in a vertical plane, the numerical deposit being much higher than the experimental one. The curve of the mass-front velocity with respect to the position of the mass-front shows a correct qualitative correspondence. The observed quantitative differences may be linked to a possible inaccuracy of the model, but are more probably related to a lack of precision in the definition of the motion of the mass-front. Indeed, the position of the mass-front has no rigorous and objective definition. In the experimental tests it was defined from an automatic algorithm on each shot of the camera (using a criterion of contrast between the released material and the plastic support), while in the numerical simulations it is based on a criterion of spatial density of the bricks (not detailed here). This difference of definition might be sufficient to explain the observed differences. Moreover, the authors of the

experiments point out a possible inaccuracy in their results due to difficulties related to low-contrast images (Manzella and Labiouse, 2009). Taking this into consideration, the numerical results seem relevant. Qualitatively, Fig. 12c provides some interesting observations about the behavior of the avalanche. One may first observe a substantial deceleration on impact of the mass-front on the horizontal plane (between 0 m and 0.2 m), followed by a zone in which the velocity decreases a little, corresponding to the accumulation of the material on the plane (between 0.2 m and 0.6 m), and finally a further substantial deceleration until the end of the motion (between 0.6 m and 0.8 m).

Fig. 13 shows a comparison between quantitative results obtained experimentally and numerically. The dimensions of the deposit (length L, runout R, width W, and height H) are very satisfactorily reproduced by the model, except for that of the height which is overestimated by 60%. The position $X_{CM}$ of the center of the mass of the deposit is correctly assessed, and both angles of propagation (travel angle related to the center of mass $\phi_{CM}$ and fahrböschung $\phi_{app}$ related to the extreme points of the deposit) are very accurately simulated. These results emphasize the ability of the model to simulate all the energy dissipation phenomena inside the mass of material, despite the fact that the parameters of contact law were derived from quite simple tests using a single brick motion.

## 4.3. Simulation of orderly piled bricks

A second experiment from (Manzella and Labiouse, 2009) is simulated. This experiment deals with 40 l of the same bricks which are positioned and piled along regular alignments. This stack is reproduced numerically and leads to a much larger number of bricks than in the first simulation (roughly 10,000 instead of 6300), due to the greater density of the initial batch of blocks for an apparently same volume. The numerical and experimental results are provided in Fig. 14, in terms of deposit morphology and of mass-front velocity. The horizontal and vertical contours of the numerical deposit show a poorer correspondence than in the case of randomly positioned bricks. Indeed, it appears that the numerical model leads to a higher, longer, and narrower deposit than the one observed experimentally, even if the global position is correct. The numerical curve of the mass-front velocity shows a similar behavior as that in the first simulation, but with greater velocities. The final deceleration and stopping of the mass are correctly reproduced, but the experimental curve does not show any brutal deceleration during the first impact of the avalanche on the horizontal plane as predicted by the numerical curve.

The differences observed of the morphology of the deposit may be explained by poor reproduction of the initial state of the batch of particles. Indeed, the numerical alignments and shapes of the bricks are absolutely perfect (in terms of regularity and parallelism), which is not the case of the experimental bricks (Fig. 10). This experimental inaccuracy is certainly understandable owing to the small differences in shapes and sizes of the bricks and the difficulty in achieving a perfect alignment of such a large number of elements. In such conditions, perfect accuracy in the numerical reproduction of the experimental initial state appears unattainable. However, one may presume that introducing a small level of perturbation in the initial alignment of the particles would improve the numerical results. To validate this assumption, the shape of the numerical bricks is slightly modified by the introduction of a bevel of 2° on three faces. The modified bricks therefore are not longer parallelepipeds, although their global dimensions are very similar to those of the original bricks. The modified bricks are stacked in the virtual box following the same alignment but with randomly varying orientations, so that the bevel introduced will not trigger any asymmetry in the flow. The numerical results are compared with the experimental ones in Fig. 15.

It appears that the modification of the shape of the bricks improves the correspondence between the numerical and experimental results.
<figure><img src="image_12.png" /><figcaption>Fig. 11. Numerical results for an avalanche of randomly poured bricks. (a) Top view and (b) perspective view.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Comparison of experimental and numerical results for the release of initially disordered bricks. (a) Horizontal and (b) longitudinal cross-sections of the deposits; (c) mass-front velocities.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Quantitative comparison between experimental and numerical results for initially disordered bricks: length *L*, runout *R*, width *W*, position of the center of the mass $X_{CM}$, travel angle $\phi_{CM}$, fahrböschung $\phi_{app}$ and height *H*.</figcaption></figure>

The comparison between Figs. 14 and 15 shows that the bevel introduced reduces the length of the deposit and increases its width, allowing a better representation of the experimental deposit. This modification also leads to a small reduction of the mass-front velocity, although the global shape of the curve remains unchanged. The differences observed between the experimental and numerical velocities may be related to the same reasons as those in the case of randomly positioned bricks.

Fig. 16 provides a quantitative comparison between the deposits obtained experimentally and numerically (with perfect and modified bricks). It is clear that the introduction of the 2° bevel substantially improves the estimation of the dimensions of the deposit (reduction in the length, in the runout and in the height, increase in the width) and in the position $X_{CM}$ of the center of mass. The travel angle and fahrböchung, which are predicted quite well in the case of perfectly shaped bricks, are not modified in any way.

## 4.4. Discussion

These results enlighten the crucial impact of the initial state of the mass of material as it travels along a slope, and therefore the shape and dimensions of the final deposit, for an apparently similar shape

and volume of the initial mass (Fig. 17). The results of the simulation with randomly positioned bricks show that the proposed model is very capable of predicting the behavior of the flow of a group of particles, if the initial state of this group has no specific order or symmetry. It is likely that this ability is also satisfactory in the case of a specific initial state (i.e. with perfect alignment), but this assertion is more difficult to validate since such an initial state is almost impossible to reproduce experimentally. However, it was shown that the artificial introduction of a geometrical perturbation of the particle shape and alignment is able to reproduce satisfactorily a specific but not perfect initial state. Nevertheless such a perturbation is quite difficult to define *a priori*.

## 5. Conclusion

A discrete element numerical model based on a dissipative contact model and realistic block shapes has been applied to the propagation of granular materials down a slope. Simulations of releases with piled and randomly packed bricks were performed and compared with results from the literature. The parameters of the contact model, which all have physical meaning, were assessed on the basis of experiments of single brick bounce on a flat surface (or on another fixed brick) for

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Comparison of experimental and numerical results for a release of perfectly shaped bricks which were initially ordered. (a) Horizontal and (b) longitudinal cross-sections of the deposits; (c) mass-front velocities.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Comparison of experimental and numerical results for a release of modified bricks initially ordered. (a) Horizontal and (b) longitudinal cross-sections of the deposits; (c) mass-front velocities.</figcaption></figure>

several angles of incidence and different initial drop heights. An optimization procedure was used to determine the most appropriate parameter set able to satisfactorily reproduce all experimental tests that were conducted. The model's ability to report the behavior of particle flow was demonstrated by comparing the numerical predictions with the experimental releases (10000 piled bricks and 6300 randomly packed bricks). The match between numerical and experimental results for the release of randomly packed bricks proved to be very satisfactory both for the flow kinematics, the front velocity, the runout or the deposit morphology. The differences observed between simulations and experimental results in the case of piled bricks are for their part attributed to a poorly controlled arrangement of the stack in the initial test case. If a slight disorder in the organization of bricks (due to the introduction of small imperfections in the brick shapes) leads to a better consistency of results, a further comparison in this case does not seem necessary given the difficulty to be faced in perfectly replicating the original experimental packing.

For the tested materials (bricks and slightly deformable plane support) the contact laws and procedures for optimizing the parameters proposed are greatly suited to account for the rebound kinematics and mechanisms of energy dissipation whether the particles are

isolated or in assembly (dissipation by friction or collision). For applications to real events such as rock avalanches, the contact laws and the optimization procedure for parameter assessment must be tested and possibly improved to incorporate more complex interaction mechanisms such as excessive penetration of blocks on the sloping ground, the movement and tearing of the bedrock or rock fractures. The advantage of the discrete model is that it enables, provided that the physical mechanisms are controlled, the integration of specific dissipative laws that are features of flow kinematics and adapted to each type of event. Analysis and backward analysis of events that took place will certainly improve the model and provide predictive results on future events.

## Acknowledgments

This research project was performed as part of the European project ALCOTRA-MASSA, with financial support from the European Funds for Regional Development (FEDER). Special thanks to V. Labiouse from the Rock Mechanics Laboratory (EPFL, Lausanne, Switzerland) for providing the experimental material that was used to assess the required numerical parameters.

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Quantitative comparison between experimental and numerical results for perfect and modified bricks initially ordered: length *L*, runout *R*, width *W*, position of the mass center *X*<sub>CM</sub>, travel angle *φ*<sub>CM</sub>, fahrböschung *φ*<sub>app</sub> and height *H*.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 17. Deposit morphology obtained by the numerical model. The initial packing had an apparent volume of 40 l and was (a) randomly poured or (b) perfectly piled.</figcaption></figure>

## Appendix A. Block motions

The algorithm for the classical discrete element method involves two stages for each rigid body i: (1) the computation of the resultant forces $F_i$ and moment $M_i$ from volume and contact forces (see Section 2.2); and (2) the time integration of Newton's second law (for translations) and Euler's equations (for rotations). This movement integration is performed by means of the velocity-Verlet scheme (Allen and Tildesley, 1989) which is a satisfactory compromise between the accuracy of the block velocities (for both translations and rotations) and memory saving.

Each body i is kinematically defined in the global framework $\mathcal{R}$ by its mass-center position $\vec{r}_i$, its mass-center velocity $\vec{v}_i$, its angular position given by a unit quaternion $\hat{q}_i$, and its angular velocity $\vec{\omega}_i$. Unit quaternions provide a convenient mathematical notation for representing orientations and rotations of objects in three dimensions. Compared to Euler angles they are simpler to compose and avoid the problem of gimbal lock. Compared to rotation matrices they are more numerically stable and use less memory. Somehow, unit quaternions can be interpreted as rotation matrices but they only hold one scalar value and one vector: $\hat{q}_i = [s_i, \vec{v}_i]$. There exists a largely developed mathematical background concerning quaternions, but all this knowledge is not absolutely necessary when dealing with rotations.

For translation motion, the velocity-Verlet algorithm can be summarized as follows for each time step $\Delta t$:

$$ \left\{ \begin{array}{l} \vec{r}_i(t + \Delta t) = \vec{r}_i(t) + \vec{v}_i(t)\Delta t + \frac{1}{2}\vec{a}_i(t)\Delta t^2 \\ \vec{v}_i(t + \Delta t) = \vec{v}_i(t) + \frac{1}{2}[\vec{a}_i(t) + \vec{a}_i(t + \Delta t)]\Delta t, \end{array} \right. \qquad (7) $$

with

$$ \vec{a}_i = \frac{F_i}{m_i} + \vec{g}, \qquad (8) $$

where $m_i$ is the mass of the body i, and $\vec{g}$ is the gravity acceleration.

For rotation motion, the velocity-Verlet algorithm is also used to determine the angular positions and velocities of the bodies. It formally reads:

$$ \left\{ \begin{array}{l} \hat{q}_i(t + \Delta t) = \hat{q}_i(t) + \hat{q}_i(t)\Delta t + \frac{1}{2}\hat{\tilde{q}}_i(t)\Delta t^2 \\ \hat{\omega}_i(t + \Delta t) = \hat{\omega}_i(t) + \frac{1}{2}[\hat{\vec{\omega}}_i(t) + \hat{\vec{\omega}}_i(t + \Delta t)]\Delta t. \end{array} \right. \qquad (9) $$

In this scheme, the first and second time derivative of the quaternion can be expressed in terms of the angular velocity vector as:

$$ \dot{\hat{q}}_i(t) = \frac{1}{2} \hat{\vec{\omega}}_i(t) * \hat{q}_i(t), \qquad (10) $$

and

$$ \ddot{\hat{q}}_i(t) = \frac{1}{2} \hat{\vec{\omega}}_i(t) * \hat{q}_i(t) + \frac{1}{4} \hat{\vec{\omega}}_i(t) * \hat{\vec{\omega}}_i(t) * \hat{q}_i(t), \qquad (11) $$

where $\hat{\vec{x}}$ denotes the quaternion $[0, \vec{x}]$, and operator * is the Hamilton product defined by

$$ [s_1, \vec{v}_1] * [s_2, \vec{v}_2] = [s_1 s_2 - \vec{v}_1 \vec{v}_2, s_1 \vec{v}_2 + s_2 \vec{v}_1 + \vec{v}_1 \times \vec{v}_2]. \qquad (12) $$

In Eq. (9), the derivative of angular velocities of each of the bodies i is obtained from Euler's equations as follows (subscripts i are removed to facilitate reading):

$$ \left\{ \begin{aligned} \dot{\omega}_1^* &= \{M_1^* - (I_2^* - I_2^*)\omega_2^*\omega_3^*\}/I_1^* \\ \dot{\omega}_2^* &= \{M_2^* - (I_1^* - I_3^*)\omega_3^*\omega_1^*\}/I_2^* \\ \dot{\omega}_3^* &= \{M_3^* - (I_2^* - I_1^*)\omega_1^*\omega_2^*\}/I_3^* \end{aligned} \right. \qquad (13) $$

where components 1, 2 and 3 are expressed in the body frame (superscripts * are added here for memory).

It is important to notice that, in Eqs. (8) and (13), $\vec{F}_i$ and $\vec{M}_i$ are obtained from the contact force laws which depend on body positions ($\vec{r}$ and $\hat{q}$) at time $t$, and their mean velocities between times $t - \Delta t$ and $t$.

## Appendix B. Identification of experimental trajectories

From the trajectories of points $A_1$ to $A_4$ obtained after synchronization and scaling, it is possible to define the movement of the whole brick after a retro-analysis study. First, one has to put into equation the position of the four points using a reduced number of parameters related to the brick motion. These parameters are then adjusted by minimizing an error function in order to obtain the best correspondence between measured and calculated positions of the four points. Impacts are obviously not instantaneous, but the time lag between two images is so much greater that their duration can be ignored.

At any time $t$, the coordinates $\vec{r}_i^c(t) = (x_i^c(t), y_i^c(t), z_i^c(t))^T$ of a point $A_i$ in the global axis system are obtained from the trajectory
<figure><img src="image_19.png" /><figcaption>Fig. 18. Coordinates of measured points A₁ to A₄ as a function of time. The optimum trajectory is obtained by minimization of the function E₁.</figcaption></figure>

of the center of mass G of the brick and from its rotation around G:

$$ \vec{r}_i^c(t) = \vec{r}_G(t) + \begin{bmatrix} a_{xx}(t) & a_{xy}(t) & a_{xz}(t) \\ a_{yx}(t) & a_{yy}(t) & a_{yz}(t) \\ a_{zx}(t) & a_{zy}(t) & a_{zz}(t) \end{bmatrix} \vec{r}_i^0. \quad (14) $$

In this expression, the coordinates of G at the moment t are given by:

$$ \vec{r}_G(t) = \left\{ \begin{array}{l} \vec{r}_G^0 + \vec{v}_G^- t - \vec{g} \quad t^2/2 \quad \text{if } t \le 0 \\ \vec{r}_G^+ + \vec{v}_G^+ t - \vec{g} \quad t^2/2 \quad \text{if } t > 0 \end{array} \right. \quad (15) $$

The matrix of rotation at a time t is given by:

$$ \left\{ \begin{aligned} a_{xx}(t) &= \cos\theta_y(t) \cos\theta_z(t) \\ a_{xy}(t) &= -\cos\theta_y(t) \sin\theta_z(t) \\ a_{xz}(t) &= \sin\theta_y(t) \\ a_{yx}(t) &= \cos\theta_x(t) \sin\theta_z(t) + \cos\theta_z(t) \sin\theta_x(t) \sin\theta_y(t) \\ a_{yy}(t) &= \cos\theta_x(t) \cos\theta_z(t) - \sin\theta_x(t) \sin\theta_y(t) \sin\theta_z(t) \\ a_{yz}(t) &= -\cos\theta_y(t) \sin\theta_x(t) \\ a_{zx}(t) &= \sin\theta_x(t) \sin\theta_z(t) - \cos\theta_x(t) \cos\theta_z(t) \sin\theta_y(t) \\ a_{zy}(t) &= \cos\theta_z(t) \sin\theta_x(t) + \cos\theta_x(t) \sin\theta_y(t) \sin\theta_z(t) \\ a_{zz}(t) &= \cos\theta_x(t) \cos\theta_y(t) \end{aligned} \right. \quad (16) $$

where

$$ \vec{\theta}(t) = \begin{cases} \vec{\theta}^0 + \vec{\omega}^- t & \text{if } t \le 0 \\ \vec{\theta}^0 + \vec{\omega}^+ t & \text{if } t > 0 \end{cases} \quad (17) $$

with

* $\vec{r}_i^0 = (x_i^0, y_i^0, z_i^0)^T$ Z, the coordinates of Aᵢ in the local frame of the brick.

* $\vec{g} = (0, 9.81, 0)^T$ Z, the gravity acceleration.

The motion of the brick before and after impact is therefore totally defined by the following 18 parameters:

* $\vec{r}_G^0 = (x_G^0, y_G^0, z_G^0)^T$ Z, the coordinates of G at t=0.

* $\vec{\theta}^0 = (\theta_x^0, \theta_y^0, \theta_z^0)^T$ Z, the rotation of the brick around G at t=0.

* $\vec{v}_G^- = (v_{Gx}^-, v_{Gy}^-, v_{Gz}^-)^T$ Z, the velocity vector at t=0⁻ (just before impact).

* $\vec{\omega}^- = (\omega_x^-, \omega_y^-, \omega_z^-)^T$ Z, the angular velocity vector at t=0⁻.

* $\vec{v}_G^+ = (v_{Gx}^+, v_{Gy}^+, v_{Gz}^+)^T$ Z, the velocity vector at t=0⁺ (just after impact).

* $\vec{\omega}^+ = (\omega_x^+, \omega_y^+, \omega_z^+)^T$ Z, the angular velocity vector at t=0⁺.

These parameters define the position and rotation of the brick at the moment of impact, and the velocity and angular velocity vectors immediately before and after impact. A given set of these parameters will provide a good description of the brick motion if the correspondence between the measured coordinates (e.g. x₁(t)) and the analytical coordinates (e.g. xᵢ(t)) is optimal. This correspondence is obtained by minimizing the following error function E₁:

$$ E_1 = \sum_{S_1} \sum_{p=1.2.4} \left\{ \begin{aligned} & [y_p(t) - y_p^c(t)]^2 + [z_p(t) - z_p^c(t)]^2 \\ & + \sum_{S_1} \sum_{p=1.2.4} \left\{ [x_p(t) - x_p^c(t)]^2 + [y_p(t) - y_p^c(t)]^2 \right\} \end{aligned} \right\} \quad (18) $$

where S₁ and S₂ refer to sets of images acquired by camera 1 and 2, respectively, over time. This minimization is performed with respect to the 18 parameters describing the motion, and with respect to coordinates $(x_1^0, y_1^0, z_1^0, x_2^0, y_2^0, z_2^0, x_3^0, y_3^0, z_3^0, x_4^0, y_4^0, z_4^0)$ of the points A₁ to A₄ in the local system of axes of the brick. Indeed, these coordinates are not precisely known because of possible uncertainties concerning the size of the brick, abrasion of the corners of the brick, and imprecise estimation of the points coordinates on images. The minimization is therefore performed on 28 variables. Including these coordinates in the minimization improves the correspondence between measured and analytical trajectories of the four points, and provides a better estimate of the brick motion. A comparison between the measured coordinates of the points and the coordinates computed with the optimal set of parameters (i.e. leading to a minimization of E₁) is proposed in Fig. 18. The correspondence is very satisfying for the 12 coordinates, and confirms that the minimization converged. This optimization procedure is used to determine with precision the brick motion for a total of six impacts (four brick-support impacts
and two brick-brick impacts), with different launch height, inclinations, and initial velocities, in order to generalize the results.

## References

Allen, M.P., Tildesley, D.J., 1989. Computer Simulation of Liquids. Clarendon Press, New York, NY, USA.

Alonso-Marroquin, F., 2008. Spheropolygons: a new method to simulate conservative and dissipative interactions between 2D complex-shaped rigid bodies. EPL (Europhysics Letters) 83 (1).

Azzoni, A., De Freitas, M., 1995. Experimentally gained parameters, decisive for rock fall analysis. Rock Mechanics and Rock Engineering 28 (2).

Banton, J., Villard, P., Jongmans, D., Scavia, C., 2009. Two-dimensional discrete element models of debris avalanches: parameterisation and the reproducibility of experimental results. Journal of Geophysical Research Earth Surface 114.

Bourgeot, J.-M., Canudas-de Wit, C., Brogiato, B., 2006. Impact shaping for double support walk: from the rocking block to the biped robot.; In: Tokhi, M., Virk, G., Hossain, M. (Eds.), Climbing and Walking Robots. Springer, Berlin Heidelberg, pp. 509–516.

Bozzolo, D., Pamini, R., 1986. Simulation of rock falls down a valley side. Acta Mechanica 63, 1–4.

Bozzolo, D., Pamini, R., Hutter, K., 10–15 July 1988. Rockfall analysis – a mathematical model and its test with field data. landslides. Proceedings of the Fifth International Symposium on Landslides, pp. 555–560.

Calvetti, F., Crosta, G.B., Tarella, M., 2000. Numerical simulation of dry granular flows: from the reproduction of small-scale experiments to the prediction of rock avalanches. Rivista Italiana di Geotecnica 2.

Campbell, C., Cleary, P., Hopkins, M., 1995. Large-scale landslide simulations – global deformation, velocities and basal friction. Journal of Geophysical Research - Solid Earth 100.

Chau, K., Wong, R., Liu, J., Wu, J., Lee, C., 1999. Shape effects on the coefficient of restitution during rockfall impacts. Proc. Of the 9th International Congress on Rock Mechanics: International Society for Rock Mechanics (ISRM), 1, pp. 541–544. Paris.

Cleary, P., Prakash, M., 2004. Discrete-element modelling and smoothed particle hydrodynamics: potential in the environmental sciences. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 362.

Cundall, P., 1987. Distinct element models of rock and soil structure. In: Brown, E.T. (Ed.), Analytical and Computational Methods in Engineering Rock Mechanics. Allen and Unwin, London, pp. 129–163.

Cundall, P., Strack, O., 1979. A discrete numerical-model for granular assemblies. Geotechnique 29 (1).

Denlinger, R., Iverson, R., 2004. Granular avalanches across irregular three-dimensional terrain: 1. theory and computation. Journal of Geophysical Research-Earth Surface 109.

Descoeudres, F., 7 nov 1997. Aspects gomcaniques des instabilités de falaises rocheuses et des chutes de blocs. Publications de la Socit Suisse de Mcanique des Sols et des Roches 135, montreux.

Evans, S., Hungr, O., 1993. The assessment of rockfall hazard at the base of talus slopes. Canadian Geotechnical Journal 30.

Favier, L., Daudon, D., Donze, F., Mazars, J., 2009. Predicting the drag coefficient of a granular flow using the discrete element method. Journal of Statistical Mechanics: Theory and Experiment 2009 (URL stacks.iop.org/JSTAT/2009/P06012).

Fornaro, M., Peila, D., Nebbia, M., 1990. Block falls on rock slopes – application of a numerical simulation program to some real cases. Proceedings of the 6th International Congress IAEG Rotterdam, NL, pp. 2173–2180.

Giani, G., Giacomini, A., Migliazza, M., Segalini, A., 2004. Experimental and theoretical studies to improve rock fall analysis and protection work design. Rock Mechanics and Rock Engineering 37 (5).

Habib, P., 20–21 mai 1976, 1977. Note sur le rebondissement des blocs rocheux. Proceedings Meeting on Rockfall Dynamics and Protective Works Effectiveness, pp. 123–125.

Hart, R., Cundall, P., Lemos, J., 1988. Formulation of a three-dimensional distinct element model – part 2: mechanical calculations for motion and interaction of a system composed of many polyhedral blocks. International Journal of Rock Mechanics Min. Sci. and Geomec. Abstr., 25.

Heidenreich, B., 2004. Small- and half-scale experimental studies of rockfall impacts on sandy slopes. Ph.D. thesis, EPFL, Lausanne, th'Â se de Doctorat.

Hungr, O., 1995. A model for the runout analysis of rapid flow slides, debris flows, and avalanches. Canadian Geotechnical Journal 32 (4).

Hungr, O., Evans, S., 2004. Entrainment of debris in rock avalanches: an analysis of a long run-out mechanism. Geological Society of America Bulletin 116 (9,10).

Jean, M., 1999. The non-smooth contact dynamics method. Computer Methods for Applied Mechanics and Engineering 177.

Johnson, K.L., 1987. Contact Mechanics. Cambridge University Press.

Linares-Guerrero, E., Goujon, C., Zenit, R., 2007. Increased mobility of bidisperse granular avalanches. Journal of Fluid Mechanics 593.

Luding, S., Tykhoniuk, R., Tomas, J., 2003. Anisotropic material behaviour in dense, cohesive-frictional powders. Chemical Engineering and Technology 26.

Mangeney-Castelnau, A., Vilotte, J., Bristeau, M., Perthame, B., Bouchut, F., Simeoni, C., Yerneni, S., 2003. Numerical modeling of avalanches based on Saint Venant equations using a kinetic scheme. Journal of Geophysical Research - Solid Earth 108.

Manzella, I., Labiouse, V., 2009. Flow experiments with gravel and blocks at small scale to investigate parameters and mechanisms involved in rock avalanches. Engineering Geology 109.

McDougall, S., Hungr, O., 2004. A model for the analysis of rapid landslide motion across three-dimensional terrain. Canadian Geotechnical Journal 41.

McDougall, S., Hungr, O., 2005. Dynamic modelling of entrainment in rapid landslides. Canadian Geotechnical Journal 42.

Moreau, J.J., 1994. Some numerical methods in multibody dynamics: application to granular material. European Journal Mechanics A/Solids 13 (4-suppl).

Oger, L., Savage, S., Corriveau, D., Sayed, M., 1998. Yield and deformation of an assembly of disks subjected to a deviatoric stress loading. Mechanics of Materials 27.

Okura, Y., Kitahara, H., Sammori, T., 2000a. Fluidization in dry landslides. Engineering Geology 56.

Okura, Y., Kitahara, H., Sammori, T., Kawanami, A., 2000b. The effects of rockfall volume on runout distance. Engineering Geology 58.

Paronuzzi, P., 2009. Rockfall-induced block propagation on a soil slope, Northern Italy. Environmental Geology 58. http://dx.doi.org/10.1007/s00254-008-1648-7.

Pfeiffer, T., Bowen, T., 1989. Computer simulation of rockfalls. Bulletin of the Association of Engineering Geologists 26 (1).

Pirulli, M., 2009. The Thurwieser rock avalanche (Italian alps): description and dynamic analysis. Engineering Geology 109 (1,2).

Pirulli, M., Mangeney, A., 2008. Results of back-analysis of the propagation of rock avalanches as a function of the assumed rheology. Rock Mechanics and Rock Engineering 41 (1).

Radjai, F., Richefeu, V., 2009. Contact dynamics as a nonsmooth discrete element method. Mechanics of Materials 41.

Salcarini, D., Conversini, P., Tamagnini, C., 2009. A comparison between numerical approaches for rockfall analysis. EGU2009-12025-1, 2009 EGU General Assembly 2009, Vol. 11, p. 12025.

Sautier, C., Labiouse, V., Pirulli, M., Scavia, C., Zhao, J., June 15–18 2010. Numerical simulation of gravel unconstrained flow experiments: a comparison between dan-3D and rash-3D codes. European Rock Mechanics Symposium EUROCK 2010, pp. 571–574.

Savage, S., Hutter, K., 1989. The motion of a finite mass of granular material down a rough incline. Journal of Fluid Mechanics 199.

AugSomfai, E., Roux, J.-N., Snoeijer, J.H., van Hecke, M., van Saarloos, W., 2005. Elastic wave propagation in confined granular systems. Physical Review E 72, 021301.

Staron, L., 2008. Mobility of long-runout rock flows: a discrete numerical investigation. Geophysical Journal International 172.

Taboada, A., Estrada, N., 2009. Rock-and-soil avalanches: theory and simulation. Journal of Geophysical Research 114. http://dx.doi.org/10.1029/2008JF001072 (URL F03004).

Thornton, C., Cummins, S., Cleary, P., 2011. An investigation of the comparative behaviour of alternative contact force models during elastic collisions. Powder Technology 210.

Tommasi, P., Campedel, P., Consorti, C., Ribacchi, R., 2008. A discontinuous approach to the numerical modelling of rock avalanches. Rock Mechanics and Rock Engineering 41 (1). http://dx.doi.org/10.1007/s00603-007-0133-z.

Tsuji, Y., Tanaka, T., Ishida, T., 1992a. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe. Powder Technology 71.

Tsuji, Y., Tanaka, T., Ishida, T., 1992b. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe. Powder Technology 71 (3), 239–250.

Urciuoli, G., 10–11 nov 1988. Sperimentazione sulla caduta di blocchi lungo un pendio nella formazione calcareo-dolomitica della penisolaconvengo sul tema: Cartografia e monitoraggio dei movimenti fransi, Bologna.

Valentino, R., Barla, G., Montrasio, L., 2008. Experimental analysis and micromechanical modelling of dry granular flow and impacts in laboratory flume tests. Rock Mechanics and Rock Engineering 41 (1).

van den Bergen, G., Oct. 2003. Collision Detection in Interactive 3D Environments (The Morgan Kaufmann Series in Interactive 3D Technology). Morgan Kaufmann.

Voellmy, A., 1955. Über die zerstörungskraft von lawinen. Schweizerische Bauzeitung 73.

Walton, O.R., Braun, R.L., 1986. Viscosity, granular temperature, and stress calculations for shearing assemblies of inelastic, frictional disks. Journal of Rheology New York 30.

Wu, S., 1985. Rockfall evaluation by computer simulation. Transportation Research Record 1031.

Wu, C., Thornton, C., Li, L., 2003. Coefficients of restitution for elastoplastic oblique impacts. Advanced Powder Technology 14.