



































060210-F538-FAP-521689-Academic Journal of Engineering and Technology Science.docx


    

 Academic Journal of Science, Engineering and Technology 
Vol.1, Issue 1; March -2023; 

https://topjournals.org/index.php/ajset; mail: topacademicjournals@gmail.com 
 

 

 

  

15| A c a d e m i c  J o u r n a l  o f  S c i e n c e ,  E n g i n e e r i n g  a n d  T e c h n o l o g y |  

https://topjournals.org/index.php/ajset 

 

Applications of Discrete Unit Method in Simulation Analysis of Loader 

Shovelling Mechanisms 

 

 
Xiaoyuan Zhang*, Jijiang Zhao, Songsheng Duan & Zijiang Wu 
School of Intelligent Engineering, Jinzhong College of Information, Jinzhong, 030800, China 

  

Abstract: Loader shovelling mechanisms are used in various industrial applications, and the resistance 

encountered during the shovel loading process greatly affects the loading efficiency. Shovel loading resistance 

is the reaction force of the material on the bucket during the shovel loading process and contains two parts: 

the resistance of the bucket insertion stage and the resistance of the lifting stage. Therefore, it is essential to 

conduct research on shovel loading resistance to improve the performance of the loader shovelling mechanism. 

This paper presents a simulation analysis of the loader shovelling mechanism using the discrete element theory 

and EDEM software. A simulation model of the loading process is established, and the forward and reverse 

shovel loading processes are simulated. The mechanical characteristics of the shovelling operation are 

analysed visually, and the obstructing effect of the material on the bucket is clarified. The shovelling resistance 

of the key parts of the bucket at different stages of the shovelling process is specifically analysed, and the 

parts where the peak shovelling resistance is located and the stages where it is located are identified. The 

influence of the shovel angle on the resistance is also analysed. The results of the simulation analysis provide 

insight into the shovel loading resistance and the performance of the loader shovelling mechanism. The study 

demonstrates that the discrete element method can be effectively used in the simulation analysis of loader 

shovelling mechanisms, providing a useful tool for the design and optimisation of the loader shovelling 

mechanism. 

In conclusion, this research sheds light on the important mechanical characteristics of the shovel loading 

process and provides a comprehensive analysis of the shovelling resistance encountered by the loader 

shovelling mechanism. The results of this study can be applied to improve the design and performance of the 

loader shovelling mechanism, which can have a positive impact on various industrial applications that rely on 

this technology. 

 

Keywords: Shovel loading resistance, discrete element theory, EDEM software, simulation model, loading 

efficiency, mechanical characteristics, obstructing effect, peak shovelling resistance, shovel angle. 

 

INTRODUCTION 

The discrete unit method is to view the medium as a set of discrete independent moving units and to build a 

mathematical model through the properties of the discrete body, treating the object of analysis as a discrete 

particle, which corresponds to the properties of the discrete body itself. EDEM software is based on the 

discrete unit method and is widely used in many fields to calculate simulation processes quickly and 

efficiently. The software allows detailed analysis of the simulation results, such as graph types, particle 

tracking, transient analysis, etc. Yang study loader shovelling mechanisms. Includes the shovelling operation 

process and the theoretical basis[1]. Based on the discrete element principle, Pang Lizhi studied the bucket 

wheel pick-up process. Taking the pick-up machine of a power plant as the research object, he applied EDEM 

software to simulate the horizontal pick-up process[2]. Wang proposed the use of RecurDyn to construct the 

loader model, and the use of EDEM to construct the material model, and carried out the coupling simulation 

mailto:topacademicjournals@gmail.com
https://kns.cnki.net/kcms2/author/detail?v=zHhSLvPVAV77oq0bSntXyki_FPuey0ORy1srzJf6gzZaWsCnovv0zYn5hf0vSYInzePdhwjQi_Y0NHju_goB74gLAFldoNvCzKQsqxlY5fY=&uniplatform=NZKPT
https://kns.cnki.net/kcms2/author/detail?v=zHhSLvPVAV5Q_Fe6vsc70EILCiZ6bRgrQAsHumDE22TrzZ6vJOtpjwEdpPuYciNtFVyTu30dHkQLXtynTeofKItgiWTNa7LXgtnAvHZsnHE=&uniplatform=NZKPT


Xiaoyuan Zhang*, Jijiang Zhao, Songsheng Duan, Zijiang Wu 
Applications of Discrete Unit Method in Simulation Analysis of Loader Shovelling Mechanisms 

 

16| A c a d e m i c  J o u r n a l  o f  S c i e n c e ,  E n g i n e e r i n g  a n d  T e c h n o l o g y |  

https://topjournals.org/index.php/ajset 

 

analysis of RecurDyn-EDEM to study the coupling effect between the loader working device and the material 

under different trajectories[3]. Zheng according to the collapsed material working conditions occurring in the 

actual operation of the material extractor, the bucket wheel excavation EDEM simulation model was 

established by studying the bucket wheel three-dimensional model, combined with the discrete unit method; 

the time course curves of the excavation resistance of multiple buckets in the material collapse process, under 

different deep burial conditions and in different directions were calculated, and the calculated values of the 

bucket wheel excavation resistance in the design specification were compared to assess the hazards of the 

maximum excavation resistance under different deep burial conditions[4]. Li through data processing and 

setting of simulation parameters, a shovel loading model was established in EDEM software for shovel 

excavation conditions, and the accuracy of the model was verified by experimental and simulation methods; 

afterwards, the shovelloading process was simulated, and the resistance of different parts of the bucket at 

different stages of shovel loading was analysed to find out where the peak of shovel loading resistance was 

located[5].The characteristics of the bulk material and the loading conditions of the loader were analyzed by 

Yu, the characteristics of the rock material were introduced, the calculation methods for calculating the 

operating resistance were summarized, and the factors influencing the size of the loading resistance were 

classified in the shovelling process, and the main factors were grasped for analysis [6].  

Bucket model  

The bucket is an important actuator for loading, transporting and discharging materials in the loader's working 

equipment. It is usually made of the front edge (sometimes equipped with bucket teeth), the bottom of the 

bucket, the circular bucket wall, the side edge, the side wall and the rear baffle welded together, and the shape 

in the transverse direction of the body remains basically unchanged, so the geometry of the bucket is 

determined by the longitudinal section size. At present, the bucket radius of gyration R is usually used as the 

basic parameter for the calculation of other parameters in the design, with the following formula.  

R=    (1)  

R: bucket radius of rotation/m; Vs : bucket flat capacity/m3 ; B0 : bucket internal measured width/m; λg: bucket 

bottom length coefficient; λz: back wall length coefficient; λk: baffle height coefficient; λr : radius of circle 

coefficient; γ: opening angle; γ1 : angle between baffle and back wall  

Bulk materials consist of dispersed particles with a bulk substance that is not normally found in solids, liquids 

and gases, and whose garments of motion follow Newton's second law [7]. Under the action of internal forces, 

the particles of an object device undergo some kind of flow and take on the characteristics of a liquid, 

eventually forming a particle flow. In our daily life there are mainly gravel, sand, coal and grains, of which 

gravel, sand and coal are of most interest. The general rock materials are granite, limestone, sandstone, shale, 

etc. The characteristic features are different and Table 1 shows the relevant physical properties of crushed 

rock [8].  

Table 1: Physical properties of aggregates  

Properties  Numerical 

values  

Properties  Numerical 

values  

Modulus of 

elasticity E  

1.5 x 108 N/m2 Resistance 

factor  

0.20  

Density  1.9 x 103 kg/m3 Coefficient 

ostatic friction  

0.90  

Friction angle  32.5°  Rolling friction 

coefficient  

0.66  

Poisson's ratio  0.35    

[ ] 
 
 
 

 
 
 

 
 

 
 
 

 
− − − + ) 

180 
0.5(1 

2 
cos 0.5 cos 

2 
1 k 0 

γ 
π 

γ 
λλγλ r Z 

S 

B 

V 

https://kns.cnki.net/kcms2/author/detail?v=zHhSLvPVAV6wX5ktKNIV0SP4DXVXKr11JFybKTsCBdfAdtyLqcucg4TFZHUBh8-TQk6T7VCasOpihblCRae-IKPb8_BMq3wzxn3-II9e960=&uniplatform=NZKPT
https://kns.cnki.net/kcms2/author/detail?v=zHhSLvPVAV6wX5ktKNIV0SP4DXVXKr11JFybKTsCBdfAdtyLqcucg4TFZHUBh8-TQk6T7VCasOpihblCRae-IKPb8_BMq3wzxn3-II9e960=&uniplatform=NZKPT
https://kns.cnki.net/kcms2/author/detail?v=zHhSLvPVAV65l8IcOd58AiciPVz23SSVkIY4rbjF3tUT-GZjGDMTC15AZKfAIcqdKTi_Rj2bvm65iOaTw30bce94WDnC5mXwGf-EdUGCEn0=&uniplatform=NZKPT


Xiaoyuan Zhang*, Jijiang Zhao, Songsheng Duan, Zijiang Wu 
Applications of Discrete Unit Method in Simulation Analysis of Loader Shovelling Mechanisms 

 

17| A c a d e m i c  J o u r n a l  o f  S c i e n c e ,  E n g i n e e r i n g  a n d  T e c h n o l o g y |  

https://topjournals.org/index.php/ajset 

 

Non-adhesive spherical particle contact forces (Hertz theory)  

The contact model is an important basis for modelling the discrete element method, its description of the 

contact behaviour between elements and the analytical calculations directly determine the magnitude of the 

forces and moments applied to the particles. Therefore when using the discrete element method for different 

objects, the contact models differ and the results vary, but all simulation models must include at least one 

basic particle-to-particle and particle-to-geometry boundary contact model [9-12]. In this paper, in order to 

simplify the simulation, the default contact model Hertz theory model set in the EDEM software is used 

uniformly, which has an efficient and accurate computational performance.  

Hertz contact theory assumes that the surfaces of particles in contact with each other are smooth and 

homogeneous, that the contact surface is small compared to the particle surface, that only elastic deformation 

occurs at the contact surface, and that the contact force is perpendicular to the contact surface [13-14]. The 

Hertz contact theory is the theoretical basis of the problem and is applicable to the elastic contact of curved 

bodies such as spheres, columns and ellipsoids, and even to the contact of micro-convex bodies between 

contact surfaces. As shown in Fig.1, two spherical particles of radii R1 and R2 are in elastic contact, and the 

normal overlap α is 

α = R1 + R2 -|r1 -r2 |>0                                (2)  

The radii of particle 1 and particle 2 are R1 , R2; r1 and r2 are the spherical position vectors of the two particles 

respectively.  

The dotted line in Figure 1 shows the location of the particle surface when no deformation is considered, and 

the contact surface between the particles is circular, then the contact surface radius a is:  

a= aR*   (3)  

The inter-particle normal forces N are:  

 1 3 

N=E R*( *) 2α 2                               (4)  

R* and E* are the effective particle radius and effective modulus of elasticity respectively.  

 1 1 1 

 = +    (5)  

 R* R R1 2 

 1 1−ν12 1−ν22 

 = +    (6)  

 E* E1 E2 

E1,v1,E2,v2 are the modulus of elasticity and Poisson's ratio of particle 1 and particle 2, respectively.  

When the increment of overlap between two contacting particles is Δα, the increment of normal force ΔN is 

calculated from (4) or (5).  

∆ =N 2 αR E* *∆ =α 2a *E∆α   (7)  

Analysis of forward shovelling resistance  

Forces on the bucket in the X, Y and Z directions  

The forces on the bucket in the X, Y and Z directions are first analysed. In the simulation, a certain number 

of gravel particles are loaded into the bucket. The bucket is first shovelled into the rubble pile and then the 

bucket filled with rubble leaves the pile. The software offers linear translational rotation, sinusoidal 

translational rotation and convey or translational rotation. Therefore, here we use linear translation and 

rotation. First, the bucket is shovelled parallel into the pile to a certain depth, then the bucket is turned over 

and lifted up. Figures 2, 3, 4 and 5 show the variation of forces in the X,Y and Z directions versus time. Where 

the full shovelling process is 15s. The horizontal axis represents the time of shovelling and the vertical axis 

represents the amount of shovelling resistance applied in each direction.  



Xiaoyuan Zhang*, Jijiang Zhao, Songsheng Duan, Zijiang Wu 
Applications of Discrete Unit Method in Simulation Analysis of Loader Shovelling Mechanisms 

 

18| A c a d e m i c  J o u r n a l  o f  S c i e n c e ,  E n g i n e e r i n g  a n d  T e c h n o l o g y |  

https://topjournals.org/index.php/ajset 

 

 
Figure 1: Particle contact deformation diagram   Figure 2: Shovel loading simulation  

Analysis of Figure 3 shows that the resistance of the bucket in the X-axis direction exceeds 9KN when the 

bucket is just shovelling into the rubble pile, there are small fluctuations in the resistance during the process 

of shovelling into the rubble, but when the bucket is filled with rubble the resistance rises rapidly to 12KN 

when the bucket is turned over and lifted off the pile, then the maximum partial force falls rapidly to the lowest 

resistance corresponding to the bucket shovelling out of the rubble pile. The X-directional force then does not 

change much. The peak force is generated in approximately 1.5s.  

 
 
Figure 3: X-directional component of shovel loading resistance  

The peak force is approximately 1.5 s. The corresponding force in the Z-axis in Figure 5 is not very large 

compared to the X and Y-axis. Y-axis is not very large. The rest of the time is very smooth. The z-directional 

force is much smaller than the x- and y-directions.  

 
Figure 4: Y-direction of the shovel load resistance  

 

2 a  

a  

2                5              8  
Time/s  

12  

 

 

11  

 

 

10  

 

 

9  

x
 - a

x
is

  f
o

rc
e/

K
N

 
 

 

2                5               8  

Time/s  

-8  

 

 
-9  

 

 
-10  

 

 
-11  

y
 - a

x
is

  f
o

rc
e/

K
N

 
 

 



Xiaoyuan Zhang*, Jijiang Zhao, Songsheng Duan, Zijiang Wu 
Applications of Discrete Unit Method in Simulation Analysis of Loader Shovelling Mechanisms 

 

19| A c a d e m i c  J o u r n a l  o f  S c i e n c e ,  E n g i n e e r i n g  a n d  T e c h n o l o g y |  

https://topjournals.org/index.php/ajset 

 

 
Figure 5: Z-directional component of shovel loading resistance  

Force analysis of the bucket base plate  

As the forces exerted on the base plate are the greatest during the shovelling process, the next discussion will 

focus on the base plate and the left and right plates of the bucket will not be analysed for the time being. 

Simplifying the actual shovellingprocess of the loader, the shovel-in phase lasts 7 seconds and the bucket 

speed is 0.2 m/s; the effective working phase lasts 4 seconds and the angular speed of the bucket is 0.15 arc/s. 

The forces exerted by the material on the bucket floor during the horizontal shovel-in phase are shown in 

Figure 6. The forces exerted by the material on the bucket during the bucket reversal lift are shown in Figure 

7.  

 
 
Figure 6: Bucket forces in the horizontal shovel-in phase  

 
 
Figure 7: Bucket stress diagram in lifting phase  

The analysis of the shovel loading resistance during the horizontal insertion phase showed that during the 

horizontal insertion phase the total force at the bottom first increased, then decreased and finally stabilised, 

3   5  7  

Time/s  

3.5  

 

3  

 

2.5  

 

2  

F
o

rc
e/

K
N

 
 

 

7               9              11  
Time/s  

2  

 

1.5  

 

1  

 

0.5  

F
o

rc
e/

K
N

  
 



Xiaoyuan Zhang*, Jijiang Zhao, Songsheng Duan, Zijiang Wu 
Applications of Discrete Unit Method in Simulation Analysis of Loader Shovelling Mechanisms 

 

20| A c a d e m i c  J o u r n a l  o f  S c i e n c e ,  E n g i n e e r i n g  a n d  T e c h n o l o g y |  

https://topjournals.org/index.php/ajset 

 

with local fluctuations during which it showed a secondary peak. The maximum force was 3.5 KN. In the 

subsequent stages, the force gradually stabilised when the amount of material loaded reached a threshold. As 

the bucket rotates upwards, the forces previously concentrated at the bottom of the bucket are gradually 

transferred to the circular wall of the bucket during this phase, so that the forces acting on the bottom of the 

bucket are significantly reduced during this phase.  

Results of reverse shovel analysis  

Forces on the bucket in the X, Y and Z directions  

Figures 8, 9 and 10 show the relationship between the total force applied to the bucket during excavation in 

the X, Y and Z directions as a function of time. Where the full shovelling process is 14s. where the horizontal 

axis indicates the digging time and the vertical axis indicates the magnitude of the digging resistance applied 

in each direction.  

 
 
Figure 8: X-directional component of the digging resistance  

The relationship between the X-directional force and time can be seen in Figure 8. The resistance of the bucket 

in the X-direction was just over 7KN when the bucket was just shovelling into the rubble pile, and its X-

directional force reached its maximum value of 9.8KN around the 8th s. Between 5.5s and 8s, the bucket force 

had a more obvious fluctuation. As the time increases, its digging resistance gradually decreases, and after 9s 

its resistance drops significantly.  

 
Figure 9: The Y-directional component of the digging resistance  

3   9  6  

Time/s  

10  
 
 

9  
 
 

8  
 
 

7  

x
 - a

x
is

  f
o

rc
e/

K
N

  
 

3               6             9  

Time/s  

-10  

 

20  

 

15  

 

-25  

Y
 - a

x
is

 f
o

rc
e/

K
N

  



Xiaoyuan Zhang*, Jijiang Zhao, Songsheng Duan, Zijiang Wu 
Applications of Discrete Unit Method in Simulation Analysis of Loader Shovelling Mechanisms 

 

21| A c a d e m i c  J o u r n a l  o f  S c i e n c e ,  E n g i n e e r i n g  a n d  T e c h n o l o g y |  

https://topjournals.org/index.php/ajset 

 

 
 
Figure 10: Z-directional component of the digging resistance  

From the results of the Y-axis force splitting simulated in Figure 9, it can be seen that the bucket's force 

splitting decreases gradually from 3 to 7.5s, reaching a minimum of -25KN, and gradually increases after 7.5s. 

After 6s, the force varies considerably, especially between 7 and 7.5s, when the trend of force decrease is 

extremely obvious.  

The results of the Z-directional forces simulated in Figure 10 show that the bucket's digging resistance 

gradually tends to increase from 3s onwards, reaching a peak of 32KN at around 5s, and then showing a 

decline to a minimum of 8KN at the initial resistance.  

Force analysis of the bucket base plate  

The following analysis of the horizontal insertion phase of the excavation resistance analysis, in edem derived 

from the bucket bottom part in the horizontal insertion phase of the force diagram as Figure 11, horizontal 

insertion phase, the bucket bottom overall force first increase and then decrease, local fluctuations. In the later 

stages, the forces are gradually stabilised as the amount of excavated material has reached its limit. As the 

bucket is turned upwards, the bucket floor, which previously bore the main part, is gradually shifted to the 

circular wall at this stage, so that the forces on the floor at this stage are on a decreasing trend, as shown in 

Figure 12.  
 

 

 
Figure 11: Bucket forces in the horizontal insertion phase  

3              6             9  

Time/s  

32  

 

 

24  

 

 

16  

 

 

8  

Z
 - a

x
is

  f
o

rc
e/

K
N

  
 

10   12  14  

Time/ 

12  
 
 

11  
 
 

10  
 
 

9  

F
o

rc
e/

K
N

  
 



Xiaoyuan Zhang*, Jijiang Zhao, Songsheng Duan, Zijiang Wu 
Applications of Discrete Unit Method in Simulation Analysis of Loader Shovelling Mechanisms 

 

22| A c a d e m i c  J o u r n a l  o f  S c i e n c e ,  E n g i n e e r i n g  a n d  T e c h n o l o g y |  

https://topjournals.org/index.php/ajset 

 

 
Figure 12: Bucket forces during the pick-up phase  

Effect of shovel entry angle on resistance  

When shovelling horizontally, the bucket is shovelled into the rubble pile at a natural horizontal ground level, 

but due to the influence of the front teeth of the bucket, the bottom plate of the bucket is tilted at a certain 

angle to the horizontal surface and cannot achieve true horizontal; while shovelling into the rubble pile at a 

certain angle, the bottom plate of the bucket must be at a certain angle to the horizontal surface. When the 

shovel entry angle is less than 11°, the simulation result is close to the calculated value. In particular, when 

the shovel entry angle is between 5° and 10°, the two almost coincide. Therefore, we analytically derived the 

process of shovelling the bucket into the rubble pile, where 7° and 9° were targeted, in order to obtain the best 

insertion angle. As shown in Figure 13, it is concluded that the shovel entry resistance increases with 

increasing shovel entry depth, and with increasing shovel entry angle. When the shovel penetration depth is 

0.3 m, the effect of the shovel penetration angle on the shovel penetration resistance is relatively small; when 

the shovel penetration depth exceeds 0.5 m, there is a significant increase in the shovel penetration resistance.  

 
Figure 13: Resistance variation curve  

Conclusion  

Based on the discrete element theory, this paper applies EDEM software to establish a simulation model of 

the loading process and realise the simulation of the forward and reverse shovelling operation process. 

Through the simulation, the shovel resistance of different parts of the bucket at different stages of the shovel 

loading process is analysed, and the parts where the peak shovel resistance is located and the stages where it 

is located are identified. The resistance of the bucket in the X, Y and Z directions was analysed in the forward 

and reverse shovelling process. According to the force analysis of the bucket, the force on the bucket floor is 

the largest, after which the force analysis of the bucket floor during forward and reverse shovelling was carried 

out to obtain the trend of force changes. Finally, the influence of the shovel entry angle on the resistance was 

analysed. It can be seen that when the shovel entry angle is 9°, the resistance increases with the depth of shovel 

entry.  

3              6             9  

Time/ 

24  

 

 
18  

 

 
12  

 

 
6  

F
o

rc
e/

K
N

  
 

 

R
es

is
ta

n
ce

/K
N

  
 

2  

 

1.5  

 

1  

 

0.5  

Bucket  penetration  

0.2              0.4             0.6  

Shovel-in angle 7°  
Shovel-in angle 9°  

 



Xiaoyuan Zhang*, Jijiang Zhao, Songsheng Duan, Zijiang Wu 
Applications of Discrete Unit Method in Simulation Analysis of Loader Shovelling Mechanisms 

 

23| A c a d e m i c  J o u r n a l  o f  S c i e n c e ,  E n g i n e e r i n g  a n d  T e c h n o l o g y |  

https://topjournals.org/index.php/ajset 

 

Acknowledgements  

This paper was supported by①“2022 Science and Technology Innovation Project for Higher Education 

Institutions in Shanxi Province (2022L659)”②“2022 Innovation and Entrepreneurship  

Project for Students in Shanxi Province(20221640)”  

References  

[1]. Yang Zikang. Research on resistance reduction shoveling strategy of loader based on EDEM [D]. Jilin 

University, 2022.  

[2]. Pang, L. C.. Research on the bucket-material interaction mechanism of bucket wheel stacker bucket 

based on discrete element method [D]. Jilin University, 2022.  

[3]. Wang Shaojie, Yin Yue, Yu Shengfeng, Hou Liang. Simulation analysis of loader coupling dynamics 

based on RecurDyn-EDEM [J]. Machine Design,2021,38(11):1-6.  

[4]. Zheng Pei, Song Wenxi, Hu Xiong. EDEM simulation of excavation resistance of bucket wheel stacker 

reclaimer under collapsed material condition [J]. Sintered pellets, 2019,44(06):50-54.  

[5]. Li R, Xu Wubin, Li Bing, Yang Xu. Research on discrete unit method for bucket shovel resistance of 

loader [J]. Journal of Guangxi University of Science and Technology,2017,28(03):77-82.  

[6]. Yu, Xuebang. Simulation study of loader bucket operation resistance based on EDEM [D]. Guangxi 

University of Science and Technology, 2018.  

[7]. Seife C. Can the laws of physics be unified? Science,2005,309(5731): 82-82.  

[8]. Yan Bo, Zhan Kai, Guo Xin, Li Hengtong, Shi Xiaojie. Simulation study of underground scraper 

shoveling process based on EDEM [J]. Nonferrous Metals (Mining part),2019, 71(06):74-77.  

[9]. Sun QC, Wang GQ. Introduction to the Mechanics of Particulate Matter. Beijing:Science Press, 2009.  

[10]. Sun Qicheng, Hou Meiying, Jin Feng. Physics and Mechanics of Particulate Matter. Beijing: Science 

Press, 2011.  

[11]. Jaeger HM, Nagel SR, Behringer RP. Granular solids, liquids, and gases. Reu Modern Phgs, 1996, 

68(4): 1259-12734.  

[12]. Kadanoff LP. Built upon sand: theoretical ideas inspiredby granular flows. Rev Modern Phys, 1999, 

71(1): 435-4445.  

[13]. De Gennes PG. Granular matter: a tentative view.RevModern Phys, 1999, 71(2): S374-S382. 

[14]Forterre Y, Pouliquen O. Flows of dense granular media.Annual Review Flauid Mechanics, 2008, 

40: 1-26.   


