Numerical Study of Rockfalls on Covered Galleries by the Discrete Element Method   Wen-Ji Shiu Ph.D. Student, Laboratoire Sols, Solides, Structures (UMR 5521) Domaine Universitaire - B.P. 53 - 38041 Grenoble cedex 9 - France e-mail: wen-jie.shiu@hmg.inpg.fr Frédéric Victor Donzé Professor, Laboratoire Sols, Solides, Structures (UMR 5521) Domaine Universitaire - B.P. 53 - 38041 Grenoble cedex 9 - France e-mail: Frédéric.Donze@geo.hmg.inpg.fr Sophie-Adélaïde Magnier Application Engineer, GEONUM, France e-mail: sam@geonum.com

Abstract

To contribute to the validity of a laboratory-based experimental formulation proposed in the design of soil-covered concrete sheds, dynamic impacts have been numerically simulated with a three dimensional  Discrete Element Method. The proposed numerical model shows the capability of the Discrete Element Method not only to be good at reproducing the experimental tests performed in the Swiss Federal Institute of Technology of Lausanne, but also to be used in a predictive manner.

Parametric studies were carried out to check the influence of local parameters. Doing so, it has been shown that the stiffness of the soil layer has a major effect on the transmitted impact force. In addition, to the use of basic elastic frictional behavior laws, a damping term must be added. Doing so, it has been seen that the influence of this damping term can dramatically change the numerical results and it must be used with caution. This study is a contribution to the founding of methods which determine the required dimensions of materials used for shock absorption when designing structures subjected to dynamic loadings.

Keywords:  Discrete Elements Method, numerical modeling, rockfall, dynamic loading, rockshed

Introduction

To prevent great damage on traffic roads due to hazardous rockfall, heavy sheltering structures have been designed. These structures are usually made of reinforced concrete slabs supported by numerous piles and covered by a shock absorbing cushion of granular material. To improve the effectiveness of these structures, laboratory-based experimental formulation have been proposed by the Rock Mechanics Laboratory of the Swiss Federal Institute of Technology of Lausanne (Montani, 1998).

These laboratory studies dealt with systematic and very wide experimental parametric tests which gave an important range of values for the dynamic loading on the surface of the concrete slab covered by the dissipative layer.

Beside these experimental tests, different numerical models were used to validate the observations and to help in the formulation of the dimensioning analytical expressions for extreme cases, which are not considered by the experiments (Calvetti et al., 1997, Donzé et al., 1999). Both of these used the discrete element approach, one using a 2D commercial code PFC^2D and the other a 3D research code SDEC (Donzé and Magnier, 1997).

In order to increse the validity of the numerical results obtained with SDEC, a new series of 3D numerical simulations have been carried out with the commercial code PFC^3D (Itasca, 2003). To do so, the numerical results are compared with the experimental results with the following protocol:

• calibrate the reference numerical model on a specific test which is the rock drop at a height of 10 meters,
  
• and use this reference model, without modifying its intrisic parameters, to predict the loading of the concrete slab and the penetration of the impacting block in the granular material for different drop heights.

Then, the numerical model is used to study the effect of the rigidity of the soil layer, the coefficient of friction between grains and the damping coefficient.

Experimental Study

The present study is based on impact experiments performed at the Rock Mechanics Laboratory of the Swiss Federal Institute of Technology of Lausanne (Labiouse et al., 1994). In these experiments, performed in a well which has a diameter of 5 m and a depth of 8 m, weights were dropped from various heights on the bottom of the well which had been covered by a compacted soil cushioning layer. The bottom of the well is a 0.9 m thick slab of reinforced concrete which is considered to be infinitely rigid. It is instrumented with 5 pressure gages. The weights were metallic cylinders with a spherical head filled with concrete. The test data used for calibration were obtained with a 500 kg weight falling vertically from a height of 10 m onto the bottom of the well covered with a 0.5 m thick soil layer. In the other laboratory experiments that will serve as comparisons, the weights were kept at 500 kg and the drop heights varied from 1 m to 10 m.

The measurements available from these tests are the tangential deceleration of the block, the pressures occurring at different places on the rigid concrete slab and the penetration of the blocks in the soil cushion. From these measurements and after some data processing, values of velocity, theoretical block penetration trajectories, dynamic loading by acceleration of the block, dynamic loading by the integration of the values recorded by the pressure gages and impact duration are obtained.

The Discrete Elements Method Used

The PFC^3D code (Particle Flow Code in 3 Dimensions, Itasca, 2003) is used to model the experimental tests. This code is based on the discrete Elements Method which uses spherical "Discrete Elements" ("DE", also called "particles"). These DE move independently of one another and interact only at contacts or interfaces between the particles. The particles are assumed to be non deformable and the behavior of the contacts is characterized by a soft contact approach. The mechanical behavior of such a system is described in terms of the movement of each particle and the inter-particle forces acting at each contact point. Newton's laws of motion provide the fundamental relationship between particle motion and the forces causing that motion (Itasca, 2003). The equations of motion applied to each elements are defined by Equations 1, 2 and 3:

where F_i is the contact force, K_i the stiffness associated to each element, with k_n in the normal direction and k_s in the tangent direction, m is the mass of each element,  and   are the translational acceleration and rotational acceleration respectively, g_i is the gravitational acceleration, M_i is the resultant moment acting on each element and I is the moment of inertia. During the calculation cycle, the force-displacement law (equation 1) is calculated first, then the new element’s position will be updated by the law of motion (equation 2, equation 3).

The slip model is defined by the friction coefficient at the contact µ [dimensionless], where µ is taken to be the minimum friction coefficient of the two contacting entities. The criterion of zero-normal strength is enforced by checking whether the overlap is less than or equal to zero. If it is the case, then both the normal and shear contact forces are set to zero. The contact is checked for slip conditions by calculating the maximum allowable shear contact force:

If  ,  the two elements slip.

Energy dissipation was also used in the numerical model. The energy involved between two interacting elements is dissipated through frictional sliding for which a Coulomb friction coefficient µ is defined. Moreover, a local non-viscous damping is available in PFC^3D, where the damping force is put together with the equation of motion such that,

 

  where F(i), M(i), and A(i) are the generalized force, mass and acceleration components respectively, and  is the damping force

 

  where α is the numerical damping ( the detailed description can be found in the PFC^3D manual, Itasca, 2003). After some pre-process numerical simulation tests, the numerical damping factor is set to 0.05 and 0.3 for the block elements and the soil elements respectively.

Set up of the numerical model

The goal is to model a structure, in which some of the macroscopic material properties (Young's modulus, Poisson's ratio, friction angle) are known. The structure’s geometry is discretized with a collection of discrete elements. To each of these elements a set of  local parameters is assigned so that the macroscopic behavior of this collection is representative of the real material. This procedure is based on the simulation of quasi-static uniaxial compression tests. A tri-axial test model is pre-developed in PFC^3D and for  a standard-sized specimen:

• a compact, polydisperse discrete element collection is generated,
• an elastic compression test is run with local elastic parameters given by the "macro-micro" relations,
• compressive axial tests are simulated to deduce the remaining local parameters.

By performing these tests, the local parameters k_n, k_s,  are set such that the global mechanical properties of the collection of discrete elements are as close as possible, to those of a soil with a 1,2 GPa Young’s modulus and a Poisson's ratio of 0,27 .
Then, when simulating the first impact test, some readjustements were needed to fit the experimental data set. This readjustment procedure was performed only once and the exact same set of parameters was used in all subsequent impact test cases to demonstrate the predicting capability of the method. The input numerical data are given in Table 1.

Table 1. Parameters of the numerical model

kn (N/m)	107
ks (N/m)	5x106
Numerical damping (soil)	0.3
Numerical damping (bloc)	0.05
Coefficient of friction (soil)	0.7°
Coefficient of friction (soil)	0.5°
DE radius (cm)	4-6

The 0.5 m soil thickness experiment is modeled by the code PFC^3D which contains two types of elements: the wall elements and distinct elements. To model the concrete bottom slab and the four edges containing the granular material, five walls elements were used. This simple way of modeling the boundary conditions is chosen, because we are just interested in computing the transmitted effort to a medium for which the deformation is neglected.

Then, the soil layer is represented by spherical elements which radii are distributed over an interval of size set between 0.04 m - 0.06 m. The size of the distinct elements soil layer is 3.4 m x 3.4 m x 0.5 m. A total of 5592 elements is obtained. The impacting block is modeled by only one spherical element, which has a radius of 0.36 m, which corresponds to a rock with a radius of a sphere of 500 kg. The initial position of the block is such that its lower part is just above the surface of the soil layer. The initial position of the center of the block is given by: , where R_b is the radius of the block and e_s is the thickness of the soil layer. The block has an initial velocity , where g is the acceleration of gravity and h is the drop height (Fig.1).

Figure 1. Numerical setup

Vertical impact test

The distribution of the contact forces inside the soil layer represented by the set of distinct elements before and after impact is shown in Figure 2. Before impact, only gravity generates the contact forces which tend to increase vertically (Fig. 2a). Then, after impact i.e. when the slab loading reaches a maximum value, the distribution of the contact force between the elements is particularly concentrated around the point of impact. This distribution which tends to decrease in the radial direction is well reproduced by the model (Fig. 2b).

The distribution of the transmitted forces on the concrete slab by the soil layer is given in Figure 3. These forces are compared to the integrated forces measured in the experiments by sensors placed on the surface of the concrete slab (Fig.4) and the results shown in Figure 5 correspond to the sensors along to the vertical direction (P3, P4, ..., P7).

The closer the sensors are from the impact location, the higher the loading force is: qualitatively, the numerical results show the same evolution (Fig.3).

Figure 2. Contact forces distribution in the distinct elements set before and after impact

Figure 3. Distribution of the loading contact forces on the bottom slab

Figure 4. Disposition of the pressure sensors for the experimental case (after Montani, 1998)

Figure 5. Loading forces on the slab integrated by the sensors for the experimental case (after Montani, 1998)

Quantitative comparison between numerical and experimental results

A 500 kg block falling from a 10 m height on an uncompacted soil layer of 0.5 m in thickness was used as the reference case. The parameters used in the numerical model are given in Table 1. As in the experimental cases, once the numerical model is calibrated for this reference case, only the drop height of the impact block is changed .

For the four different drop heights, 1, 4, 7 and 10 m, the penetration of the block in the soil layer is directly recorded by tracking the vertical position of the block (Fig. 6) during its impact process. The reference vertical position is defined as the position where the first contact takes place between the impact block and the soil layer. Despite the fluctuation of the characteristic sizes of the distinct elements, which range from 0.08 to 0.12 m, the values of penetration obtained by numerical simulations are comparable with those estimated in the experiments.

Figure 6. Impact block penetration for drop height of (a) 10 m, (b) 7 m, (c) 4 m, (d) 1 m

The loading of the concrete slab is numerically measured by summing the forces applied on the bottom wall element by the spherical distinct elements in contact. The results are shown in Figure 7. These numerical results are compared with the experimental results on the same graph. The results are similar except for the 1 m drop height. However, Montani (Montani et al., 1997), notes that for a low drop height, the dispersion of the results is important. Then, the results dealing with the block velocity are also compared (see Figure 8). Here again the agreement with the experimental data is good. Finally, the block acceleration is recorded (see Figure 9). The numerical results are obtained by deriving the velocity data. Note that, even if the numerical velocity results are quite good, they are slightly higher than the experimental results at the time preceeding the immobilization of the block (Fig. 8). The derivatives of the velocity data increases this difference, and the numerical results of block acceleration are thus higher than the experimental results, the difference between the peak values of the two curves is about 15 to 20 %.

Thus, the good agreement between numerical and experimental results confirms the preliminary results obtained with the numerical code SDEC. It seems that once the model is calibrated and validated with a reference case (here, the 10 m drop height), the model can successfully predict the behaviors of the structures for various drop heights. It is then important to analyze the effect of some parameters used, such as the dissipative factor.

Figure 7. Slab loading for (a) 10 m, (b) 7 m, (c) 4 m, (d) 1m height drop

Figure 8. Block velocity for (a) 10 m, (b) 7 m, (c) 4 m, (d) 1 m height drop

Figure 9. Block acceleration for (A) 10 m, (b) 7 m, (c) 4 m, (d) 1 m height drop

Parametric study

Two parameters are involved in the contact law between the soil layer distinct elements: the stiffness and the coefficient of friction. The effects of these two parameters were previously studied (Donzé et al., 1999) and a validation of these effects is undertaken here.

Another parameter, which seems to play an important role, is the coefficient of the numerical damping acting on the inertial term (Itasca, 2003). Its physical meaning is not clear, but using purely frictional contact with an elastic component is not enough to dissipate the energy transmitted to the set of distinct elements. Thus, if no velocity effect is taken into account, this numerical damping needs to be considered. Let's check its influence on the results.

Influence of the stiffness

The stiffness is one of the most significant parameter in the numerical model, but as of now, there does not exist a well established direct relationship between the local and the global stiffness for a disordered medium. Using the reference formulation proposed by S. Hentz (Hentz et al.,  2004), the values leading to a good agreement with the experimental data are: for the normal stiffness k_n = 10⁷ N/m, and for the tangential stiffness k_s = 5x10⁶N/m.

Figure 10. Influence of the stiffness on the penetration, the acceleration of the block and loading of the slab

The influence of the stiffness on the penetration, the acceleration of the block and loading of the slab is shown in Figure 10. This parametric study has been carried out for two different tests in order to check the repeatability of the simulation. To analyze the influence of the stiffness, a one dimensional analogy can be made when the block falls on a spring with various values of rigidity. The larger the stiffness, the less the block penetrates. By increasing the stiffness, the impact energy of the block can be transmitted more easily to the concrete slab, which indues an increase of the slab loading. There is also a relation between the stiffness of the soil layer and the block acceleration, according to the equation . The larger the stiffness, the larger the pulse of the block, which implies an increase of the acceleration.

According to Figure 10, the contribution of the stiffness of the soil layer with respect to the penetration and the acceleration of the block and to the loading of the slab is really important.

Influence of the friction in the soil layer

The coefficient of friction plays an important role in the contact force as soon as the shear component becomes important. During the impact, when the friction coefficient is low, the distinct elements can move easily, which induces an increase of the penetration of the block in the soil with a decrease of its acceleration (Fig. 11). When the friction coefficient becomes important, the transmitted force increases if in a shearing mode which induces an increase of the loading of the slab (Fig. 12). However, the increase of the loading slab is not comparable to the one observed with the increase of the stiffness value. This shows that for this type of loading, the influence of the friction coefficient remains small because shearing is not the main process of the force transmission. In the present model, µ = 0.7.

Figure 11. Relation between Fn and Fs for the Mohr Coulomb criterion

Figure 12. Influence of the friction coefficient of soil layer

Influence of the numerical damping

We are interested in the influence of the numerical damping of the soil, while insisting on the fact that this inertial damping is used to compensate the lack of "representativity" of the model to local dissipation. This parameter is different from the two previous parameters (stiffness and friction coefficient) because it is not directly defined in the local behavior law. Moreover, the physical meaning of this damping parameter is not obvious. The damping used in code PFC^3D is a non-viscous damping, which is expressed in terms of a force in the motion equation, where it directly removes  a part of the global force applied to the distinct elements by mean of a damping coefficient (Itasca, 2003). For example, if this damping coefficient is a = 0.4, the global force applied to the elements decreases by 40 %, thus the contact force is weaker. Several tests have been carried out to simulate the best soil behavior, and the resulting value is a = 0.3 ~ 0.4.

Figure 13. Influence of the damping coefficient for the soil layer

Figure 14. Forces applied to the block during the contact.

According to Figure 13, when the damping coefficient acting on the contact law between the distinct elements representing the soil layer increases, the block penetrates less and the loading of the slab decreases, because there is less force transmitted to the soil. As shown in Figure 14, by increasing the coefficient, the apparent resulting force Fr due to the soil layer increases, (note that the gravity force of the block is kept constant), inducing an increase of the acceleration of the block.

Thus, even if this type of damping does not have a physical meaning, it is a parameter which can deeply modify the results obtained by the numerical model.

The influences of the parameters are summarized in the following table:

Table 2. Summary of the influence of the parameters

	Damping of the soil	Friction coefficient of the soil	Stiffness of the soil
Penetration of the block
Acceleration of the block
Slab Loading

conclusion

A numerical model using PFC^3D has been set up to simulate the block fall experiments on a concrete slab covered by a soil layer which were run at the LMS laboratory (EPFL). The results of this study show that after the calibration of the model for a reference case (corresponding to a drop height of 10 m), this model could predict the behavior for various drop heights even for parameter values inot available to the experimentation.

Parametric studies were made for the damping coefficient, the stiffness and the friction coefficient defined for the distinct elements representing the soil layer. The results show the influence of these parameters on the penetration and the acceleration of the block and the slab loading. It has been seen that:

• the friction coefficient does not have an important influence because the loading process is mainly compressive.
• At the opposite, the local stiffness has a major effect on the value of the transmitted impact force.

The damping coefficient used has not only a significant influence on the numerical results, but it must be used to get quantitative results. It has been seen that using a simple elastic-frictional local behavior law with spherical discrete element does not dissipate the impact energy sufficiently. However, the inertial based formulation proposed here has a very unclear physical meaning. An other option would be to use a dissipative term based on viscous damping formulation (Zhang and Withen, 1996, Leszczynski, 2004), but once again, its physical meaning remains unclear. A possible update in the local description would also be the use of a more realistic geometry to represent the granular material, for example the use of tetrahedra. It would involve the transfer of a moment between elements, with a possible local plastic limit, which could contribute to dissipate the impact energy. Finally, this dynamicaldamping problem requires a deep investigation.

A three-dimensional model is advantageous because all the parameters are directly available from the experiments and a large variety of parameters can be studied. However, the idenfication of the local parameters can be difficult. Moreover, using a 3D Discrete Element Model at the structure scale, is still expensive from a computational point of view. Despite that, the predicitve results obtained are quite promising.

References

1. Calvetti, F., R. Genchi, L. Nesta, and R. Nova (1997) Numerical simulation of rock block impacts on rigid sheds, S. Pietruszczak, G. N. Pande Eds., Numerical Models in Geomechanics, NUMOG VI, 635-640, Montreal, CA., Balkema Rotterdam Publishers.

2. Donzé, F. V., S. A. Magnier, S. Montani, F. Descoeudres. (1999) Numerical study of rock block impacts on soil-covered sheds by a discrete element method / Etude numérique de la sollicitation dynamique des galeries de protection lors des chutes de blocs par la méthode des éléments discrets, 5ème Conférence Nationale en Génie Parasismique, AFPS'99, Cachan.
   

3. Donzé, F. V., and S.-A. Magnier (1997) «Spherical Discrete Element Code» In: Discrete Element Project Report no. 2. GEOTOP, Université du Québec à Montréal (http://geo.hmg.inpg.fr/frederic/publications_FVD.html).

4. Hentz, S., L. Daudeville, F.V. Donzé (2004) Identification and Validation of a Discrete Element Model for Concrete,  ASCE Journal of Engineering Mechanics, ASCE Journal of Engineering Mechanics  vol:130 No 6 pp.709-719.
   

5. Itasca Consulting Group, Inc. (2003) PFC3D (Particle Flow Code in 3 Dimensions), Version 3.0. 656 pp.

6. Leszynski, J.S. (2004) Using the fractional interaction law to model the impact dynamics in arbitrary form of multiparticle collisions, Physical Review B, E70,051315.

7. Labiouse, V., F. Descoeudres, S. Montani, and C.-A. Schmidhalter (1994) Etude experimentale de la chute de blocs rocheux sur une dalle en béton armé recouverte par des matériaux amortissants, Revue Française de Géotechnique, vol.69, pp. 41-62.

8. Montani, S. (1998) Sollicitation dynamique de la couverture des galeries de protection lors de chutes de blocs, Ph. D. Thesis, Ecole Polytechnique Fédérale de Lausanne.

9. Montani, S., F. Descoeudres, and K.M. Bucher (1997) Numerical analysis of rock blocks impacting a rock shed covered by a soil layer, S. Pietruszczak, G. N. Pande Eds., Numerical Models in Geomechanics, NUMOG VI.

10. Zhang, D., and W.J. Whiten (1996) The calculation of contact forces between particles using spring and damping models, Powder Technlogy vol:88, pp. 59-64.

	© ejge, 2006