Effects of sizes and shapes of gravel particles on sediment transports and bed variations in a numerical movable-bed channel

Advances in Water Resources 72 (2014) 84–96

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Effects of sizes and shapes of gravel particles on sediment transports and bed variations in a numerical movable-bed channel S. Fukuoka ⇑, T. Fukuda, T. Uchida Research and Development Initiative, Chuo University, 1-13-27 Kasuga, Bunkyo-ku, Tokyo 112-8551, Japan

a r t i c l e

i n f o

Article history: Available online 13 June 2014 Keywords: Numerical movable-bed channel Three-dimensional computation Sediment transport Gravel shape Hydrodynamic force Channel-bed configuration

a b s t r a c t The paper presents a numerical movable-bed channel capable of simulating three-dimensional motions of flows and gravel particles in different shapes. At first, the numerical channel was tested against results of fixed-bed channel experiments in which gravel particles were transported. Simulated particle motions were validated in comparison with those in the laboratory experiment. Next, numerical movable-bed experiments with sphere particles and gravel particles were conducted. The results of these experiments clearly elucidated the difference in motion between the large and the small particles, effects of shapes of gravel particles on sediment-transport rates, and hydrodynamic forces and contact forces at incipient motion and at settling. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Natural gravel-bed rivers display a large variety of particles sizes and shapes which often include cobbles. Although the sediment-transport rate in gravel-bed rivers has been discussed [1–8], the effect of sizes and shapes of gravel particles on sediment transport and river-bed variation have not been clarified sufficiently. The sediment-transport process in gravel-bed rivers during floods is explained that, at first small particles on the bed surface move, and next large gravel particles are exposed, resisting hydrodynamic forces primarily [1]. Exposed large particles do not move easily in large floods, but roll intermittently. However, studies on sediment transport for mixed-size sediment [3,9] in the laboratory experiments have been conducted primarily without large particles, such as cobbles. Such sediment-transport mechanism seems to differ from that of gravel-bed rivers. The bed-variation analysis using sediment-discharge formulas for mixed-size sediment [9] cannot explain adequately the variations of the bed levels and grain size distributions on the bed surface in gravel-bed rivers because the sediment-discharge formulas cannot evaluate the sheltering effect of large particles on small particles and hydrodynamic forces on projecting large particles [4,5]. Bed-surface patches caused by formation of large-particle clusters affect sediment transport in gravel-bed rivers [6–8]. Hendrick et al. [6] made field observations of formation of large-particle clusters on a gravel bed in several flow events. They demonstrated ⇑ Corresponding author. Tel./fax: +81 3 3817 1625. E-mail address: [email protected] (S. Fukuoka). http://dx.doi.org/10.1016/j.advwatres.2014.05.013 0309-1708/Ó 2014 Elsevier Ltd. All rights reserved.

that large-particle clusters were formed adjacent to riffles related to river-bed topographies. Piedra et al. [7] made experiments with mixed-size movable-bed, and investigated the relationship between tractive forces and ratio of large-gravel cluster occupying area to the bed surfaces. However, laboratory experiments and field observations have difficulty in investigating details of gravel-transport mechanism and arrangements of gravel particles on the bed due to the hydrodynamic forces. In recent years, the progress in computing performances has made possible to analyze the transport mechanism of individual particle motions by a Lagrangian method [10–12]. In most of these analyses, hydrodynamic forces on particles are calculated using drag coefficient [11,12]; therefore, evaluations of hydrodynamic forces remain inaccurate. A numerical method for solid–liquid multiphase flow has been developed in which hydrodynamic forces are calculated directly by calculating three-dimensional fluid motion around a solid object using finer computational cells than the object size [13,14]. Harada et al. [15] conducted large eddy simulation of solid–liquid multi-phase flow for drifting particles on a movable-bed of spheres in three different sizes and investigated vertical sorting processes in oscillatory flows. However, characteristic structural arrangements of particles in gravel-bed rivers, such as imbrications [16], indicate that motions of particles near the bed surface are strongly affected by the sizes and the shapes of particles. This natural image suggests that the investigation of particle motions in gravel-bed rivers should take irregular particle shapes into consideration. This paper presents numerical movable-bed channel composing different particle sizes and shapes where three-dimensional fluid

85

S. Fukuoka et al. / Advances in Water Resources 72 (2014) 84–96

motion was simulated in an Eulerian computational method for solid–liquid multiphase flow and gravel particles motion by a Lagrangian method. We use gravel particles with different sizes and shapes made of 8–10 small spheres superposed, as shown in Fig. 1 [17]. Eulerian descriptions of particle motions [18–20], especially descriptions of depth-integrated sediment transport [18,19], are practical for simulating sediment transport over a large area. However, the Eulerian description for the motion of particles of various sizes and shapes requires clear understanding of interactions among particle and particles themselves, flow and particles and the bed composition. We believe that flow and sediment-particle interaction dynamics of the Eulerian approach remains unsolved. This paper attempts to evaluate the dynamic interactions of particles themselves and particle–fluid which consider the effect of particle-size distributions, shapes of particles, and the bed structures. At first, the developed numerical channel was applied to the results of a fixed-bed channel experiment (hereinafter referred to as the real-scale experiment) in which motions of real-scale gravel particles in streams were measured [21]. Simulated particle motions were validated by comparing with those in the physical experiment. Next, two numerical experiments were conducted using the numerical movable-bed channel where mixed-size sphere particles and modeled gravel particles were used as the bed materials. The flow and particle motion and hydrodynamic forces on particles were predicted and visualized. The effects of size and shape of particles on motion of gravel particles, sorting mechanism of particles on the bed surface and flow structures near the bed were analyzed.

Particle

Different colors indicate different density Fig. 2. Concept of fluid motion analysis.

Sij ¼

  1 @ui @uj þ 2 @xj @xi

ð3Þ

m ¼ l=q

ð4Þ qffiffiffiffiffiffiffiffiffiffiffiffi

mt ¼ ðC S DÞ2 2Sij Sij

ð5Þ

where ui: ith 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, q: density, l: dynamic viscosity, gi: gravitational acceleration, mt: SGS turbulent viscosity, CS: Smagorinsky constant, D: computational grid size. Physical property / (density q, dynamic viscosity l) and averaged velocity ui were calculated as volume-averaged values and mass-averaged values, respectively, as follows:

/ ¼ a/s þ ð1  aÞ/f ;

/f ¼ f /l þ ð1  f Þ/g

ð6Þ

2. Numerical model

ui ¼ faqs usi þ ð1  aÞqf ufi g=q

In the present numerical model, fluid motions were simulated with the Eulerian method [14] and gravel motions were simulated with the Lagrangian method. To take into account the effect of solid phase on liquid phase, fluid motions were simulated by the governing equations of multiphase flows, as illustrated in Fig. 2. Fluid dynamic forces on particles were evaluated directly by integrating the forces on a particle region in the multiphase flow. Gravel particles with various shapes are made by the superposition of several small spheres, as shown in Fig. 1 [17]. Motions of the particles were simulated as rigid bodies.

where f: fluid volume fraction in fluid calculation cell, a: solid volume fraction in fluid calculation cell, suffixes l, s and g denotes liquid phase, solid phase, and gas phase, respectively. uf is fluid velocity, and us is solid velocity. Continuity equation and momentum equation were solved by the SMAC scheme with staggered grids.

ð7Þ

2.2. Evaluation of solid rigidity in fluid calculation

2.1. Governing equations of the fluid motion

The solid rigidity was evaluated by simulating particle motion as rigid body and setting average density q and average velocity ui with Eqs. (6) and (7). Solid velocity usi in Eq. (7) was calculated by the following equation:

The fluid motion was simulated by one-fluid model for the solid–liquid multiphase flow using the Smagorinsky model as the subgrid turbulence model:

us ¼ r_ G þ x  r f

@ui ¼0 @xi

ð1Þ

@ui @ui 1 @P @ þ uj ¼ gi  þ f2ðm þ mt ÞSij g @t @xj q @xi @xj

ð2Þ

ð8Þ

where rG: position vector of the center of gravity of the particle, x: angular velocity vector of particle, rf: position vector measured from particle’s center of gravity to the calculation point of us. Dot notation denotes the time derivative. The solid volume fraction in fluid calculation cell a in Eqs. (6) and (7) was calculated by dividing a fluid calculation cell into sub-cells and counting the number of

Particle surface Sub-cell

Fluid-computation cell Fig. 1. A model gravel particle constructed with multiple spheres.

Fig. 3. Evaluation of particle volume with sub-cells.

86

S. Fukuoka et al. / Advances in Water Resources 72 (2014) 84–96

sub-cells that were included in a gravel particle, as illustrated in Fig. 3.

coordinates system. In these coordinates transformation between global coordinates system and local coordinates system, the quaternion was used instead of transformation matrix R [14].

2.3. Evaluation of water-surface variation 2.6. Evaluation of fluid forces The water-surface variation was evaluated with continuity equation of fluid volume fraction f based on the VOF method [22].

@f @fui þ ¼0 @t @xi

ð9Þ

Gas motion was not calculated to reduce the computation time and fluid motion was calculated under the zero-pressure condition applied at the water surface as the boundary condition.

Fluid forces on a particle was evaluated by integrating forces on a particle region in the multiphase flow [14].

F f ;i ¼

Nf ;i ¼

 Z  @P @  þq f2ðm þ mt ÞSij g dX @xi @xj Xs 

Z

eijk rf ;j 

Xs

2.4. Modeling of particles in various shapes Particles of various sizes and shapes were made by the superposition of several small spheres to detect correctly contact and collision points between particles. It was difficult to calculate shapes of lapped part of spheres geometrically; therefore, the evaluation of rigid body properties such as mass, center of gravity and tensor of momentum inertia, were of importance. We calculated rigid-body properties of particles by using sufficiently small cells, and numerical integration considered volumes and points of cells included in individual particles, as shown in Fig. 4. Tensors of momentum inertia Ir of local coordinate system fixed on the gravel particles were not changed in motion; therefore, rigid-body properties of particles were calculated only once at the beginning of the simulation for every particles, respectively. 2.5. Governing equations of the particle motion in flows Particle motions were calculated by the momentum equations and the angular momentum equations of rigid bodies:

M €rG ¼ Mg þ F f þ F c

ð10Þ

n

x_ r ¼ I 1 R1 ðN f þ N c Þ  xr  I r xr r

o

ð11Þ

where bold face letters indicate vector tensor and matrix, M: mass of particle, rG: position of center of gravity, g: gravitational acceleration, F: forces on particles surfaces, N: torque on particles, x: angular velocity, R: transformation matrix (from the local coordinates system to the global coordinates system), and I: tensor of momentum inertia. Suffixes f and c: fluid force and contact force, respectively, suffix r: components in local coordinates systems of each particle. In solving angular momentum equations, next time step angular velocities in local coordinates system xr were first solved by Eq. (11). Then, second, angular velocities in global coordinates system x were calculated with transformation matrices R. Next time step transformation matrices R were set in consideration of the rotation in time step Dt with the angular velocities x. Next time step small sphere positions are set with next time step transformation matrices R, and small sphere positions in local

0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0

Fig. 4. Cells image for calculation of rigid body properties.

 @P @ þq f2ðm þ mt ÞSkl g dX @xk @xl

ð12Þ

ð13Þ

where Ff,i: ith component of fluid forces, Nf,i: ith component of the torque caused by fluid forces, rf,i: position vector from the center of gravity of the particle to the fluid calculation cells, Xs: area occupied by a particle, and eijk Levi–Civita symbol. 2.7. Evaluation of contact force Contact forces acting on particles were calculated by the contact detection of each small sphere composing a gravel particle, as shown in Fig. 5. Contact forces and torques on centers of gravity of particles were calculated by the summation of the contact forces calculated concerning each small sphere:

Fc ¼

X

F cp;n ; N c ¼

X

rcp;n  F cp;n

ð14Þ

where, Fc, and Nc: contact force and torque acting on the center of gravity of a particle, Fcp,n: contact force on each sphere, rcp,n: position vector from particle’s center of gravity to the contact point. Contact forces between particles were calculated by the Distinct Elements Method [23]. The spring constants k and coefficients of dashpot c were calculated by Eqs. (15)–(19) [24].

(   2 )13 4 r1 r2 E kn ¼ en 9 r1 þ r2 1  pos2

ð15Þ

ks 1 ¼ kn 2ð1 þ posÞ

ð16Þ

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m1 m2 cn ¼ 2h kn m1 þ m2

ð17Þ

pffiffiffiffiffi c s ¼ c n s0

ð18Þ

ln b h ¼  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2 þ ðln bÞ2

ð19Þ

s0 ¼

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 and two orthogonal directions on the tangential plane, respectively, and b: coefficient of restitution. The spring constants k was derived from the relationship between

Fcp ,1

rcp ,1

Contact detection and calculation of contact forces are conducted with each small sphere composing a gravel particle

Center of gravity of a particle Fig. 5. Evaluation of contact forces.

S. Fukuoka et al. / Advances in Water Resources 72 (2014) 84–96

87

the normal force and displacement of Hertzian contact theory [24]. Dashpot coefficient h (Eq. (17)) was derived from damping oscillation of two objects (m1, m2). h is related to the coefficient of restitution b [25] and takes the value 0.11 corresponding to the value of the coefficient of restitution, 0.7, in this calculation. Values of E, pos, and h were selected in consideration of the comparison between particle motions in the numerical channel and those in the real-scale experiment described in Section 4.

shapes of gravel particles. Figs. 7a and 7b shows that gravel particles share the same plane surfaces with others, although sphere particles met with others at a point. Therefore, it can be explained that gravel particles being about to slide or rotate were supported tightly by sharing the plane surfaces with surrounding particles. We call this the engagement effect of gravel particles.

3. Shape and mechanical property of gravel particles

4.1. Condition of numerical analyses

3.1. Shape of gravel particle

The validation of the developed numerical channel was checked by applying measured three-dimensional motions of flow and particles of the real-scale experiment in a fixed-bed channel conducted by Fukuoka et al. [21]. Moving trajectories of particles in the experiment were measured by the image analysis of video footage through the glass wall section. Figs. 8 and 9 show the experimental channel that was 45 m long and a bed slide of 1:20 and its cross section, respectively. Two experiments were conducted. The first test measured the motion of every single particle in streams and the second that of a cloud of particles. The developed numerical channel was applied to a single particle moving in the real-scale experiment using no. 1 gravel particle (shown in Fig. 6) and a spherical particle. The effects of the particle shapes on the single-particle motions were estimated by comparing trajectories of the gravel particle and the spherical particle in the numerical channel, and those of gravel particles in the laboratory experiment. Their results are presented in Section 4.2. The developed numerical bed was applied to the cloud of particles. The validation of particle motions under the intensive dynamic interactions between particles and flow was checked by the comparison between simulated and experimental vertical distributions of particle velocities and those of probabilities of particles in existence. The detailed results are shown in Section 4.3. We applied our numerical model to particle motions over the entire length (38 m long) of the laboratory flume. In the numerical analyses, the coordinates system was defined as the x axis in the longitudinal direction and the z axis in the vertical direction of the channel bed, as defined in Fig. 8. The numerical fixed-bed channel was packed with spheres with different radii and heights varying in the range of 0.02 m in order to reproduce the roughness of the experimental channel bed. A discharge of 0.5 m3/s was given at the upstream end of the channel, and a zero-pressure condition was given at the downstream end. In the application to the cloud of moving particles, gravel particles of five different sizes (25, 35, 50, 75, and 105 mm) were supplied to the numerical channel so that particle-size distributions of the numerical channel reproduced those of the laboratory experiments, as shown in Fig. 10. Twenty different types of particles (5 sizes  4 shapes) were used in the simulation. A total of 360 kg particles for 4 s was supplied at the upstream end of the simulation channel according to the laboratory experimental condition. The parameters used in these simulations are shown in Table 3.

Four different shapes of gravel particles were prepared to investigate their effects on motions in flows. Gravel particles were made to simulate shapes of real gravel particles shown in Fig. 6. Diameters of particles in different shape were given as diameters of spheres with the same volume. To investigate properties of particles shapes, three orthogonal axial lengths a, b and c (lengths of the longest, the intermediate, and the shortest axes of individual particles, respectively) were measured where the longest axis was first determined, then, the shortest axis was determined in the plane perpendicular to the longest axis and lastly the intermediate axis perpendicular direction to two other directions was determined. Table 1 shows a, b, and c of gravel particles nondimensionalized by their diameters, and Corey shape factor c/(ab)1/2. 3.2. Effects of particle shape on mechanical property of bed materials In order to understand effects of particle shapes on the mechanical property as granular media, angles of repose of spheres, and gravel particles in water were investigated. Particles size distributions and simulation parameters were the same as those in the numerical movable-bed experiment to be mentioned in Section 5. We packed particles in the numerical channels (length: 5 m, width: 1 m, and depth: 1.5 m) so that the length, the width, and the thickness of the particle group were 1 m each using the wall to prevent the particle group from sliding. Then, the wall was removed instantly. Figs. 7a and 7b shows the angles of repose of spheres and gravel particles. The angle of repose of the submerged spheres was u = 22° (tan u = 0.4), and that of the gravel particles was u = 31° (tan u = 0.6). The test result demonstrated that angle of repose (tan u) of the submerged gravel model particles was 1.5 times (0.6/0.4) larger than that of the spheres due to the effect of

1

2

3

4

1

2

3

4

4. Application to real-scale experiments

4.2. Single particle movement

Fig. 6. Comparison of shapes of simulated and real gravel particles.

Table 1 Lengths of a, b and c axes of model gravel particles. Shape no.

1

2

3

4

a b c Shape factor

1.26 0.98 0.88 0.79

1.29 1.06 0.81 0.69

1.49 0.89 0.76 0.66

1.36 0.99 0.78 0.67

Simulated trajectories of a single particle transported by the flow in the numerical channel are shown in Figs. 11 and 12. Motions of a gravel particle and a sphere were displayed on the x–z plane and the x–y plane, respectively. These figures also show representative trajectories of a single particle measured in the realscale experiments. The simulated gravel particle saltated more frequently, and moved downstream widely across the channel. On the other hand, the simulated sphere moved downstream just along the centerline of the channel where the lowest part of the flume

S. Fukuoka et al. / Advances in Water Resources 72 (2014) 84–96

z (m)

88

φ =22 ˚, tanφ =0.4

0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

y (m)

Cumulative percentage (%)

Fig. 9. Cross section of experimental flume.

Fig. 7a. Angle of repose of spheres in water.

φ =31 ˚, tan φ =0.6

100

Exp. Cal.

80 60 40 20 0 10

100

1000

Particle size (mm) Fig. 10. Particle-size distributions.

4.3. A cloud of particles movement Analyses of the numerical channel were also applied to the results of movements of a cloud of particles in the real-scale experiment [21]. A snapshot of simulated particle’s movement in the stream is illustrated in Fig. 13. In the experiment, all of the supplied particles reached the downstream end of the flume in about 20 s; however, in the simulated results, some small particles remained at the supplying section even after 20 s. Takiguchi et al. [26] clarified that the fluid-computation cell size smaller than 1/ 8 of a particle size was required for an accurate simulation of the drag force acting on a fixed sphere accompanied by unsteady vortex shedding in the uniform flow. The fluid-computation cell size was 10 mm in this simulation (the fluid-computation cell size was 1/2.5 of the smallest particles size), as shown in Table 2. This fluid-computation cell size was larger than the suitable size for calculating fluid motions around small particles. It was difficult to make fluid-computation cell sizes smaller because of large computational loads required. The larger particles could move over small particles by receiving proper hydrodynamic forces because the fluid-computation cell size was suitable. They were then transported downstream. Momentum exchanges of flows acting on deposited particles in the supplying section were underestimated because relatively large fluid-computation cell sizes were used for small particles. Therefore, a portion of small particles was left within the supplying section. We investigated mutual interactions between moving particles and flows at sufficiently downstream sections. Vertical distributions of the experimental and the numerical moving particle velocities are compared in Fig. 14. This figure also shows a simulated

Fig. 7b. Angle of repose of gravel particles in water.

Measuring section of moving particles in the experiment (Glass wall section)

45m 20m

z

4m

y x

Fig. 8. Real-scale laboratory experimental flume and its coordinates system.

cross section was located, as shown in Fig. 9. The difference in trajectories between gravel particles and spheres arose from the fact that gravel particles had larger frontal areas than those of sphere particles with the same volume, and collided with the flume bottom in irregular manners because of different shapes. It, therefore, can be concluded that the shape of particles was an important factor for particle motions in streams.

Distance from Channel Bottom Z(m)

Exp. (D=20-30 mm) Exp. (D=90-120 mm)

Exp. (D=30-40 mm) Cal. Gravel (D=105 mm)

Exp. (D=40-60 mm) Cal. Sphere (D=105 mm)

Exp. (D=60-90 mm)

0.25 0.20 0.15 0.10 0.05 0.00 15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

Longitudinal Distance X (m)

Fig. 11. Trajectories of a single particle movement on x–z plane.

34

35

36

37

38

39

40

89

S. Fukuoka et al. / Advances in Water Resources 72 (2014) 84–96

Lateral Distance Y(m)

Exp. (D=20-30 mm) Exp. (D=90-120 mm)

Exp. (D=30-40 mm) Cal. Gravel (D=105 mm)

Exp. (D=40-60 mm) Cal. Sphere (D=105 mm)

Exp. (D=60-90 mm)

0.20 0.15 0.10 0.05 0.00 -0.05 -0.10 -0.15 -0.20 15

16

17

18

19

20

21

22

23

24 25 26 27 28 29 Longitudinal Distance X (m)

30

31

32

33

34

35

36

37

38

39

40

Fig. 12. Trajectories of a single particle movement on x–y plane.

Water Surface Flow

24m 26m Measuring Section

28m

Yellow ····· : D =105 mm Green ········· : D = 75 mm Light Blue·· : D = 50 mm Red ············ : D = 35 mm Blue ··········· : D = 25 mm The flume wall is not shown in the figure.

Fig. 13. A snapshot of the simulated particles in the measuring section.

Table 2 Parameters used in numerical simulation.

Dx, Dy, Dz: fluid-computation cell size for a single particle movement simulation Dx, Dy, Dz: fluid-computation cell size for a cloud of particles movement simulation Dt: time step for fluid-motion simulation Dt0 : time step for particle-motion simulation qw: density of water qs: density of particles lw: viscosity coefficient of water ls: viscosity coefficient of particles CS: Smagorinsky coefficient E: elastic modulus pos: Poisson’s ratio h: Dashpot coefficient

0.02

m

0.01

m 4

1.0  10 1.0  106 1000 2650 8.9  104 8.9  104 0.173 5.0  1010 0.33 0.11

s s kg/m3 kg/m3 Pa  s Pa  s Pa

flow-velocity distribution. The simulated particle velocities were almost the same as the experimental values. Vertical distributions of existence probabilities of particles are shown in Fig. 15. The simulated existence probabilities of particles for D < 50 mm were larger near the bottom and smaller at higher elevations than those of the experiment. However, simulated probabilities for D > 60 mm had good agreements with the probabilities of the

Distance from channel bottom z (m)

0.25 D = 20 - 30 mm D = 30 - 40 mm

0.20

D = 40 - 60 mm

Exp.

D = 60 - 90 mm

0.15

D = 90 - 120 mm Flow Velocity

0.10

D = 25 mm D = 35 mm

Cal.

D = 50 mm

0.05

D = 75 mm D = 105 mm

0.00 0

1

2

3

4

Velocity (m/s) Fig. 14. Vertical distributions of particle and flow velocities.

experiment. The investigation of the vertical distributions of particle velocities and existence probabilities showed that it is possible to analyze the movement mechanism of particles and flows near the surface of gravel bed using the developed numerical movable-bed model.

5. Numerical movable-bed channel experiment 5.1. Numerical experiment conditions Two cases of the mixed-size movable-bed experiment were conducted in the developed numerical channel where spheres and gravel particles were used as the bed materials. Particle movement mechanism, sediment-transport rates and bed structures are discussed on particles in different sizes and shapes. The numerical movable-bed channel and definition of coordinates system used in the numerical analysis are shown in Fig. 16. The numerical channel designed was 15-m long, 1-m wide with a bed slope of 1:20 by considering flow conditions that were large enough to transport gravel particles, sufficiently long for bed forms to form, and accommodated the computational load. A water discharge of 0.5 m3/s was supplied at the upstream end of the channel, and the zero pressure condition was set at the downstream end, which were the same in the real-scale experiments. Gravel particles in four different shapes, as illustrated in Fig. 6, were prepared for bed materials of the numerical movable-bed channel. The investigation of the real-scale experiment revealed that the reproducibility of numerical simulations for smaller particle motions was not satisfactory. Therefore, the particle size distributions in the numerical movable-bed experiment were made of five different sizes of particles (40, 50, 70, 90, and 120 mm) so that the smallest size of particles was greater than that of the real-scale experiment. These particles were packed randomly to form the channel bed. The particle size distribution of the numerical channel bed thus produced is shown in Fig. 17. The same number and the volume of particles discharging out from the end of the channel were introduced back simultaneously to a section (x = 1–2 m) near the upstream end of the channel. Numerical analyses started at the same time as the beginning of the water supply. The starting time of the numerical experiment was when the discharge became almost constant.

Exp. ( D = 30 - 60 mm ) Cal ( D = 35, 50 mm )

0.20 0.15 0.10 0.05 0.00 0

10

20

30

Existence probability of particles (%)

0.25 Exp. ( D = 60 - 90 mm ) Cal. ( D = 75 mm )

0.20 0.15 0.10 0.05 0.00 0

10

20

30

Existence probability of particles (%)

Distance from the channel bottom z (m)

0.25

Distance from the channel bottom z (m)

S. Fukuoka et al. / Advances in Water Resources 72 (2014) 84–96

Distance from the channel bottom z (m)

90

0.25 Exp. ( D = 90 -120 mm ) Cal. ( D = 105 mm)

0.20 0.15 0.10 0.05 0.00 0

10

20

30

Existence probability of particles (%)

Fig. 15. Vertical distributions of existence probabilities of particles.

The numerical experiments lasted for 200 and 400 s for the sphere bed and the gravel bed, respectively. In fluid and gravel-particle movement calculations Dt were set at 5.0  104 and 5.0  106 s, respectively. The fluid-computation cell size was 0.01 m, and other parameters used were the same as the real-scale experiment, as shown in Table 2. 5.2. Bed variation by the flow Average flow depths, average flow velocities, and Froude numbers used at the final stage of the experiments in both sphere and gravel beds are shown in Table 3. There were practically no difference in values of both cases, and flows were supercritical. Time variations in the longitudinal profiles of laterally averaged-bed levels and water-surface profiles in sphere bed and gravel bed are shown in Figs. 18 and 19, respectively. These profiles display the values averaged over two seconds. In both cases, the water surface profiles varied in the same phase as the bed waves. Antidunes moved in the upstream direction; however, the developing speed of antidunes in the sphere bed was faster than that of the gravel bed. Antidunes appeared after about t = 150 s in the sphere bed, but they formed at about t = 300 s in the gravel bed. Sphere particles were found to move easily than gravel particles. Time variations of the sphere and the gravel-bed surfaces are shown in Figs. 20 and 21, respectively. Time variations of the sphere-bed and the gravel-bed levels are shown in Figs. 22 and 23, respectively. Figs. 20 and 21 demonstrate clearly sorting of particles

15m

0m

Table 3 Comparison of flow conditions for sphere and gravel particles. Case

Average flow depth (m)

Average flow velocity (m/s)

Froude number

Measurement time (s)

Sphere Gravel

0.273 0.274

1.83 1.83

1.12 1.12

t = 200 t = 400

and formation of large-particle clusters on the bed surface. Antidunes developed over the whole channel length, and large particles were found on the crests of antidunes and small particles were at the troughs. The height of antidunes of the gravel bed at t = 400 s was greater than that of the sphere bed at t = 200 s, as can be seen in Figs. 18 through 23. The reason was that spheres were easy to move due to the smaller engagement effect than gravel particles, and they deposited at troughs of antidunes more easily than gravel particles. 5.3. The process of large-particle cluster formation Trajectories of gravel particles visualized by emitting the light on moving gravel particles at every 1/10 s are illustrated in Fig. 24. Small particles were not able to climb over large projecting particles, but they took detours. Large particles were nearly in rectilinear motions over projecting particles, resulting in receiving greater flow velocities. Large particles, therefore, were not able to stop moving on the plane bed, but at projecting beds as high as diameters of large particles. They had a tendency to form largeparticle clusters on the bed. Relatively deep water passages were formed among large-particle clusters of the bed surface, and small particles meandered along deeper water passages. Characteristic motions of large particles and small particles resulted in the particle sorting of the gravel-bed surface.

10m

1m

0.3m

Slope of the initial movable bed: 1:20

15m

Cumulative percentage (%)

Fig. 16. Numerical movable-bed channel and coordinates system.

100 80 60 40 20 0 10

100

1000

Particle size (mm)

Fig. 17. Particle-size distributions of the numerical channel bed.

5.4. Vertical velocity distributions of flows and particles near the movable bed Streamlines of representative flows near the bed surface where large-particle clusters were formed at x = 8 m and where no clusters were formed at x = 12 m are shown in Fig. 25. Vertical distributions of the time-averaged flow velocities in the xdirection measured at the center of the channel at sections x = 8 m and x = 12 m are shown in Fig. 26. Velocities were averaged laterally over 0.15 m from the channel centerline which was about the largest particle size. Fig. 26 also shows time-averaged velocities of gravel particles moving over the lateral range of 0.15 m and log-law velocity distributions depicted by assuming ks = 70 and 240 mm. Values of ks = 70 mm and ks = 240 mm were equivalent to D60 and twice the maximum particle diameter, respectively. A shear velocity of 0.29 m/s used was the average value taken from x = 5 to x = 13 m. The location of x = 8 m

91

S. Fukuoka et al. / Advances in Water Resources 72 (2014) 84–96 0.80 Water level (50S) Bed level (50S)

0.70

Water level (100S) Bed level (100S)

Water level (150S) Bed level (150S)

Water level (200S) Bed level (200S)

z(m)

0.60 0.50

Initial bed level

0.40 0.30 0.20 0.10 2

3

4

5

6

7

8

9

10

11

12

13

14

15 x(m)

Fig. 18. Longitudinal profiles of averaged-water surface and bed levels in sphere-bed case.

0.80 0.70 z(m)

0.60

Water level (100S)

Water level (200S)

Water level (300S)

Water level (400S)

Bed level (100S)

Bed level (200S)

Bed level (300S)

Bed level (400S)

0.50

Initial bed level

0.40 0.30 0.20 0.10 2

3

4

5

6

7

8

9

10

11

12

13

14

15 x(m)

Fig. 19. Longitudinal profiles of averaged-water surface and bed levels in gravel-bed case.

120 mm

90 mm

70 mm 50 mm 40 mm

5m Before the experiment 10 m

5m t = 100 s 15 m 10 m

5m t = 200 s

15 m 10 m

15 m Fig. 20. Time variation of sphere-bed surface.

was the section where clusters of large particles were formed. Vertical velocity distributions did not follow the log law. On the other hand, where large-particle clusters were not formed, velocity distributions of the flow appeared to follow approximately the log law (ks = 70 mm). From these investigations, we

found that at the location forming large-particle clusters, velocity gradually increased from the bed surface in the vertical direction at a scale of a large-particle size because there were a number of large pores in the clusters, and large-particle clusters brought great resistance to flow.

92

S. Fukuoka et al. / Advances in Water Resources 72 (2014) 84–96

120 mm

90 mm 70 mm 50 mm 40 mm

5m Before the experiment 10 m

5m t = 100 s

15 m 10 m

5m t = 200 s 15 m 10 m

5m t = 400 s 15 m 10 m

15 m

t = 100 s

y (m)

t = 200 s

y (m)

Before the experiment

y (m)

Fig. 21. Time variation of gravel-bed surface.

Bed level distance from channel bottom z (m)

x (m)

y (m)

t = 200 s

t = 400 s

Bed level distance from channel bottom z (m)

x (m)

y (m)

t = 100 s

y (m)

Before the experiment

y (m)

Fig. 22. Time variation of sphere-bed level.

Fig. 23. Time variation of gravel-bed level.

93

S. Fukuoka et al. / Advances in Water Resources 72 (2014) 84–96

(b) Trajectories of small particles (D = 40mm, and 50mm)

(a) Trajectories of all particles.

Fig. 24. Trajectories of moving gravel particles.

5.5. Shear-stress distributions acting on bed surfaces

sf ;zx We converted hydrodynamic forces on particles to shear stresses on the x–y plane and investigated the relationship among the shear stresses, uneven bed-surface shapes, and sizes of bed-material particles. The shear stress exerted over particles by flow was obtained by integrating the x-directional hydrodynamic forces from the channel bottom to the top of a moving particle and dividing by the area of the x–y integrating area, as shown in Eq. (20).

o R R R n @P a  @x1 þ q @x@ j f2ðm þ mt ÞS1j g dxdydz RR ¼ dxdy

Shear-stress distributions of the hydrodynamic forces at t = 100 s and t = 400 s on the gravel bed that were calculated using Eq. (20) are shown in Fig. 27. The shear-stress distribution of the large hydrodynamic forces at t = 100 s is seen to distribute sparsely on the bed surface where large particles rested on the bed surface as shown in Fig. 21. On the other hand at t = 400 s, high shear stresses were concentrated and located at regular intervals in the

Flow

Flow

At x = 8 m

At x = 12 m 2.0 m/s

1.5 m/s

1.0 m/s

0.5 m/s

Flow velocities

z (m)

z (m)

Fig. 25. Streamlines of flows within 0.1 m from the gravel-bed surface.

0.60 0.55 0.50 0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10

ks=240mm ks=70mm Log law velocity distributions

0

0.5

1

1.5

2

2.5

Velocity in x direction (m/s)

At x = 8 m

3

3.5

ð20Þ

0.60 0.55 0.50 0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10

Log law velocity distributions ks =240mm ks=70mm

Flow velocities

0

0.5

1

1.5

2

2.5

3

Velocity in x direction (m/s)

At x = 12 m

Fig. 26. Vertical velocity distributions of flows and particles.

3.5

Particle velocity D = 120 mm D = 90 mm D = 70 mm D = 50 mm D = 40 mm

94

S. Fukuoka et al. / Advances in Water Resources 72 (2014) 84–96

t = 100 s y (m)

Shear stress, acting on gravel bed (Pa) x (m)

y (m)

t = 400 s

x (m)

0.8 0.6

20

Hydrodynamic force

Velocity Contact force

0.4

0

0.0 -0.2 -0.4

0 0.4 -0.4 0 0.4 -0.4 0 0.4 -0.4 0 0.4 time (s) time (s) time (s) time (s) Small particles Small particles Large particles Large particles (D = 40 mm, 50 mm, 70 mm) (D = 40 mm, 50 mm, 70 mm) (D = 90 mm, 120 mm) (D = 90 mm, 120 mm) Start of particle motion

-20

nondimensionalized hydrodynamic force and contact force in x direction

x - direction velocity (m/s)

Fig. 27. Distributions of shear stresses, sf,zx, acting on the gravel bed.

Rest of particle motion

0.8 0.6

20

Hydrodynamic force Contact force

0.4 Velocity

0

0.0 -0.2 -0.4

0 0.4 -0.4 0 0.4 -0.4 0 0.4 -0.4 0 0.4 time (s) time (s) time (s) time (s) Small particles Large particles Small particles Large particles (D = 40 mm, 50 mm, 70 mm) (D = 90 mm, 120 mm) (D = 40 mm, 50 mm, 70 mm) (D = 90 mm, 120 mm) Start of particle motion

-20

nondimensionalized hydrodynamic force and contact force in x direction

x - direction velocity (m/s)

Fig. 28. Time variations of particle velocities, hydrodynamic forces and contact forces when particles were about to move and settle in the sphere-particle bed.

Rest of particle motion

Sediment transport rate ×10-3 (m2/s)

D=120mm D=50mm

D=90mm D=40mm

D=70mm Total

5.0 4.0 3.0 2.0 1.0 0.0 0

50

100

150

200

time ( s ) Fig. 30. Sediment-transport rates of sphere particles in different sizes.

channel where antidunes were formed. Bed topographies of antidunes resulting from large-particle clusters apparently produced a larger resistance to channel flows. 5.6. Forces acting on particles at incipient motion and at settling Time variations of x-direction particle velocities, x-direction contact forces, and hydrodynamic forces on sphere particles that

Number of transported particles ( 1/s )

Fig. 29. Time variations of particle velocities, hydrodynamic forces and contact forces when particles were about to move and settle in the gravel bed.

Total D=70mm

D=40mm D=90mm

D=50mm D=120mm

20 15 10 5 0 0

50

100 time ( s )

150

200

Fig. 31. Number of sphere particles leaving from downstream end of the channel.

were nondimensionalized by the x-component of submerged weight of a sphere are shown in Fig. 28. These data are ensemble means of data over t = 180 s to t = 200 s. Similarly, ensemble means for gravel particles over t = 380 s to t = 400 s are shown in Fig. 29. The time zeros in Figs. 28 and 29 were chosen at the time when particles started to move and came to a rest, respectively. It is clearly seen in Figs. 28 and 29 that when particles started to move, hydrodynamic forces on gravel particles were much greater than those on spheres because gravel particles must overcome

95

S. Fukuoka et al. / Advances in Water Resources 72 (2014) 84–96

Sediment transport rate ×10-3 (m2/s)

D=120mm

D=90mm

D=70mm

D=50mm

D=40mm

Total

5.0 4.0 3.0 2.0 1.0 0.0

0

50

100

150

200 time ( s )

250

300

350

400

Number of transported particles ( 1/s )

Fig. 32. Sediment-transport rates of gravel particles in different sizes.

Total

D=40mm

D=50mm

D=70mm

D=90mm

D=120mm

20 15 10 5 0 0

50

100

150

200 time ( s )

250

300

350

400

Fig. 33. Number of gravel particles leaving from downstream end of the channel.

engagement effects present between particles. On the other hand, gravel particles accelerated faster than sphere particles due to larger projection areas against the flow. Hydrodynamic forces on small particles decreased after settling. On the other hand, hydrodynamic forces on large particles did not decrease after settling. This difference arose from the fact that smaller particles were able to stop at small velocity regions behind stationary large particles, while larger particles were able to stay at high elevations and to receive larger hydrodynamic forces. Therefore, larger particles could not stop easily but rest by the collision or support of surrounding large particles. 5.7. Sediment-transport rate Computed sediment-transport rates of spheres in different sizes and the number of spherical particles leaving from the downstream end of the channel per unit time are shown in Figs. 30 and 31, respectively. Those of the gravel particles are shown in Figs. 32 and 33, respectively. The experiments of the numerical movable-bed channel demonstrated that sediment-transport rates of larger particles were greater than smaller particles in both spheres and gravel particles. However, sediment-transport rates of smaller particles were greater compared to those of large particles in a field measurement conducted in a gravel-bed river (Joganji River in Japan) where the particle-size distribution of the bed materials was widely distributed [27]. The difference is explained by following reasons: particles smaller than 40 mm which were transported easily by the flow, and rollers for the motion of larger particles were not included in the numerical movable-bed channel. This numerical movable-bed channel proved to simulate appropriately the motion of particles in different sizes and shapes. Therefore, computation results of sediment-transport rates were correct as long as the present particle-size distribution of bed materials was used. Further numerical computations of sediment-transport rates will be

conducted for particle-size distributions which include smaller particles. The total sediment-transport rate of all particles decreased compared to that of the early stage as shown in Figs. 30 and 32. It was considered that after antidunes formed on the bed, particles became difficult to move over antidunes. The total sediment-transport rate of spheres was greater than that of the gravel particles, as shown in Figs. 30 and 32, while the number of transported particles smaller than 50 mm of gravel particles was greater than that of spheres, as shown in Figs. 31 and 33. This suggests that gravel particles make a difference in transport rates of particles in different sizes smaller. Gravel particles starting to move force surrounding particles to move together due to engagement effects among particles. The average sediment-transport rate of spheres per unit width after development of the antidunes was 1.8  103 (m2/s) at t = 150 s, and that of gravel particles 1.0  103 m2/s at t = 300 s. The unit-width sediment-transport rate of gravel particle was 0.56 (1.0  103/1.8  103) times smaller than that of spheres. This demonstrates that the effect of particle shapes on sediment-transport rate was significant. Note that the ratio of the angle of repose of submerged spheres to that of submerged gravel particles was 0.67 (0.4/0.6) and was close to the ratio of sedimenttransport rate of 0.56. Therefore, the angle of repose in submerged particles could be used as a rough estimate of the effect of particle shapes on the total sediment-transport rate.

6. Conclusions A numerical movable-bed channel capable of simulating threedimensional motions of flows and gravel particles in different sizes was developed. At first, the numerical channel was tested against real-scale experiments of gravel particle motion in streams. Next, numerical movable-bed channel experiments were conducted to investigate movement and interaction of flow and gravel particles

96

S. Fukuoka et al. / Advances in Water Resources 72 (2014) 84–96

near the movable-bed surface. Main conclusions drawn from the present investigation are described as follows: 1. The numerical channel developed in the present study provided us with good agreements with the results of the real-scale experiment with regard to trajectories of particles, vertical distributions of particle velocities, and existences probabilities of particles. 2. The numerical movable-bed channel experiment drew a new methodology to understand three-dimensional motions of flow and gravel particles near the bed surface, effects of shapes and sizes of particles on the movement and settling of gravel particles, sorting mechanism of the bed materials, and sedimenttransport rates. 3. Larger particles were found to be lifted higher in elevation as they roll and saltate, receiving larger flow velocities, and they did not settle easily. Larger particles had a tendency to rest by contact forces with surrounding large particles and to form large-particle clusters. Smaller particles moved along deeper water passages formed among large-particle clusters of the bed surface. This is an important mechanism for sorting of gravel particles on the bed surface. In the supercritical flow, clusters of large particles developed into antidunes where larger particles gathered around the crest of the antidunes and smaller particles around the trough. 4. The sediment-transport rate of gravel particles was roughly a half of that of spheres in the particle-size distributions used in the present experiments, indicating that effects of particle shapes on sediment-transport rate are significant. 5. Gravel particles starting to move force surrounding particles to move together due to engagement effects among particles. They had a tendency to make the differences in the transport rates of each particle size smaller.

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

References [1] Fukuoka S, Abe T. Mechanism of low water channel formation and the role of grain-size distribution in gravel-bed rivers. In: 10th international symposium on river sediment, Moscow Russia; 2007. [2] Parker G. Surface-based bedload transport relation for gravel rivers. J Hydraul Res 1990;28(4):417–36. http://dx.doi.org/10.1080/00221689009499058. [3] Wilcock PR, Crowe JC. Surface-based transport model for mixed-size sediment. J Hydraul Eng 2003;129:120–8. http://dx.doi.org/10.1061/(ASCE)07339429(2003)129:2(120). [4] Osada K, Fukuoka S, Analysis of two-dimensional riverbed variation in channels with stony beds. In: 11th International Symposium on River sediment, Stellenbosch, South Africa, CD-ROM; 2010. [5] Osada K, Fukuoka S. Two-dimensional riverbed variation analysis method focused on the mechanism of sediment transport and the bed surface unevenness in stony-bed rivers. J JSCE, B1 Hydraul Eng 2012;13:339–44. http://dx.doi.org/10.2208/jscejhe.68.1. [6] Hendrick RR, Ely LL, Papanicolaou AN. The role of hydrologic processes and geomorphology on the morphology and evolution of sediment clusters in

[22]

[23] [24]

[25]

[26]

[27]

gravel-bed rivers. Geomorphology 2010;114:483–96. http://dx.doi.org/ 10.1016/j.geomorph.2009.07.018. Piedra MM, Haynes H, Hoey TB. The spatial distribution of coarse surface grains and the stability of gravel river beds. Sedimentology 2012;59:1014–29. http://dx.doi.org/10.1111/j.1365-3091.2011.01290.x. Nelson PA, Bellugi D, Dietrich WE. Delineation of river bed-surface patches by clustering high-resolution spatial grain size data. Geomorphology 2014;205:102–19. http://dx.doi.org/10.1016/j.geomorph.2012.06.008. Ashida K, Michiue M. Study on hydraulic resistance and bed-load transport rate in alluvial streams. Proc JSCE 1972;206:59–69. http://dx.doi.org/10.2208/ jscej1969.1972.206_59. Gotoh H, Sakai T. Numerical simulation of sheetflow as granular material. J Waterway Port Coastal Ocean Eng 1997;123(6):329–36. http://dx.doi.org/ 10.1061/(ASCE)0733-950X(1997)123:6(329). Schmeeckle MW, Nelson JM. Direct numerical simulation of bedload transport using a local, dynamic boundary condition. Sedimentology 2003;50:279–301. http://dx.doi.org/10.1046/j.1365-3091.2003.00555.x. Calantoni J, Thaxton CS. Simple power law for transport ratio with bimodal distributions of coarse sediments under waves. J Geophys Res 2008;113: C03003. http://dx.doi.org/10.1029/2007JC004237. Kajihima T, Takiguchi S, Hamasaki H, Miyake Y. Turbulence structure of particle-laden flow in a vertical plane channel due to vortex shedding. JSME Int J, Ser B 2001;44(4):526–35. http://dx.doi.org/10.2208/jscejb.64.128. Ushijima S, Fukutani A, Makino O. Prediction method for movements with collisions of arbitrary-shaped objects in 3D free-surface flows. JSCE J B 2008;64(2):128–38. http://dx.doi.org/10.2208/jscejb.64.128. Harada E, Tsuruta N, Gotoh H. Large eddy simulation for vertical sorting process in sheet-flow regime. J JSCE, B2 Costal Eng 2011;67(2):I471–5. http:// dx.doi.org/10.2208/kaigan.67.I_471. Millane RP, Weir MI, Smart GM. Automated analysis of imbrication and flow direction in alluvial sediments using laser-scan data. J Sediment Res 2006;76:1049–55. http://dx.doi.org/ 10.2110/jsr.2006.098. Matsushima T, Katagiri J, Uesugi K, Tsuchiyama A, Nakano T. 3D shape characterization and image-based DEM simulation of the lunar soil simulant FJS-1. J Aerosp Eng 2009;22(1):15–23. http://dx.doi.org/10.1061/(ASCE)08931321(2009)22:1(15). Greco M, Iervolino M, Leopardi A, Vacca A. A two-phase model for fast geomorphic shallow flows. Int J Sediment Res 2012;27:409–25. http:// dx.doi.org/10.1016/S1001-6279(13)60001-3. Evangelista S, Altinakar M, Di Cristo C, Leopardi A. Simulation of dam-break waves on movable beds using a multi-stage centered scheme. Int J Sediment Res 2013;28:269–84. http://dx.doi.org/10.1016/S1001-6279(13)60039-6. Penko AM, Calatoni J, Rodriguez-Abudo S, Foster DL, Slinn DN. Threedimensional mixture simulations of flow over dynamic rippled beds. J Geophys Res Ocean 2013;118:1543–55. http://dx.doi.org/10.1002/jgrc.20120. Fukuoka S, Watanabe A, Shinohara Y, Yamashita S, Saitou K. Movement mechanism for a mass of gravel flow with high speed and estimation of channel bed abrasion. Proc Adv River Eng JSCE 2005;16:291–6. Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 1981;39:201–25. http://dx.doi.org/10.1016/00219991(81)90145-5. Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Geotechnique 1979;29(1):47–65. http://dx.doi.org/10.1680/geot.1979. 29.1.47. Tsuji Y, Tanaka T, Ishida T. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe. Powder Technol 1992;71:239–50. http://dx.doi.org/10.1016/0032-5910(92)88030-L. Tsuji Y, Kawaguchi T, Tanaka T. Discrete particle simulation of twodimensional fluidized bed. Powder Technol 1993;77:79–87. http:// dx.doi.org/10.1016/0032-5910(93)85010-7. Takiguchi S, Kajihima T, Miyake Y. Numerical scheme to resolve the interaction between solid particles and fluid turbulence. JSME Int J, Ser B 1999;42(3):411–8. http://dx.doi.org/10.1299/jsmeb.42.411. Soyama K, Ohkuma Y, Hatanaka Y, Asano F, Fukuoka S. Study on observation technique and assessment of bed load discharge. Adv River Eng 2011;17:5–10.