Evaluation of the Impact Force of Dry Granular Flow onto Rock Shed

In the design of rock sheds for the mitigation of risk due to rapid and long landslides, a crucial role is played by the evaluation of the impact force exerted by the flowing mass on the rock sheds. This paper is focused on the influencing factors of the impact force of dry granular flow onto rock shed and in particular on the evaluation of the maximum impact force. The coupled DEM-FEM model calibrated with small-scale physical experiment is used to simulate the movement of dry granular flow coupled with impact forces on the rock-shed. Based on the numerical results, three key stages were identified of impact process, namely startup streams slippery, impact and pile-up. The maximum impact force increases linearly with bulk density, and the maximum impact force exhibits a power law dependence on the impact height and slop angle respectively. The sensitivities of bulk density, impact height, and slope angle on the maximum impact force are: 1.0, 0.496, and 2.32 respectively in the benchmark model. The parameters with high sensitivity should be given priority in the design of the rock shed. The results obtained from this study are useful for facilitating design of shed against dry granular flow.


Introduction
The mountainous areas in Southwest China, where the topography is complex, feature seriously weathered rock masses and have experienced a large number of landslides. These landslides provide abundant materials for the initiation of dry granular flow. When the earthquake or heavy rainfall occurs, the landslide will collide with the mountain and slide down the slope, leading to a dry granular flow (Chen & Zhang, 1994;Zhu, Wang, & Tang, 2000). The formed dry granular flow, with high speed and large displacement, can greatly threaten the safety of the residents and the smooth traffic flow. For instance, the Wenchuan earthquake caused more than 15,000 geo-hazards in the form of dry granular flows which resulted in about 20,000 deaths in 2008 (Yin, Wang, & Sun, 2009). Countermeasures have been made to minimize the dry granular flow's risk to downstream residential areas or transportation routes. There are mainly two kinds of protection structures that used to minimize this hazard: active ones (like nets) and passive ones. As active ones is hard to carry out because of avalanches' potential source area is difficult to figure out, engineers and researchers usually choose the passive ones. Rock sheds are regarded as passive protection structures and are widely used to protect against mountain hazards such as dry granular flow, due to its unique edge in terms of low construction cost and strong constructability in complex areas (Pei, Liu, & Wang, 2016;Kawahara & Muro, 2006;Mommessin, Perrotin, & Ma, 2012). Most rock sheds are made of concrete, and have a shock-absorbing layer such as sands on top of the structure (Montani, Descoeudres, & Labiouse, 1996;Kishi & Konno, 2003).
Up to now, the design of such structures takes into account only the impact of an individual rock block (Kishi & Konno, 2003;Calvetti, 2011;Delhomme, Mommessin, & Mougin, 2005;Wang, Zhou, & Luo, 2017). Though the standardized design of rock sheds under a single block impact has accumulated rich engineering experience (Montani, Descoeudres, & Labiouse, 1996;Kawahara & Muro, 2006), it cannot be applied to the design of rock sheds impacted by dry granular flow, due to totally different dynamic mechanical characteristics. To date, No firm guidelines built upon sounded theoretical basis are available for the design of rock shed impacted by dry granular flow. Therefore, further researches on the dynamic behavior of a rock shed impacted by a dry granular flow are urgently required. However, due to the estimation of the impact force exerted by dry granular flow is a prerequisite parameter for shed design, so it is necessary to study the influencing factors of the impact force, so as to provide references for facilitating design of shed against dry granular flow.
Physical modeling has been widely used in geotechnical engineering research because of its excellent controllability in testing conditions and good reliability of testing results. For instance, using indoor experimental methods, Jiang et al. investigated the impact of dry granular flow against a rigid retaining wall by calculating the impact force (Jiang & Towhata, 2013;Jiang, Zhao, & Towhata, 2015). Jiang designed a set of experiments to investigate the impact mechanism of dry granular flow against a curved rock shed (Jiang, Wang, & Son, 2018).
Thus, the research results can serve as a significant reference for practice engineering. A quantitative analysis of impact force is eagerly needed in the design.
In engineering practices, several semi-empirical methods have been used to es-  (Shen, Zhao, & Zhao, 2018). These studies are useful for providing us with ideas.
Nevertheless, these available methods still have the difficulties in estimating the impact force of dry granular flow on rock sheds. This is because each method was obtained in specific impacting and boundary conditions with strong assumptions, such that they cannot be generalized for wider applications. In addition, these methods fail to consider the influence of dry granular flow-rock shed coupled interaction. In the present study, a robust numerical tool is a good choice. As the dry granular flow is a collection composed of a large number of discrete particles, DEM is an effective method for studying dry granular flow. Lo et al. used the PFC-3D software to study the maximum impact of the rock shed suffered by the dry granular flow (Lo, Lee, & Lin, 2016). Bi et al. studied the optimization of buffer layer under the impact of dry granular flow by two-dimensional discrete element software, and obtained the optimal thickness of the buffer layer (Bi, He, & Li, 2016). However, the discrete element method is not suitable for investigations of disaster-structure coupled interaction.
The above researches generally focus on the research of a single factor of uncoupled impact force. Unfortunately, the coupled dynamic interaction between dry granular flow and a rock shed is very complicated because it depends on the kinematics of dry granular flow (like solid mass and velocity), the stiffness and geometrical characteristics of the rock shed. So a quantitative analysis of these conditions on impact force is eagerly needed in the design.
The DEM has been widely used for numerical modeling of rock avalanches (Lo, Lee, & Lin, 2016;Bi, He, & Li, 2016;Cundall, 2008). It is an appropriate tool for modeling rock avalanches because of the discrete nature of materials involved in these phenomena. On the other hand, the FEM, based on continuum mechanics theory, has been well developed. Stress-strain development path and failures of elements are easy to simulate by FEM. Therefore FEM is a highly suitable method to model the rock shed (Albaba, Lambert, & Kneib, 2017

DEM Modeling of Particle System
The DEM is employed to model the particle system in a dry granular flow. It is assumed that the particles are all elastic soft spheres of different sizes. The use of spherical particles in DEM simulations will inevitably lead to a soil structure different from that of real natural soils with a reduced granular internal friction.
However, through careful model calibrations, an assembly of spherical particles with proper mechanical and physical properties can still be used to simulate the behavior of debris flows. This setting simplifies the complexity of dry granular flow but is also able to deliver a realistic simulation of the interaction between obstacle and flow (Bi, He, & Li, 2016;Cundall, 2008;Karajan, Han, & Teng).

Governing Equation
The motion of discrete elements is governed by the second Newton's law, and there are one or more forces acting on each element. The distribution and evolution of the system are described through the motion and state change of each element in the system (Cundall, 2008;Albaba, Lambert, & Kneib, 2017;Karajan, Han, & Teng, 2014). For element i: where g is the gravitational acceleration. i m , i u  , i I and i θ  are the mass, translational acceleration, rotary inertia and rotational acceleration of element i respectively. n,ik f , t,ik f and ik T are the normal contact force, the tangential contact force and the torques of element i acted by its neighboring element k respectively. ik T can be obtained by formula , and ik l is the arm vector of the force to the center of element i.

Evaluation of Contact Forces
Particles in the simulations are interacting with a linear spring-dashpot contact (LSD) law with Coulomb failure criterion, which is simple and computationally efficient compared to Hertz contact model (Karajan, Han, & Teng, 2014). The contact model of two particles is shown in Figure 1. The overlap δ of two particles is calculated as follows: where i r and k r are the radius of particles i and k respectively. i x and k x are the position vector of particles i and k respectively.
The normal contact force n,ik f between interacting particles is calculated as follows: where n k , n c , δ  , and n are the normal spring stiffness, normal damping coefficient, the relative normal velocity and the unit normal displacement vector respectively.
The tangential contact forces t,ik f between interacting particles are calculated as follows: where t k , t c , t δ and µ are the tangential spring stiffness, tangential damping coefficient, the incremental tangential displacement, and the friction coefficient respectively. t k is taken as 2/7 n k (Albaba, Lambert, & Kneib, 2017;Karajan, Han, & Teng, 2014). n k is calculated as follows (Karajan, Han, & Teng, 2014): where k n is a stiffness proportionality constant. i κ and k κ are the bulk modulus of particle i and k respectively. E and ν are the elastic modulus and poisson ratio of particle respectively. n c and t c are calculated as follows (Karajan, Han, & Teng, 2014): where n η and t η are the normal damping ratio and tangential damping ratio of particles respectively.

Coupled DEM-FEM Model
The coupled governing equations are given by Equation (7). The first and second conditions refer to the governing equations of DEM. The final condition gives the governing equation of FEM.
n, t, n, t, 1 1 and ij T are the normal contact force, tangential contact force and the torques of discrete element i acted by its neighboring finite element j respectively. ij T can be obtained by formula , and ij l is the arm vector of the force. M, C, and K are the mass matrix, damping matrix and stiffness matrix of system respectively. X is the displacement of finite element node. a f and b f are the external force vector of finite elements and the contact force vector of between finite elements and discrete elements respectively.
The interaction between contact surfaces is handled following the penalty method. As before-mentioned, the combined finite-discrete element method proposed in this paper is focused at dynamic simulation, and the Central Difference Method (CDM) is employed to solve Equation (7). Since CDM is conditional convergence, the step must satisfy the numerical stability conditions. Both DEM and FEM adopt the conditional stable central difference method, and their coupling requires that their integrals must be synchronized, which requires both to adopt the same time step under the same calculation framework. The time step DEM-FEM t ∆ takes the smaller value of both (Karajan, Han, & Teng, 2014).
c is the material sound speed. β is the scaling coefficient of time step length. m is the particle mass. spring K is the contact spring stiffness of particles, and min L is the minimum finite element size.

Numerical Model
Due to the existence of inter-particle porosity, the density of sand is set to 2800 kg/m 3 through the numerical volume test (Adrian Jensen, Kirk Fraser & George Laird, 2014), so that the bulk density of the initial debris deposition could be guaranteed to be 1350 kg/m 3 . This test allows the analyst to adjust the bulk density. The normal damping ratio between particles is set to 0.7. The tangential damping ratio is set to 0.4. The stiffness proportionality constant is set to 0.01. These three values were obtained by trial and error, so that the overall numerical results of debris dynamics can match the experimental observations in the model validation process. The particle Young's modulus and Poisson's ratio are set according to the commonly used values in numerical simulations of granular medium, as listed in Table 2. The initial debris deposition is composed of an assembly of 4968 randomly distributed spherical particles. As a channel, the side wall and bottom wall have little influence on the test, and so they are modeled as rigid wall. In the experiment, the load cells upon impact have a very small normal strain, and so the retaining wall is simulated by elastic wall. The material parameters are shown in Table 2

Model Validation
It can be observed that the numerical results can match well the experimental measurements (see Figure 3). Both of them are very similar in granular deposition shape, the length error of the granular deposition shape along the bottom direction of the sand trough is 14.3%, and the length error of the granular deposition shape along the height direction of the barrier baffle is 10%. A static pressure dead zone is formed at the baffle (see Figure 3(a)). In particular, it is apparent that the numerical simulation result can represent the general trend of the impact force evolution. The maximum impact force of test and simulation is 788.6 N and 824.4 N respectively, and the error is no more than 4.54% (see Figure 3(b)). The numerical simulation can capture the characteristics of the peak force observed in experiments.

Parametric Study
A comparison between the experimental and the numerical results indicates that the coupled DEM-FEM experiment adopted in this study can well simulate the

Problem Geometry and Numerical Model
The model of rock shed impacted by a dry granular flow is established by placing the retaining wall at the horizontal plane. For the simplified model, the rectangularly shaped debris flow material may not have the same impact energy compared with the actual case; however, this paper aims to study the regular variation of impact energy for qualitative analysis rather than quantitative examination. Simplifying the model makes the analysis simpler and easier. The geometric scheme adopted for this study is shown in Figure 4, which identifies the key parameters, including the bulk density of dry granular flow ( ρ ), impact height of dry granular flow (H) and slope angle ( θ ). In Table 3, numerical values are assigned to allgeometric parameters used here. In this study, the scale shed model aims to quantitatively study the variation of impact force with these key parameters. Because this paper aims to study the forces on the slab of the shed structure, only the slab is taken for coupled analysis. The slab is fixed with four corners. The slab here is simulated with elastic shell elements. Their material parameters  Table 2. The buffer layer covering the slab is composed of gravel cushion layer. The gravel parameter studies using DEM have been extensively conducted by many scholars (Chaplot, Walter, & Curmi, 2000;Takahara & Miura, 1998). Takahara et al. used the DEM to simulate the gravel material and found that DEM reflected the physical properties better (Takahara & Miura, 1998). In this paper, the particle diameter is randomly and uniformly set as 8 -10 mm. For convenience, other material parameters are the same as dry granular flow.

Impact Process
In this section, the general features of granular flow impacting on a rock shed of ρ = 1350 kg/m 3 , H = 2.5 m and θ = 60˚ are illustrated. Figure 5 shows the evolution of granular profiles during the impacting process. Based on the numerical results, the evolutions of granular flow deposition can be divided into three key stages, namely the startup streams slippery, impact and pile-up. In the startup streams slippery stage (0 -0.78 s), it can be observed that, with the increase of time, the speed difference between front and tail increases. The granular flow consists of three parts: front, middle and tail, and the flow pattern of granular flow is basically formed. When t = 0.78 s, the velocities of the three parts of the particle flow are 6.4 m/s, 3.8 m/s and 2.6 m/s in the front, middle and tail respectively. In the impact stage (0.78 -1.1 s), it can be observed that, the cushion layer is highly deformed. The buffer layer dissipates the kinetic energy of the granular flow and reduces the impact force. In the pile-up stage (0.78 -1.1 s), the granular flow particles begin to slow down, and the shape of the granular flow becomes gradually convex with respect to the shedslab. The final deposition is shown in Figure 5 (t = 3.0 s). Only a small portion of the granular flows stays on the shed slab, and most of them will eventually accumulate on the ground after sliding out of the shedslab. The granular flow is basically in a stable state and accumulates in the dead zone, which can act as a buffer layer against the next flow shock. Figure 6 shows the effect of bulk density on impact force. It can be seen from Figure 6(a) that particles of different densities reach their maximum impact force at almost the same time, and the whole impact process last about 0.5 s. Figure 6(b) gives a strong linear correlation between maximum impact force and bulk density. A popular formula for estimation of the maximum impact force is based on the well-known hydrodynamic model which shows a positive linear correlation between the maximum impact force and the density (Kwan, 2012). This also indicates the common character of granular impact: the maximum impact force has a strong linear correlation with the bulk density. Figure 7 shows the effect of impact height on impact force. As shown in Figure 7(a), the maximum impact force increases with the increase of impact height. Because of the increasing of the impact height, the gravitational potential Liu et al. energy of dry granular flow increases, leading to the increasing of the velocity of the dry granular flow front. Figure 7(b) gives a strong power function correlation between the maximum impact force and the impact height. The rate of increment about the maximum force declines as the impact height ascends.  Figure 8 shows the effect of slope angle on impact force. As shown in Figure   8(a), as the slope angle increases, the maximum impact force increases and the impact duration decreases. This is because under the same impact height, the overall potential energy of the dry granular flow is the same. The larger slope angle, the smaller the velocity difference between the front and tail, the shorter flow length of the flow, and the larger average flow depth. On the other hand, with the increasing of slope angle, the friction energy consumption reduces and the front velocity of the flow increases, leading to the greater impact force at the shed slab. Figure 8(b) gives a strong power function correlation between maximum impact force and slope angle. The rate of increment about the maximum force ascends as the slope angle ascends.

Sensitivity Analysis
The sensitivity analysis method in the system analysis can be used to determine the main influence parameters. The relationship between the parameters of the dry granular flow and the maximum impact force is respectively fitted, and Equation (9) is obtained.

R E T R A C T E D
According to the sensitivity calculation formula 10 (Zhang & Zhun, 1993): is the sensitivity function of sensitivity parameter k a , which is the fitting function in this paper. * k a is the reference value of sensitivity parameters. * U is the value of sensitivity function when . The most important factor affecting the maximum impact force is slope angle, followed by bulk density, and the impact height has the least effect on the maximum impact force. The parameters with higher sensitivity value should be considered in the rock shed design.

Conclusion and Discussion
The impact of a dry granular flow on a rock shed has been investigated by a novel numerical framework. A coupled DEM-FEM approach is employed in this framework. The coupled DEM-FEM model was successfully verified by indoor test results. As such, the validated model was then employed to investigate the evaluation of the maximum impact force of dry granular flow against rock shed under different influencing factors. The key findings from this study are summarized as follows: Based on the numerical modeling, three key stages during impact process, namely the startup streams slippery, impact and pile-up were identified.
Certain outcomes were discussed with particular emphasis on the influences of bulk density, impact height, and slop angle on the impact forces exerted on the rock shed. The maximum impact force increases linearly with bulk density.
The maximum impact force increases in the form of a power law with the increases of the impact height, and its power index is less than 1. The maximum impact force increases in the form of a power law with the increases of the slope angle, and the power index is greater than 1. The sensitivities of bulk density, impact height, and slope angle on the maximum impact force are: 1.0, 0.496,

R E T R A
C T E D Journal of Geoscience and Environment Protection and 2.32 respectively in the benchmark model. The parameters with high sensitivity should be given priority in the design of the rock shed.
However, this paper remains a rather preliminary pilot study. Only the maximum impact force of the dry granular flow on the rock shed is taken as the dependent variable, and the response to the internal force of rock shed needs to be further modeled and analyzed. On the other hand, the effects of the buffer layer and initial shape of granular flow on the maximum impact force need to be further study.