E-proceedings of the 36th IAHR World Congress
28 June – 3 July, 2015, The Hague, the Netherlands

EFFECTS OF PARTICLE SHAPES AND SIZES ON VERTICAL SORTINGSIN NUMERICAL EXPERIMENTS OF SLOPE FAILURE IN THE WATER
HAZUMU TADOKORO
(1)

(1)

& SHOJI FUKUOKA(2)

Graduate School of Science and Engineering, Chuo University, Kasuga 1-13-27, Bunkyo-ku, Tokyo 112-8551, Japan,
e-mail : danhazumu@civil.chuo-u.ac.jp(2)

Research and Development Initiative, Chuo University, Kasuga 1-13-27, Bunkyo-ku, Tokyo 112-8551, Japan,
e-mail : sfuku@tamacc.chuo-u.ac.jp

ABSTRACT
Gravel-bed rivers are composed of particles of various shapes and sizes. The shapes and sizes of particles in the gravelbed river have great influences on sediment transport and river bed variation. Numerical methods of sediment transportand river bed variation in gravel-bed rivers must consider effects of particles shapes properly. This paper mainly
demonstrates the effects of particle shapes on vertical sorting by comparing motions of spheres and gravels in numericalslope failure experiments in the water. The numerical experimental results demonstrated that the angle of repose of gravel
particles was larger than that of spheres due to engagement effects of gravel particles. Moreover, slope failures of both
upper part and lower part of sphere-particles stratum progressed gradually. On the other hand, the slope failure of lower
part of gravel particles progressed more slowly compared to sphere particles due to large engagement effects among
particles. Then upper particles of gravel fell down in a relatively long time. In both cases of spheres and gravels, large
particles tended to move parallel to the channel bed, but small particles tended to drop in the lower part of the channel.
Consequently, the vertical sorting made particles on the slope surface larger. Particles on the surface of the sphere slope
became larger than those of gravels. Moreover, the porosity of gravels in the moving layer was larger than that of spheres.
Keywords: particle shape, vertical sorting, gravel-bed rivers, volume concentration, numerical slope failure experiment
1.

INTRODUCTION

Gravel-bed rivers are composed of widely distributed particles in size. The sediment-transport process in gravel-bed rivers
during floods is explained as follows: First, small particles on the river-bed move. Then, large particles are exposed, resist
against hydrodynamic forces primary and reduce entrainment of small particles. Those process bring stable river bed(Fukuoka et al., 2007). Osada et al. (2012) developed a numerical method of bed variation analysis in gravel-bed rivers
considering resistance by exposed large particles. Parker (1986) showed that large particles tend to move in high level onriverbeds when all size particle can move, and riverbeds became stable structures against streams due to sorting ofparticles on riverbed surfaces. In this manner, sorting of particles has a large influence on sediment transport rate and bed
variation in gravel-bed rivers.
Vertical sorting of riverbed are strongly affected by not only sizes but also shapes of particles. Matsushima et al. (2009)
developed a numerical model capable of simulating motion of non-spherical particles which superposed small spheres as
particle models and computed contact forces by the distinct element method. They showed that shear strength of nonspherical particles were larger in comparison with spheres. In general, we have assumed that spheres are representativeof particles in derivation of sediment transport formulas, for the simplification of the phenomenon.
Effects of particle shapes on sediment transport and bed variation have not been clarified, because it is difficult to observe
particle motion and sediment transport rates in gravel-bed river during floods. Fukuda et al. (2012) developed a numericalmovable-bed channel which can simulate three-dimensional motions of flows and gravel particles in different shapes andsizes to investigate movement of gravel particles and sorting of the riverbed particles in streams. Fukuoka et al. (2014)
evaluated effects of particle shapes on sediment transport and sorting structures by comparing particle motions between
spheres and gravels in a numerical movable-bed channel. In this study, we conducted numerical slope failure experimentsof spheres and gravels in the water using the model of Fukuda et al. (2012). It is necessary to understand the Lagrangianparticle movement of different particles in size and shape to predict Eulerian motion of particle group and structure ofriverbed. Therefore we discuss Lagrangian particle motions from the particles slope failure experiments in the water. The
effects of particle shapes and sizes on vertical sorting were investigated from slope failure experiments using spheres andgravels. The slope failure experiment provides the comparison of average trajectories, changes of particle size distribution
on the slope surface, failure process of the slope and the porosity of particles.

1E-proceedings of the 36th IAHR World Congress,
28 June – 3 July, 2015, The Hague, the Netherlands

2.

NUMERICAL METHOD

Particle motions were simulated in the Lagrangian method as the rigid-body and fluid motions were simulated in the
Eulerian method. To take into account the effect of solid phase on liquid phase, fluid motions were simulated by thegoverning equations of multiphase flows. Particle contacts were computed by distinct element method. Fluid dynamic
forces on particles were computed directly by integrating the forces on a particle region in the multiphase flow. Particle
motions were simulated by the momentum equation and the angular momentum equations of rigid bodies.
2.1 Governing equations of fluid motion
We use one-fluid model for the solid-liquid multiphase flow to simulate fluid motions in the Eulerian method, using the
Smagorinsky model as the subgrid turbulence model:

ui
0
x i

[1]

ui
u
1 P

 u j i  gi 2   t Sijt
x j xi x j



Sij 



1  ui u j 
2  x j xi   
 t  CS 2 2Sij Sij

[2]

[3]
[4]

[5]

where ui : i-th component of averaged velocity given by the weight of the mass within a fluid cell, P : sum of the pressure
and isotropic component of SGS stress,  : density,  : dynamic viscosity, g i : gravitational acceleration, t : SGS
turbulent viscosity, C s : Smagorinsky constant,  : computational grid size.
Physical property  (density  , dynamic viscosity  ) and averaged velocity ui were calculated as volume-averaged
values and mass-averaged values, respectively, as follows:  s  1   f,

 f  fl  1  f gui  susi  1    f u fi  

[6]
[7]

where f : fluid volume fraction in fluid computation cell,  : solid volume fraction in fluid computation cell, suffixes l , s
and g denotes liquid phase, solid phase, and gas phase, respectively. u f is fluid velocity, and u s is solid velocity.
Continuity equation and momentum equation were solved by the SMAC scheme with staggered grids.
The water-surface variation was evaluated with continuity equation of fluid volume fraction f in fluid computation cells.

f fu i0t x i2[8]E-proceedings of the 36th IAHR World Congress
28 June – 3 July, 2015, The Hague, the Netherlands

2.2

Governing equations of particle motion in streams

We use the momentum equations and the angular momentum equations of rigid bodies to simulate particle motions in the
Lagrangian method :MrG  Mg  Ff  Fc

[9]



1
ω r  I r R 1 N f  N c   ωr  I r ωr



[10]

where rG :translational acceleration vector of particles, suffix f , c : fluid force, contact force between particles,
respectively, suffix r :components in the local coordinate systems of each particles, R :translational matrix from the global
coordinate system to the local coordinate system and I :tensor of momentum inertia. After the computation of rotational
velocity in local coordinate system  r in Eq. [10], rotational velocity  in the global coordinate system is solved bycoordinate transformation. In computation of coordinate transformation from the local coordinate system of the particles to
the global coordinate system, the quaternion is used instead of translational matrix R (Usijima et al., 2008).
2.3

Computation of fluid force

Fluid forces on a particle were computed by integrating forces on a particle region in the multiphase flow.

 P


2   t Sij dF f ,i   s x j xi


[11]

 P


2   t Skl dN f ,i   s  ijk rf , j 

xl
 xk


[12]

where Ff , i : i-th component of fluid force, N f , i : i-th component of torque exerted by fluid force, rf ,i : position vector fromthe gravity center of particle to the fluid computation cells,  s : an area occupied by particles and  ijk : Levi-Civita symbol.
2.4

Computation of Contact ForceContact forces acting on particles were computed by the contact detection of each small sphere composing a gravel
particle. Contact forces and torques on centers of gravity of particles were calculated by the summation of the contact
forces computed concerning each small sphere:

Fc   Fcp , n,

N c   rcp, n  Fcp, n[13]where Fc , N c : contact force and torque acting on the center of gravity of a particle, Fcp,n :contact force acting on eachsphere, rcp, n : position vector from the gravity center of particle to the contact point.
Contact forces between particles were computed by the distinct elements method (P.A. Cundall et al., 1979). The spring
constants k and coefficients of dashpot c were calculated by Eq. [14] through [18] (Tsuji Y., 1992).
1

 4  r r  E 2  3 ekn    1 2 2  n
 9  r1  r2  1  pos  

s0 

ks
1

k n 21  pos 

[14]

[15]

3

E-proceedings of the 36th IAHR World Congress,
28 June – 3 July, 2015, The Hague, the Netherlands

cn  2 h

m1m2kn
m1  m2[16]

c s  c n s0
h

[17]

ln b

[18]

  ln b 22

where e : spring force, E : elastic modulus, pos : Poisson’s ratio, r1 , r2 : radii of two contacting spheres, m1 , and m2 :
masses of contacting spheres suffixes n and s : a component of direction from center of spheres to the contact point andtwo orthogonal directions on the tangential plane, respectively, and b : coefficient of restitution. The spring constants k
was derived from the relationship between the normal force and displacement of Hertzian contact theory (Tsuji.Y , 1992).
Dashpot coefficient h was derived from damping oscillation of two objects ( m1 , m2 ). h is related to the coefficient ofrestitution b and takes the value 0.11 corresponding to the value of the coefficient of restitution, 0.7, in this calculation.
The parameters used in the simulation are shown in table 1.

Table 1. Parameters used in simulation.

ΔｘΔｙΔｚ ： Fluid calculation cell sizeΔt : Time step for fluid calculation
ρ w : Density of waterρ s : Density of gravel
μ w : Dynamic viscosity of water
μ s : Dynamic viscosity of gravel
Δt' : Time step for DEM calculation
E : Elastic moduls
pos : Poissom's ratio
h : Coefficient of dashpots

3.

0.01
5.0×10 -4
1,000
2,650
8.9×10 -4
8.9×10 -4
5.0×10 -6
5.0×10 10
0.33
0.11

m
s
kg/m3
kg/m3
Pa・s
Pa・s
s
Pa
-

NUMERICAL SLOPE FAILURE EXPERIMENTS IN THE WATER

3.1 Experimental conditions
We carried out two numerical slope failure experiments by using spheres and gravels. Those numerical experiments used
five particle sizes in both spheres and gravels, and four particle shapes in gravels. Gravel particles were made to imitateshapes of real gravel particles by superposing several small spheres without gap space between small ones shown in
Figure 1. Particle size distributions (40mm:blue, 50mm:red, 70mm:light blue, 90mm:green and 120mm:yellow) of spheres
and gravels are shown in Figure 2. Diameters of gravels were defined as diameters of spheres having the same volume.
Three orthogonal axial lengths a, b and c (lengths of the longest, the intermediate, and the shortest axes of each particles)
are shown in table 2. The numbers of four different shapes of gravels in each size were the same. We packed particles
into the numerical channels (length: 8 m, width: 1 m and depth: 3 m) so that the length and the height of the particle groupwere 2m respectively, and width of particle group was 1m. A wall prevents the particle group from failure. Particles arepacked as follows: At first, 1 set of particles was packed into the numerical channels, the particle size distribution of which
agreed with that shown in Figure 2. After packing 1set of particles, next set of particles was packed. We continued theoperation until the last set. In packing, particles were dropped at random into the channel. To check structures of packed
particles of spheres and gravels, vertical distributions of volume concentrations (porosity) were measured, as shown in
Figure 3. Those concentrations were calculated in a section (0.6m square in the x-y plane) within every 0.01m in thevertical direction. Values of vertical distribution of volume concentrations of both spheres and gravels were the same andnearly 0.7 (porosity 0.3). Slope failures experiments were carried out under the condition of the depth of the water 3m byremoving the wall instantly.

4E-proceedings of the 36th IAHR World Congress
28 June – 3 July, 2015, The Hague, the Netherlands

SphereFigure 1. Shape of gravels particles.

Table 2. Length of a, b and c axes of gravel particles.

通過重量100分率(%)

Cumulative percentage (%)

shape No.
a: longest
b: intermidatec: shortest

sphere
1
1
1

1
1.26
0.98
0.88

2
1.29
1.06
0.81

3
1.49
0.89
0.76

4
1.36
0.99
0.78

100
80
60
40
20
0
10

100
粒径(mm)

1000

Grain Size (mm)Distance fromthe channel bottom (m)

合計
Total

D=40mm

D=90mm

2.00

D=120mm

Sphere

1.60
1.20
0.80
0.40
0.00
0.0

0.2

0.4

0.6

Volume concentrations

0.8

1.0

Distance fromthe channel bottom (m)Figure 2. Particle size distribution.

2.00

Gravel

1.60
1.20
0.80
0.40
0.00
0.0

0.2

0.4

0.6

0.8

1.0

Volume concentrations

Figure 3. Vertical distributions of volume concentrations.

5E-proceedings of the 36th IAHR World Congress,
28 June – 3 July, 2015, The Hague, the Netherlands

3.2 Experimental results
3.2.1 Angle of repose of submerged particle group
The angles of repose (φ) of submerged sphere particles and gravel particles in the center of the channel after slope failure
are shown in Figure 4. The angle of repose of sphere particles was φ =22˚ (tanφ= 0.4), and that of gravel particles was φ =31˚ (tanφ = 0.6). The experimental result demonstrated that tanφ of the submerged gravel particles was 1.5 times (0.6/0.4)
larger than that of sphere particles due to the effect of shapes of gravel particles. A difference of particle shapes result in a
noticeable difference in angle of repose. We discuss mechanism making the difference between these angles of repose
from Figure 4. This figure shows that gravel particles shared a flat surface with others, however sphere particles contactedwith others at a point. Gravels have a large surface area by its shape deviating from a sphere. Their contact surfaces were
nearly flat. For this reason, when gravels move, gravels have a tendency to rotate with surrounding gravels by sharing
contact surfaces. Then neighbouring gravels resist movements of gravel particles. We call this the engagement effect ofgravels.
The effect of the height of packed particles on the angle of repose was checked by slope failure experiments under thecondition that height, length and width of packed particles were 1m respectively. The tip of travel distances of spheres(2.5m) and gravels (2.1m) became small in comparison with the previous experiments, but the tanφ (sphere: 0.4, gravels:
0.6) were almost unchanged.

Sphereφ＝22°
tanφ＝0.4

Figure 4a. Angle of repose of spheres in the water.

Gravelφ＝31°
tanφ＝0.6

Figure 4b. Angle of repose of gravels in the water.

3.2.2 Stopping process in slope failure
Failure processes of particles are demonstrated with coordinate axes in Figure 5. When a wall was removed instantly, theslope failure of sphere particles progressed smoothly with time from the upper part to the lower part, and formed a stableslope. On the other hand, the slope failure of gravel particles is different from spheres. The lower part of gravel particles
progressed more slowly due to large engagement effects among particles compared to spheres. Then, upper particles ofgravel particles fell down. It took 9 seconds in gravel case and took 7 seconds in sphere case until stable condition of
failure process. The tip of particles stopped at first, then following particles on slope surface gradually stopped in bothspheres and gravels. The travel distance of tip of spheres (5.0m) was larger than that of gravels (4.3m). The difference ofparticle shapes cause difference in failure process and travel distance of tip particles.
Conditions of channel bottom may have an influence on angle of repose φ. Therefore, we changed bed condition fromfixed to movable bed. Particles were set with a height of 0.2m over the whole channel. Then, the tip travel distances
became small for spheres (4.8m) and gravels (4.0m), because failure particles contact with the movable bed particles, andlet them move. However, angles of repose of spheres (tanφ= 0.4) and gravels (tanφ= 0.6) hardly changed. Thus, tanφ isinvariant under different conditions between fixed channel and movable bed channel.

6E-proceedings of the 36th IAHR World Congress
28 June – 3 July, 2015, The Hague, the Netherlands

Sphere1.5s

Gravel

1.5s

3m

4m

2m

z
y
x
0m

1m

2m

3m

4m

5m

0m

1m

2m

2.5s

5m

5.5s

Tip part
0m

1m

2m

3m

4m

5m

0m

1m

2m

3m

3.5s

1m

2m

3m

4m

5m

7.0s

Tip part
0m

4m

Tip part
5m

0m

1m

2m

3m

7.0s

4m

5m

9.0s

Tip part

0m

1m

2m

3m

4m

5m

Figure 5a. Time process of sphere failure particles.

0m

1m

2m

3m

4m

5m

Figure 5b. Time process of gravel failure particles.

3.2.3 Vertical sorting in slope failures
Many large yellow and green particles in both spheres and gravels were seen on the surface after slope failures, as shownin Figure 6. Many small ones were located in lower part, as shown in Figure 5. Large particles tended to more parallel tothe channel bed, but small particles tended to drop in the lower part of the channel. Consequently, the vertical sorting
made particles of the surface larger. The longer travel distances of particles are, the more small particles fall into the gaps
made by moving particles. Sphere and gravel size distributions at the surface after the failure are shown in Figure 7.
Particle size distributions were measured at x=1-3m where particles moved actively. In the downstream section (x > 3m),
the motion of tip particle groups were greatly influenced by the channel bottom. Most of the particles at the x < 1m did not
move. Therefore, the measurement was carried out at 1m < x < 3m. The average sphere particle size (dm=81mm) on
surface after the failure was larger than initially packed one (dm=70mm). Similarly, gravel particle size (dm=75mm) after the
failure is larger than initial one (dm=70mm).

7E-proceedings of the 36th IAHR World Congress,
28 June – 3 July, 2015, The Hague, the Netherlands

SphereGravel

0m

0m
2m

2m

4m

4m

Figure 6a. Slope surface of sphere.

Figure 6b. Slope surface of gravel.

30

Sphere
20
10
0
40

50

Particle size distribution on slope surface
after slope failure

Weight percentage (%)

Weight percentage (%)

Particle size distribution
before slope failure

70

90

Diameter (mm)120

Figure 7a. Sphere size distribution after failure.

30

Gravel
20
10
0
40

50

70

90

Diameter (mm)120

Figure 7b. Gravel size distribution after failure.

Figure 8 shows time changes in the ratios of average vertical travel distances ∆z of particle group to average horizontal
travel distances ∆x of particle group in the slope failure. We chased particles of D=40mm (blue) and D=120mm (yellow)
within x=0.8-1.2m, y=-0.3-0.3m and z=1.6-2.0m during the process of the slope failure. Both spheres and gravels formed aslope generally from start of failure to approximately 4 seconds. There were the difference in ∆z/∆x between small
particles and large particles. This reason is that large particles move parallel to the channel bed, but small particles fall into
gaps made by moving particles. Time of the large difference of spheres ∆z/∆x between small particles and large particles
is longer than that of gravels.
Particle volume concentrations in the center of the channel after slope failure were investigated in section shown in Figure
9. The vertical distributions of the volume concentration are shown in Figure 10. Volume concentrations of large sphereparticles on the surface are greater compared to gravels, as shown in Figures 6, 7 and Figure 10. It is one reason thatengagement effect of gravels is strong and do not move actively. The other reason is that the longer travel distance ofparticles are easy to receive particles sort. We defined 0.3m from surface as particle moving layer in the slope failure. The
porosity in the moving layer was larger than that before the failure, as shown in Figure 3 and Figure 10. There are two
reasons. One reason is that vertical sorting occurred on surface, large particles were seen there mainly, and surface
porosity became large. The other reason is that moving of particles made gap between them and small particles fall into
there. Additionally increase in porosity on spheres (x=2.0m) surface within z=0.8-1.2m is gentler than that of gravels(x=2.0m), as shown in figure 10. Due to the engagement effect, gravel particles can stop in the high porosity state.
However, about tip particles, vertical porosity distributions of sphere (x=4.0m) and gravel (x=3.5m) are similar, because
particle height at the tip is nearly equivalent to particle size. Moreover, we confirm that tip particles consist of many large
particles, porosity become large, as shown in Figure 10. It is easy for large particles to gather on the tip. The reason is that
due to vertical sorting, volume concentrations of surface become small in comparison with lower part, and so large particle
in surface is easy to move in comparison with small particle in lower part.

8

E-proceedings of the 36th IAHR World Congress
28 June – 3 July, 2015, The Hague, the Netherlands

0

2

4

t (s)

6

8

0

10

t(s)

4

6

8

10

0

0
‐0.2

‐0.2

Sphere

‐0.4

‐0.4

‐0.6

‐0.6

Δz/Δx

Δz/Δx

2

‐0.8
‐1
d=40mmd=120mm‐1.2

Gravel‐0.8‐1d=40mmd=120mm

‐1.2

‐1.4

‐1.4

Figure 8. Time changes in the ratio of vertical travel distance ∆z and horizontal travel distance ∆x of spheres and gravels.

▽water surface
①

②

2m

4m

0m

Sphere

8m

③

④

2m

3.5m

Gravel8m

Figure 9. Vertical distribution of volume concentrations.

Total
合計

Section ① x=2m

0.80
0.60
0.40

0.40

0.20
0.00
0.0

0.2

0.4

0.6

0.8

Volume concentrations

1.0

0.0

0.2

0.4

0.6

0.8

1.0

Volume concentrations
SphereD=120mm

1.20

Section ② x=4m

0.80

0.00

D=90mmDistance from
水路底面からの高さ(m)channel bottomsz （ｍ）

Distance fromchannel bottoms (m)1.20

D=40mmSection ③ x=2m

0.80

Section ④ x=3.5m

0.60

0.80
0.40

0.40
0.20

0.00

0.00

0.0

0.2

0.4

0.6

0.8

1.0

Volume concentrations
体積割合

0.0

0.2

0.4

0.6

0.8

1.0

Volume concentrations
GravelFigure 10. Vertical distributions of volume concentrations after the failure.

4.

CONCLUSIONSWe discussed the effects of particle shapes and sizes on vertical sorting, angle of repose, changes of particle size
distribution on the slope surface, shapes of slope failure and the porosity of particles in slope failure. Because gravels havedifferent shapes compared to spheres, differences arise between them in resistance among particles, restraint of verticalsorting and vertical distribution of porosity in surface. Main conclusions drawn from the investigation are described asfollows:
1. The angle of repose of the submerged sphere particles was φ =22˚ (tanφ= 0.4), and that of the gravel particles was φ= 31˚ (tanφ = 0.6). The experimental result demonstrated that tanφ of the submerged gravel particles was 1.5 times(0.6/0.4) larger than that of the sphere particles due to the effect of shapes of gravel particles.

9E-proceedings of the 36th IAHR World Congress,
28 June – 3 July, 2015, The Hague, the Netherlands

2. The slope failure of sphere particles progressed smoothly with time from the upper part to the lower part, and formed astable slope. On the other hand, the slope failure of gravel particles is different from spheres. The lower part of gravel
particles progressed more slowly due to large engagement effects among particles compared to spheres. Then, upper
particles of gravel particles fell down.
3. The average sphere particle size (dm=81mm) on surface after the failure was larger than initially packed one(dm=70mm). Similarly, gravel particle size (dm=75mm) after the failure is larger than initial one (dm=70mm). We
investigated time changes in the ratios of average vertical travel distances ∆z of particle group to average horizontal
travel distances ∆x of particle group in the slope failure. The ratio ∆z/∆x of large particles is larger compared to small
particles in both sphere and gravel slope failures. Time of the large difference of spheres ∆z/∆x between small
particles and large particles is longer than that of gravels.
4. Porosities in the moving layer around 0.3m from the surface after failure was larger than that before the failure. Due tothe engagement effect, gravel particles can stop in the high porosity condition in comparison with sphere particles.
REFERENCES
Cundall, P.A. and Strack, O.D.L. (1979). A discrete numerical model for granular assemblies, Geotechnique, Vol.29, No.1,
pp.47-65.
Fukuoka, S. and Abe, T. (2007). Mechanism of low water channel formation and the role of grain-size distribution in gravelbed rivers, 10th International Symposium on River Sedimentation, Moscow, Russia.
Fukuda, T., Fukuoka, S. and Uchida, T. (2012). Mechanism of stones moving in streams and three-dimensional numericalmodeling. Journal of JSCE, B1 (Hydraulic Engineering), Vol.68, No.4, pp.I_937-I942, in Japanese.
Fukuoka, S., Fukuda, T. and Uchida, T. (2014). Effects of sizes and shapes of gravel particles on sediment transports andbed variations in a numerical movable-bed channel., Advances in Water Resources, Volume 72, p. 84-96.
Matsushima, T., Katagiri, J., Uesugi, K., Tsuchiyama, A. and Nakano, T. (2009). 3D Shape Characterization and ImageBased DEM Simulation of the Lunar Soil Simulant FJS-1, J.Aerospace Engineering, Vol. 22, No. 1, pp. 15-23.
Osada, K. and Fukuoka, S. (2012). Two-dimensional riverbed variation analysis method focused on the mechanism of
sediment transport and the bed surface unevenness in stony-bed rivers, Journal of JSCE, B1 (Hydraulic Engineering),
Vol.13, pp.339-344, in Japanese.
Parker, G. (1986). On armoring, Journal of JSCE, No. 375, pp. 17-27, in Japanese.
Tsuji, Y., Tanaka, T. Ishida, T. (1992). Lagrangian numerical simulation of plug flow of cohesionless particles in ahorizontal pipe (1992) Powder Technol. 71; 239-250 http://dx.doi.org/10.1016/0032-5910(92)88030-L.
Ushijima, S., Fukutani, A. and Makino, O. (2008). Prediction method for movements with collisions of arbitrary-shaped
objects in 3D free-surface flows, JSCE Journal B, Vol.64 No.2 pp.128-138, in Japanese.

10