Home Search Collections Journals About Contact us My IOPscience Predicting the drag coefficient of a granular flow using the discrete element method This content has been downloaded from IOPscience. Please scroll down to see the full text. J. Stat. Mech. (2009) P06012 (http://iopscience.iop.org/1742-5468/2009/06/P06012) View the table of contents for this issue, or go to the journal homepage for more Download details: IP Address: 211.138.121.36 This content was downloaded on 30/09/2013 at 06:52 Please note that terms and conditions apply. J ournal of Statistical Mechanics: Theory and Experiment An IOP and SISSA journal Lionel Favier, Dominique Daudon, Frédéric-Victor Donzé and Jacky Mazars 3S-R Laboratory, UMR 5521, Domaine Universitaire, BP 53, 38041 Grenoble Cedex 9, France E-mail: lionel.favier@hmg.inpg.fr, dominique.daudon@hmg.inpg.fr, frederic.donze@hmg.inpg.fr and jacky.mazars@inpg.fr Received 30 October 2008 Accepted 10 February 2009 Published 16 June 2009 Online at stacks.iop.org/JSTAT/2009/P06012 doi:10.1088/1742-5468/2009/06/P06012 Abstract. Passive-protection structures, against snow avalanches, are designed with rough estimates of the drag coefficient, depending straightforwardly on the obstacle’s geometry. In this paper, assuming an avalanche as a dry granular flow, a numerical model is presented for studying the influence on the drag coefficient of both the shape and size of an obstacle impacted by a granular flow. Smallscale laboratory experiments were conducted to validate the numerical model. During the experiments, velocity profiles were estimated using the particle image velocimetry method applied to the pictures recorded by a fast camera, which focused perpendicular to the lateral wall. Flow thickness variations along the slope were estimated by post-processing the images recorded by another fast camera, which films a laser line projected on the free surface flow. From the lateral motion of that line, the thickness could be determined. The granular impact force was measured through an instrumented obstacle positioned at the lower end of a canal. A 3D numerical model, based on the discrete element method using the YADE code, was set up to reproduce the experimental configuration. The law of contact between discrete elements involved elastic components (normal and tangential stiffnesses) and dissipative components (a normal restitution coefficient and a friction coefficient based on the Coulomb friction law). The model was validated by comparisons with both the experimental flow characteristics (velocity profiles and thicknesses) and the impact load history. Once validated, the numerical model was used to investigate the contribution of the height and c 2009 IOP Publishing Ltd and SISSA 1742-5468/09/P06012+14$30.00 J. Stat. Mech. (2009) P06012 Predicting the drag coefficient of a granular flow using the discrete element method Predicting the drag coefficient of a granular flow using the discrete element method shape of the obstacle to the drag coefficient. Finally, results are discussed and compared with ones from other studies. Keywords: avalanches (experiment), discrete fluid models Contents 2 2. Experimental set-up 3 3. Numerical model description 3.1. Force–displacement law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. Micromechanical parameters . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3. Simulation methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 6 8 9 4. Model validation 9 5. Obstacle shape numerical influence on the drag coefficient 10 5.1. Size of the influence zone . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5.2. Mean velocity and solid fraction computation . . . . . . . . . . . . . . . . 11 5.3. Evolution of Cd . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 6. Conclusions and perspectives 13 Acknowledgments 13 References 14 1. Introduction The design of passive-protection structures, against snow avalanches, and generally against fast flows under gravity, is limited because of the lack of knowledge on the exerted pressure history. This leads to the use of common calculation methods, which do not correctly take into account the duration and intensity of this specific dynamic phenomenon [1]. Basically, the design of these structures uses the drag force (NF EN Eurocode, 2003), defined as Fd = 1 2 ρ Ux2 A × Cd , (1) where ρ is an estimation of the flowing snow density, Ux its velocity along the main slope, Cd the drag coefficient depending on the flow regime and obstacle geometry, and A the reference surface. Theoretically, Fd is the load applied by a perfect fluid in stationary flow on an immersed obstacle. Consequently, using equation (1) leads to a wide range of values for Cd , whose influence is a complex combination of flow regime [2], obstacle surface shape and size [3]. However, the Swiss historical guideline applied to snow avalanches [4] recommends that Cd varies in the range of 1–2 for small immersed obstacles. The first historical studies to investigate snow avalanche load against obstacles were based on theoretical analyses [5]: the load was found to be proportional to the squared velocity. Later, full-scale experiments were carried out [6, 7], and the results indicated that the avalanche load was dependent on both the squared velocity and the density in doi:10.1088/1742-5468/2009/06/P06012 2 J. Stat. Mech. (2009) P06012 1. Introduction Predicting the drag coefficient of a granular flow using the discrete element method 2. Experimental set-up The experimental canal (figure 1) is inclined 43◦ from the horizontal direction. It is 20 cm wide and 2 m long, divided into two 1 m parts: the upper one for the tank and the lower doi:10.1088/1742-5468/2009/06/P06012 3 J. Stat. Mech. (2009) P06012 the deposition zone. Afterwards, research was conducted to estimate the proportionality coefficient associated with this dependence. Because of the low quantity of data generally obtained by full-scale experiments, various researchers investigated this problem using small-scale experiments. These granular flows which hit obstacles of various sizes and shapes [3] highlighted their important effect on the drag coefficient. In the previously mentioned study, Cd was found to be 30% lower in the case of cylindrical obstacles compared to a square geometry. Moreover, Cd could be increased by a factor of two when the ratio between the obstacle height and the flow thickness was increased from 2 to 5. Some studies have quantified the dependence of Cd on the flow regime, which can be characterized by the Froude number, defined as Ux Fr = √ , (2) g hf where g is the gravity constant, hf is the flow thickness and Ux is the flow velocity. Fr is a dimensionless number, comparing inertial and gravitational forces within an avalanche flow. A recent study [8] reported that Fr < 1 for wet–dense avalanches and that Fr is in the range of 1–6 for dry–dense ones. Experimental studies in the laboratory on immersed objects in granular flow for Fr < 1 found Cd to be a function of Fr −2 [9, 10], which leads to a non-velocity-dependent load. Recent full-scale experiments [2] characterized Cd as a decreasing function of the Froude number with an empirical formula Cd = 10.8 Fr −1.3 for 0.8 < Fr < 6. Other full-scale measurements involving a 20 m high tubular instrumented pylon impacted by avalanches in various regimes [8], showed that the measured pressures are velocity independent if Fr < 1, but velocity dependent if 1 < Fr < 6. These studies show an increased interest in the drag coefficient. There are common results, such as Cd being a function of Fr −2 , leading to a velocity-independent load for a Froude number less than 1. Nevertheless, the flow’s velocity, density and thickness used to compute Cd are estimated at the free surface. The drag coefficient has to be related to the flow’s characteristics measured directly upstream from the obstacle, but in experiments, whether in the laboratory or in situ, inside flow velocities cannot be measured and density has to be estimated. The purpose of this study was therefore to propose numerical tools, based on discrete element modelling, for predicting the real drag coefficient. To do so, we first had to design small-scale laboratory experiments to calibrate the numerical model, which was used afterwards to predict the drag coefficient as a function of the impacted obstacle. In this paper, we first describe the experimental configuration, the measurement methods, and the numerical model based on the discrete element method (DEM). Then the model is validated by comparing experimental and numerical results for the velocity profiles, thicknesses and granular flow action. Finally, we numerically show two geometrical influences on the drag coefficient: the ratio between the obstacle height and the flow thickness, and the obstacle shape. Then, these results are compared with previous similar studies. Predicting the drag coefficient of a granular flow using the discrete element method one for the flow zone. The 15 l tank is filled with granular material before each experiment. The aperture of a manual gate releases the material, which then flows down under the gravity effect. A flat obstacle was added at the lower end of the canal to measure the granular impact force history. The obstacle height is ho = 4 cm, its width is lo = 4 cm. The granular material is a set of nearly spherical and quasi-monodisperse glass beads, whose mean radius is Rm = 2.466 mm with a scattering of δR = 2%. Two fast cameras gave flow pictures (figure 1), used afterwards to measure the velocity profiles (camera B) and thicknesses (camera A). There are several longitudinal positions for the cameras: with Ds the distance between the focused part of the canal and the obstacle, the ratio Ds /lo can be equal to 2.5, 5, 7.5, 10, 12.5. The obstacle is equipped with a load cell inside to obtain the normal force data, and a strain gauge at the back to obtain both normal and tangential data. The thickness is estimated at the centre of the canal, outside any effect from the lateral walls. A laser plane is projected downwards on the canal (figure 2). The method can be split into two steps. First, a calibrated piece of wood, made up of three 1 cm thick slices, is put down on the basal surface. As the laser plane meets the different slices of the piece, the line formed has a lateral displacement, which depends on the thickness of the piece. Second, during the flow, this relation is used to determine the thickness from the laser line displacement. The velocity profiles are estimated through the lateral wall or the basal wall. The particle image velocimetry (PIV) method is applied to n consecutive pictures, here taken from camera B. Each picture is discretized with analysis windows, which can overlap. Their size depends on the bead velocity, the position of the camera and the focal length. In these experiments, the size was found empirically. After the grid construction, we obtained mean velocities and standard deviations from each analysis window. doi:10.1088/1742-5468/2009/06/P06012 4 J. Stat. Mech. (2009) P06012 Figure 1. Frontal view and lateral view of the experimental set-up. The obstacle height is ho = 4 cm, its width is lo = 4 cm. The red points indicate various longitudinal positions for both cameras: 2.5 < Ds /lo < 12.5 with 2.5 steps. h(t) 3 cm 2 cm 1 cm 0 cm Camera A 3 cm 2 cm 1 cm 0 cm Predicting the drag coefficient of a granular flow using the discrete element method Laser plan Calibration Line displacement Figure 2. Central flow thickness measurement method. From the laser line lateral displacement, the thickness variations can be deduced. Figure 3. Experimental preprocessing and preliminary results. (a) Moving average applied to raw data from the load cell; the standard deviation is about 15%. (b) Similarity between load cell and strain gauge averaged curves. (c) Reproducibility of the experiments: arithmetic average on the force history for eight similar experiments. The experimental results using which we described the measurement methods are shown in section 4, where experimental and numerical results are compared. Here the preprocessing applied to the raw data and some preliminary results are detailed. First, because of the mean size of the granular beads, which are fairly large in relation to the obstacle side length (eight times the mean diameter), the data points of the impacted force against the obstacle are scattered. In order to use them later for comparison with numerical data points, we applied a moving average to signals (figure 3(a)). The corresponding standard deviation is about 15% around the averaged signal. The average curve can be split into three parts. The first and the last correspond, respectively, to a gradual increase and decrease in the granular action. The central part is characterized by a steady force, which is used to make comparisons. Second, comparing load cell and strain gauge signals for one experiment, it can be observed that they are temporally close (figure 3(b)). Therefore, we consider that the tangential action against the obstacle can be neglected; consequently the strain gauge is no longer used. Third, figure 3(c) highlights the good reproducibility of these experiments. This graph represents various statistical doi:10.1088/1742-5468/2009/06/P06012 5 J. Stat. Mech. (2009) P06012 Laser line moving Predicting the drag coefficient of a granular flow using the discrete element method properties of eight averaged force histories: the arithmetic average, which is equal to 15 N, and the standard deviation, which is not more than 3.5% for the entire history. 3. Numerical model description The experimental set-up was numerically reproduced. The canal geometry and the obstacle is the same as that in figure 1, the bead’s mean radius is also Rm = 2.466 mm and the granular material’s initial volume and mass are also similar. The model was designed using the numerical code called YADE [11, 12], based on DEM [13]. The algorithm during a time step is the following. From the previous step, the positions and velocities for every element are known. Consequently, one can deduce overlaps for every contact, from which the resultant forces and moments can be computed from the force–displacement law. Then, solving Newton’s second law of motion for every element, velocities and positions are updated and the cycle repeated on the next time step. 3.1. Force–displacement law When an overlap occurs between two discrete elements A and B, the force–displacement law computation gives a relation between the relative displacement and velocity, and the local contact force applied to both elements. This contact force is basically broken down (figure 4(b)) as the sum of a component (here Fn ) in the normal direction, and a component (here Fs ) in the tangential direction. Fn is computed from the Walton model [14], shown in figure 4(a). The normal force depends on the sign of (VnA − VnB ), with VnA and VnB the respective normal velocities of A and B, but also on the contact history. During contact, an element first undergoes an elastic response1 , computed from Fn = kn1 · δn , (3) 1 Even if a Hertzian spring were to be more correct for elastic contact, a linear spring is used. It is simpler and leads to faster computations; moreover results are found to be similar for the two springs [15]. doi:10.1088/1742-5468/2009/06/P06012 6 J. Stat. Mech. (2009) P06012 Figure 4. (a) Local normal contact force from the Walton model. (b) Two grains contact with overlap δn and local forces. Predicting the drag coefficient of a granular flow using the discrete element method Table 1. Micromechanical parameters associated with the force–displacement law. (a) Given local parameters. (b) Contact parameters: en and ϕ are given, kn1 and ks1 are computed from the local parameters and the mean radius Rm = 2.466 mm (equations (4) and (6)). (a) Grain Obstacle Wall E (N m−2 ) ν 107 0.21 107 0.34 107 · (1/26) 0.31 Contact grain/ . . . Grain Obstacle Wall −1 en 25 000 0.5 25 000 0.5 964 0.5 ks1 (N m−1 ) ϕ (deg) 5250 30 6490 18 566 19 (b) Normal contact Tangential contact kn1 (N m ) where δn is the normal overlap, and kn1 the first normal contact stiffness defined as kn1 = 2 · knA knB , knA + knB (4) where knA and knB are the local stiffnesses associated with elements A and B. Moreover, knA = E A · RA with RA the radius of the element A and E A a local stiffness parameter. A similar definition could be given for knB . Actually, kn1 is the harmonic average of the two local stiffnesses. In figure 4(a), Fnmax depends on kn1 , the beads’ mass and the relative normal velocity at the beginning of the contact. When Fn reaches Fnmax , the force–displacement curve enters the hysteretic path controlled by a harder coefficient kn2 , computed from kn1 and en , the normal restitution coefficient (given in table 1), which is the ratio between initial and final relative normal velocities in the contact history. From this definition, one can prove from the equations of motion that [14]  kn1 , (5) en = kn2 from which kn2 is deduced. The tangential behaviour follows the regularized Coulomb friction law (figure 5). The tangential stiffness is defined as ksi = 2 · ν A knA · ν B knB × ξ, ν A knA + ν B knB (6) where ν A and ν B are the Poisson coefficients of associated real materials, and ξ a numerical coefficient depending on the contact response. When i = 1, the contact is in the elastic phase so ξ = 1. When i = 2, the contact is in the hysteretic phase, and ξ = kn2 /kn1 to keep the same ratio between normal and tangential stiffnesses. Then the tangent force is computed from Fs = ksi · δs , where δs is the tangential relative incremental displacement. doi:10.1088/1742-5468/2009/06/P06012 7 J. Stat. Mech. (2009) P06012 Element Predicting the drag coefficient of a granular flow using the discrete element method The maximum tangent force Fs,max is defined as Fs,max = Fn · tan(ϕ). (7) Following Coulomb’s friction law, the tangential contact force is in the end the minimum value between Fs,max and Fs . 3.2. Micromechanical parameters The final values of the contact micromechanical parameters, like the values of the local parameters E and ν, for the materials are gathered in table 1. For the normal contact stiffness, the computed value is kn1 = O(105 ), which is representative of a smooth material, in comparison to glass. This is due to the time step—the smaller it is, the higher the stiffness is—which prohibits any large-scale study. A few studies have already shown agreement between simulations and experiments [16, 17] for similar bead sizes and contact stiffnesses. A large number of beads are required for simulation (135 000 beads), so we chose to decrease the stiffness. The ratio between the stiffnesses of different elements is equal to the corresponding real material’s Young modulus ratio. We performed simulations with various restitution coefficients without observing significant effects on the granular flow properties and the load, except when en is close to 0 or 1. Similar conclusions on the velocity profile were related for granular flows down an inclined plane [18]. For the obstacle loading, we can explain this non-influence with the load history curve for two elements: whether we use low or high  stiffness, the maximum normal force is different, but the impulse remains the same ( Fn (t) dt). Therefore, the average energy transmitted to the obstacle must be the same, which leads to similar averages of the load history in the two cases. The friction angle between a grain and the obstacle is ϕobstacle = 18◦ . Between a grain and one of the canal’s walls, it is ϕwall = 19◦ . These values both correspond to the results that we found in our laboratory from direct shear test measurements carried out on the real materials. The friction angle between two grains is ϕgrain = 30◦ , which is equal to the internal friction angle that we found experimentally during triaxial shear test measurements on glass bead samples with a height diameter ratio of 2. During these doi:10.1088/1742-5468/2009/06/P06012 8 J. Stat. Mech. (2009) P06012 Figure 5. Regularized Coulomb friction law for computing local tangential force. Predicting the drag coefficient of a granular flow using the discrete element method laboratory tests, samples were loaded in a quasi-static way. It would have been more correct to use dynamic friction angles, but they are complicated to measure, so we kept the static values. Moreover, we will see in the next section that the experimental and numerical results are close, so the static and dynamic friction angles must be close too. 3.3. Simulation methodology Before flow simulations, we performed once a gravity deposition of the granular material. The deposition process was over when the global volume of the sample was similar to that of the experimental sample, which corresponds to a low residual kinetic energy. 4. Model validation The numerical model was validated by comparing numerical and experimental parameters. Comparisons were made for some of the flow characteristics, lateral wall velocity profiles (top of figure 6) and thicknesses (bottom left of figure 6) along the slope, and the history of the load exerted on the obstacle (bottom right of figure 6). doi:10.1088/1742-5468/2009/06/P06012 9 J. Stat. Mech. (2009) P06012 Figure 6. Characteristics used for the numerical model validation: (top) lateral velocity profiles, grains considered: yes for the first and third, no for the second; (bottom left) central thickness along the slope; (bottom right) load history. Predicting the drag coefficient of a granular flow using the discrete element method 5. Obstacle shape numerical influence on the drag coefficient In this final section, we use the validated model to characterize the effect of various obstacle sizes and shapes on the drag coefficient. We define the upstream border of the influence zone created by the interaction between obstacle and flow as Dc . For obstacle sizes and shapes and their effect on Cd , we assume that Dc is still slightly similar. The first obstacle shape has its surface facing the flow. The second shape shows two perpendicular sides turned with a 45◦ angle. We also studied the influence of the ratio between the obstacle height ho and the flow thickness htop obtained during the steady regime, in the case of both immersed obstacles (ho ≤ htop ) and emerged obstacles (ho ≥ htop ). doi:10.1088/1742-5468/2009/06/P06012 10 J. Stat. Mech. (2009) P06012 The numerical method used to compute the flow thickness for a couple (t, Ds ) (Ds being the distance from the obstacle) is the following. First, a window in the (X, Z) plan is defined, centred on X = Ds and Z = 10 cm (middle line between the two lateral walls); this window is about 8 × 8 grains square. Thus, we consider in that window the position of the highest grain which is in contact with at least another grain, to avoid any mistake with grains in saltation above the free surface flow. Each thickness profile history can be split into three parts. We are interested in the central one, which corresponds to a steady phase where the thickness is nearly constant during three-fifths of the history. The thickness in the steady phase is used for comparisons, and called htop . The numerical and experimental htop profiles along the slope are close. A slight increase in the thickness can be observed from Ds /lo = 5 to Ds /lo = 2.5, which is to be understood as the effect of the obstacle on the granular flow. The velocity profiles are computed in the steady regime. To calculate the numerical wall profiles, only grains in the vicinity of the wall are considered. Therefore, a specific distance Dout from the lateral wall must be defined: if the grain’s centre position is beyond Dout , it is not taken into account in the mean velocity computation. We assume Dout to be four-fifths of the mean diameter. That method decreases the number of grains on which mean velocities are computed, which explains the scattering of the numerical velocity profiles. We note the effect of the basal wall on the velocity profiles for Ds /lo = 2.5 or Ds /lo = 5: they are concave and the top velocity is around 20% higher than the bottom velocity. Numerical velocities are not more than 10% overestimated in relation to the experiments. The Dout value can explain this overestimation; however, it would be nearly impossible to modify it. Indeed, we do not really know Dout for the experiments, and decreasing it would also lead to a decrease in the grain number, already low, on which the mean velocity is calculated. The numerical normal impact load against the obstacle is merely the resultant force computed from all the contacts between a grain and the obstacle surface. Note that we applied a moving average to the raw data, as in the experiments. The numerical model predicts the load history quite well, whether in the steady regime or in the transient regime. On the basis of the previous comparisons, the numerical model was estimated to be in sufficient agreement with the experiments and so could be used as a predictive numerical tool. Consequently, in the following section, we will use simulations to interpret flow parameters interior to the flow during the steady regime, for the computation of the drag coefficient. Predicting the drag coefficient of a granular flow using the discrete element method 5.1. Size of the influence zone The flow characteristics used in computing the drag coefficient, such as velocity and density, must be computed upstream of the impacted obstacle and upstream of its influence on the flow. Therefore, the first thing to do is to estimate Dc . Figure 7 shows basal lateral velocity profiles versus the dimensionless number Z/lo , for various values of Ds /lo and computed during the steady regime. We can observe that the flow starts to slow down in the central part of the curve when Ds /lo < 3.75, so we consider that Dc = 3.75 lo . 5.2. Mean velocity and solid fraction computation Figure 8 shows both the central velocity Ux and central solid fraction Φ profiles around Dc (the computation window is about 8 × 8 grains square). We note that the central velocities are about 20% higher than the lateral velocities (figure 6). To compute the central solid fraction profile, the altitude is discretized with a dy step. The solid fraction associated with an altitude i · dy, where i is an integer, is equal to the ratio between the solid volume intersecting the corresponding plane, which is a surface, and the total area of the computation window. Periodic undulations can be observed along the solid fraction curve. Actually, grains are regularly arranged so the period is equal to the grain’s mean size. We define the macroscopic density ρ of the flow as the glass density ρglass (equal to 2500 kg m−3 ) multiplied by the mean solid fraction Φ. Therefore, from equation (1), the drag coefficient is newly defined as Cd = Fd , (1/2) ρglass Φ Ux 2 A (8) where Φ and Ux  are respectively the mean velocity and the mean solid fraction (figure 8), computed from the base up to the reference height. Φ, Ux  and A are given for the distance Dc = 3.75 lo from the obstacle. For the reference surface, A equals doi:10.1088/1742-5468/2009/06/P06012 11 J. Stat. Mech. (2009) P06012 Figure 7. Basal lateral velocity profiles for various values of Ds /lo . The velocity starts to decrease when Ds /lo < 3.75. The lateral position of the obstacle has been added. Predicting the drag coefficient of a granular flow using the discrete element method Figure 9. Obstacle influence on Cd . Effect of the obstacle height ho normalized with the flow thickness htop , and the obstacle shape (shown on the right). lo · href , with href the reference height, defined as the minimum between the flow thickness htop and the obstacle height ho . 5.3. Evolution of Cd Results of the influence study are shown in figure 9. Note that the Froude number, as defined in equation (2), is Fr ≈ 5.5. For a similar Froude number, in [2], the empirical formula Cd = 10.8 Fr −1.3 found from full-scale experiments gives Cd = 1.18. The order of magnitude is therefore similar in our study, where Cd is in the range 1.5–2.5. The two curves are qualitatively close; Cd computed with the frontal obstacle is on average 30% higher than with the 45◦ obstacle. The interaction between the obstacle and the granular flow leads to the formation of an influence zone upstream from the obstacle, where the grain kinetic energy is as lowered as grains are closer to the obstacle. Actually, grains are more easily deviated sidewards in the case of the 45◦ obstacle, so their energy is less transmitted to the obstacle, which lowers the force that it is subjected to. doi:10.1088/1742-5468/2009/06/P06012 12 J. Stat. Mech. (2009) P06012 Figure 8. Velocity profile (left) and solid fraction profile (right) along the altitude, computed at the position Dc . Predicting the drag coefficient of a granular flow using the discrete element method 6. Conclusions and perspectives We developed a numerical model in order to predict the drag coefficient from internal characteristics of a dry granular flow. The model, based on the discrete element method, was validated with small-scale experiments on granular flow impacting a frontal obstacle, conducted in a laboratory canal. After validation, from comparisons between numerical and experimental flow characteristics and load history, we considered the model as a predictive tool. The model can take into account internal characteristics in the computation of Cd , which is generally not possible during small-scale or in situ experiments. Then, modifying the shape and the size of the numerical obstacle, we carried out a parametric study on the geometrical influence of the obstacle on the drag coefficient. Even if they are not exactly similar, because of different measurement methods and types of study, the results are in good agreement with the order of magnitude reported in the literature. The evolution of the drag coefficient for obstacle’s height between 0.5 and 3.5 times the flow thickness is similar whatever the obstacle shape. The drag coefficient is 30% higher in the case of a frontal obstacle. Moreover, both curves may be split into two parts: either the obstacle is immersed within the flow, or it is emergent out of the flow. In the first situation, the drag coefficient increases as the obstacle height decreases. In the second situation, the drag coefficient increases in the range 1–2 and becomes constant beyond this range. Those behaviours have been explained with global considerations on the grain kinetics, but in the future it would be worthwhile to explain them through the analysis of the contact network evolution around the obstacle. The drag coefficient associated with the impact forces of natural snow avalanches is influenced by a number of parameters, other than the shape effect analysed herein. The numerical tool that we developed may be used to predict the effect of the Froude number on the drag coefficient. However, our study investigated dry and somewhat ideal granular material, with simpler geometry and mechanical properties than for snow involved in real avalanches. Therefore, it would be interesting to modify the force–displacement law by implementing a cohesion component, to get closer to the material properties of flowing snow. Acknowledgments We would like to acknowledge the LOCIE laboratory where we carried out the experiments, both the ANR (National Research Agency) project OPALE, and the ‘Cluster environment’, in the Rhônes-Alpes region, for their financial support. doi:10.1088/1742-5468/2009/06/P06012 13 J. Stat. Mech. (2009) P06012 Whatever the obstacle’s shape, the curve may be split in two parts. The first one where Fr ≤ 1 for the immersed obstacle, the second one where Fr ≥ 1 for emerged ones. The left part of each curve increases as ho /htop decreases. It corresponds to an increase of the thickness of the layer above the obstacle, which leads to a lower possibility for the grains to flow over it, so they are more canalized in the direction of the obstacle, leading to an increase in its loading. The right part of each curve strongly increases as ho /htop increases in the range 1–2. For obstacle heights exceeding about two flow thicknesses, the drag coefficient no longer varies. This is due to the impossibility for grains to reach the obstacle surface above two flow thicknesses, because of their previous sidewards deviation caused by the obstacle. For the tallest obstacle, Cd is about 30% more than that for those which have the same height as the flow thickness, as also reported by [3]. Predicting the drag coefficient of a granular flow using the discrete element method References doi:10.1088/1742-5468/2009/06/P06012 14 J. Stat. Mech. (2009) P06012 [1] Berthet-Rambaud P, Structures rigides soumises aux avalanches et chutes de blocs: modélisation du comportement mécanique et caractérisation de l’interaction ‘phénomène-ouvrage’ , 2004 PhD Thesis Université Joseph Fourier [2] Thibert E, Baroudi D, Limam A and Berthet-Rambaud P, Avalanche impact pressure of an instrumented structure, 2008 Cold Reg. Sci. Technol. 54 206 [3] Hauksson S, Pagliardi M, Barbolini M and Jóhanesson T, Laboratory measurements of impact forces of supercritical granular flow against mast-like obstacles, 2007 Cold Reg. Sci. Technol. 49 54 [4] Salm B, Burkard A and Gubler H U, Calcul des avalanches: une méthode pour le praticien avec des exemples, 1990 Traduction Ancey C, Rapport du SLF n47 Davos (Switzerland) [5] Mellor M, Avalanches, 1968 Cold Regions Science and Engineering part III, section A3d (Hanover, NH: US Army CRREL) [6] Schaerer P A, Observations of avalanche impact pressures, 1973 Advances in North American Avalanche Technology, Symp. USDA Forest Service General Technical Report RM-3 pp 51–4 [7] Schaerer P A and Salway A A, Seismic and impact pressure monitoring of flowing avalanches, 1980 J. Glaciol. 26 179 [8] Sovilla B, Schaer M, Kern M and Bartelt P, Impact pressures and flow regimes in dense snow avalanches observed at the Vallée de la Sionne test site, 2008 J. Geophys. Res. 113 F01010 [9] Chehata D, Zenit R and Wassgren C R, Dense granular flow around an immersed cylinder , 2003 Phys. Fluids 15 1622 [10] Albert R, Pfeifer M A, Barabàsi A and Schiffer P, Slow drag in a granular medium, 1999 Phys. Rev. Lett. 82 205 [11] Kozicki J and Donzé F V, A new open-source software developed for numerical simulations using discrete modeling methods, 2008 Comput. Methods Appl. Mech. Eng. 197 4429 [12] YADE Open Source Discrete Element Code, 2004 http://yade.wikia.com/wiki/Yade [13] Cundall P A and Strack O D L, A discrete numerical model for granular assemblies, 1979 Géotechnique 29 47 [14] Walton O R and Braun R L, Viscosity, granular-temperature, and stress calculations for shearing assemblies of inelastic, frictional disks, 1986 J. Rheol. 30 949 [15] Silbert L E, Ertas D, Grest S, Halsey T C, Levine D and Plimpton S J, Granular flow down an inclined plane: Bagnold scaling and rheology, 2001 Phys. Rev. E 64 051302 [16] Hanes D M and Walton O R, Simulations and physical measurements of glass spheres flowing down a bumpy incline, 2000 Powder Technol. 109 133 [17] Taberlet N, Écoulements gravitaires de matériaux granulaires, 2005 PhD Thesis Université de Rennes 1 [18] Midi G D R, On dense granular flows, 2004 Eur. Phys. J. E 14 341