Numerical Movable Bed Channels Required for Investigationof Various Sizes and Shapes Particle Motions in Gravel Bed Rivers
Tomoo Fukuda1 and Shoji Fukuoka2
1

Research and Development Initiative, Chuo University
Kasuga 1-13-27, Bunkyo-ku, Tokyo 112-8551, JAPANt-fukuda@tamacc.chuo-u.ac.jp

2

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

1. INTRODUCTION
Sorting of particles and imbrication on the bed surface with particles of large variety of sizes
and shapes are seen in gravel-bed rivers and the mechanism of sediment transport is more
complicated than that in sand-gravel bed rivers. To discuss the mechanism of gravel particle
motions in streams, the Lagrangian numerical approaches are important because motions of
the flow and particles near the bed surface during floods can be hardly measured in rapid
flows. Most numerical studies in the Lagrangian method have used spheres and been carriedout in short channels (dozens of times the particle size) applying the periodic boundary
conditions. These approaches have not demonstrated important mechanisms of particle
motions of various shapes and development mechanisms of the bed structures such as sortingsof particles and bed waves. We developed a numerical movable bed channel capable ofsimulating three-dimensional motions of flow and particles of various sizes and shapes [1].
Irregular shape particles were made by the combination of several small spheres (as shown inFigure 1). This paper demonstrates the numerical movable bed channels required forinvestigations of flow and particle motions in gravel bed rivers with particles of various sizes
and shapes and presents the results of the numerical movable bed experiments.
2. NUMERICAL METHOD
2.1 Numerical method of fluid motions
Fluid motions are simulated with the governing equations of the solid-liquid multiphase flowwhere the regions of gravel particles are considered as a different density fluid [2]. We
applied the Smagorinsky model as the sub grid turbulence model:
ui
0
xi

(1)

ui
u1 P
2  t Sij 
 u j i  gi tx j xi x j(2) t  CS  2 2Sij Sij(3)where ui : i-direction averaged velocity including
gravel particle regions in a fluid calculation cell, ρ:density, P: sum of pressure and isotropic component- 223 -

1234Figure 1. Different shapes of particles.

5612 131417116 15 2 15 47826 71110 9 8334(a) Connecting cells numberswith sphere numbers

56712 14131112 1 2 5
3 473158685 1112 6 1314 7 1515 8 1617

16 171049

1 12 2 35(b) Spheres sorted by cells
number3 6 7 4 810(c) Maximum and minimum sorted
sphere number in each cell

Figure 2 Relationship between cell number and sphere number.
Table 1 Connecting gravel particle number withsphere number.

(a) List for gravel particles
referring to sphere number.

1526NSCG

SphereNo.

Sphere
No.

12345
1234
5
6

1234567891011

1234
5678
91011Gravel
Particle
No.
1
1
1
1
12
2
2222···MAXNSCG

Table 2 Relationship between sorted numberand sphere number.
Sorted No.
123Sphere No.
15165···2.2 Modeling of particles of various sizes andshapes

GravelParticleNo.

(b) List for spheres referring
to gravel number.

···of SGS stress, g: gravitational acceleration, Sij:
strain rate tensor, ν: kinematic viscosity, νt:
SGS turbulent viscosity, Cs: Smagorinsky
constant, ∆: calculation grid size. Physicalproperty (density ρ, dynamic viscosity μ) andaveraged velocity ui are calculated asvolume-averaged values and mass-averagedvalues, respectively. Fluid forces acting on aparticle are estimated with integrating forces
acting on a particle region in the multiphaseflow. Water surface variations are estimatedby calculating advection of water volume
fraction in fluid calculation cells with theDonor Cell Method.

···········

Table 3 List of contact candidate pair spheres.
Particles of various sizes and shapes are madePair list
ForcesContact
Listby the superposition of several small spheres as
acting on(A)Smaller(B)Larger
PointNo.
sphere A
Sphere No.
Sphere No.
shown in Figure 1. The motion of gravel
(a)
(b)
(c)
(d)
(e)
1
515
F1
X1
particles is simulated with the momentum
22
11F2
X2
equations and the angular momentum
equations of rigid bodies. Rigid body properties
Table 4 Relationship between sphere numberof gravel particles (mass, gravity center andand list number.
SphereNo.
List No.
tensor of momentum inertia) are calculated by1
numerical integrations in consideration of22
positions and volumes of sufficiently small
51cells (1.0×10-9m3) included in spheres
112
composing gravel particles. Gravel particles
151
refer to small spheres composing them andspheres do to gravel particles in this simulation,respectively. Gravel particles store small spheres composing them with Table 1(a) by
MAXNSCG (the maximum number of spheres composing a gravel particle) and sphere
numbers corresponding to each NSCG (number of spheres composing a gravel particle). Small
spheres store the gravel particle number composing themselves with Table 1(b).

2.3 Data flow in calculation of contact forces
The contact force and the motion of large numbers of particles in the numerical movable bed
have to be calculated with the high speed. Nishiura and Sakaguchi [3] have developed theefficient algorithms of calculations of contact forces and motions of spheres on shared- 224 -Percentage ( % )Slope of the initial movable bed: 1:2015 mmemory systems. We applied the PGS method
120mm 90mm
70mm
(Contact candidate Paring of Grid based Sorted
50mm40mm
particle index) and the FSC method (ForceSimulation referencing the Contact pair table) Before the experimentdeveloped by Nishiura and Sakaguchi to themotions of gravel particles of various sizes and1mshapes. The general method for detections of0.3 mspheres’ contacts stores spheres into the near t = 400 scells. It is, however, difficult to estimate
maximum numbers of spheres to be in a cell
Figure 3. Variation of the bed surface.
before the simulation in the case of spheres of100
various sizes (shown in Figure 2(a)). This makes
80
60
the determination of array size of cells for40
storing spheres difficult. The PGS method sorts
200
the sphere number by the cell number shown in
10100
1000
Fig 2 and Table 2. This makes it possible toParticle size ( mm )store all spheres in the cells by storing the Figure 4. Particle size distributions.
minimum sphere number and the maximum sphere number of spheres to be in the cell after
the sorting. Therefore, the PGS method is effective for contact detection of spheres of various
sizes in the numerical movable bed channel. Each sphere searches neighboring cells with
Table 2 and makes the list of contact candidate pairs (Table 3 (a),(b),(c)) and Table 4
connecting the sphere number with the number of the list of contact candidate pairs (Table3).
Each sphere refers to Table 1(b) and does not add overlapping spheres composing samegravel particles to the list of contact candidate pairs (Table 3). The FSC method calculates thecontact forces Fn of contacting two spheres by referring to the list of contact candidate pairs(Table 3). To calculate torques acting on gravel particles of various sizes and shapes, we store
contact points Xn with the Table 3. In the calculation of total contact forces Fg and torques Tg
acting on a gravel particle, the gravel particle call spheres composing the gravel particles byreferring to Table 1(a). The spheres refer to the contact pair list Table 3 by using Table 4 andcalculate torques Tn by using contact forces Fn and contact points Xn stored in Table 3 and thegravity center of the gravel particle. Then the contact forces Fn and torques Tn are added tothe forces Fg and torques Tg acting on the gravel particles. The routine summing up forces andtorques of gravel particles is conducted with the loop of the gravel particle number in Table
1(a) and is able to be parallelized on a shared memory system without data conflict.

3. NUMERICAL MOVABLE BED EXPERIMENTFigure 3 shows the numerical movable bed channel. 20 types of particles (5 in sizes (40 mm,50 mm, 70 mm, 90 mm, 120 mm) × 4 in shapes shown in Figure 1) were packed randomly inthe channel of 15m long so that particle size distributions in the numerical channel agreed
with the particle size distributions in Figure 4. The discharge of 0.5m3/s is given at theupstream end. The numerical experiment was lasted for 400 s. This numerical movable bed
experiment was carried out by 33,470 gravel particles, 301,231 spheres composing gravel
particles and 24,319,980 fluid calculation cells (∆x=∆y=∆z= 0.01 m).
4. RESULTS OF THE NUMERICAL MOVABLE BED EXPERIMENT
Figure 3 shows movable bed surfaces before and after the experiment. Large particles formedclusters and particle sortings occurred on the bed surface. Large particle clusters madeprotruding sections on the bed surface. The strong relationship was found between the sorting- 225 -D=120mm

D=90mm

D=70mm

D=50mm

D=40mmTotal

Sediment transport
rate (×10-3 m2/s )on the bed surface and the bed topography.
4.0Thewatersurface
variationof3.0super-critical flow was simulated stably by
2.0
the relatively simple estimation method.
1.0
The developed numerical movable bed
0.0
050
100
150
200250300
350
400
channel provides sound explanations fort(s)Figure 5. Temporal variation of sediment transport rate.
the motion of particles of various sizes andshapes and antidune formations. Figure 5 shows sediment transport rate of every particle size
calculated at the downstream end. Sediment transport rate of smaller particles were great in
the field experiments conducted in a gravel bed river (The Joganji river in Japan ) [4].
However, the sediment transport rates of larger particles are greater in the present numericalmovable bed channel. The difference seems to be caused by the following reasons: the first,smaller particles than 40 mm are not included in the numerical movable bed channel which
are transported easily and have a role of rollers for the motion of larger particles, the second,
the fluid calculation cell size (one fourth of the smallest particles size) is relatively large forsimulating fluid motions around smaller particles. The present numerical movable bed
experiment took about 3 months by the calculation using shared memory systems (64core:Intel Xeon E7-8837 (2.67GHz) 8core×8). It is too much computation time.
5. CONCLUSIONSThe present paper discussed results derived from the numerical movable bed channel
experiments with particles of various sizes and shapes. Following main conclusions areobtained.
1)The numerical movable bed channel is able to simulate reasonably the motion of particles
of various sizes and shapes, the sorting of particles at the bed surface and the motion of theflow on the bed waves.
2)The numerical channel with relatively large gravel particles showed that sediment transport
rates of larger particles were greater than smaller particles. This seems not to be seen in theordinary gravel bed rivers. The reason is because particles smaller than 40 mm are notincluded in the numerical channel and the fluid calculation cell size (one fourth of thesmallest particles size) is larger for simulating fluid motions around small particles. The
present numerical experiment took about three months. Large parallel computations on
distributed memory systems are indispensable for higher performance of the computation
with smaller particles and smaller fluid calculation cells.
References[1] T. Fukuda, S. Fukuoka and T. Uchida: Three-dimensional motions in numerical movable
bed channel with particles of various shapes and sizes, 12th International Symposium onRiver Sedimentation, Kyoto, Japan.
[2] S. Ushijima, A. Fukutani and O. Makino: Prediction method for movements withcollisions of arbitrary-shaped objects in 3D free-surface flows. JSCE Journal B 64(2),pp.128-138, 2008. in Japanese.
[3] D. Nishiura and H. Sakaguchi: Parallel-vector algorithms for particle simulations onshared-memory multiprocessors, J. Computational physics, No.230, pp.1923-1938. 2011.
[4] K. Soyama, Y. Ohkuma, Y. Hatanaka, F. Asano and S. Fukuoka: Study on observation
technique and assessment of bed load discharge, Advances in river engineering, Vol. 17,
pp.5-10, 2011, in Japanese.

- 226 -